operatorordering.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290
  1. """Functions for reordering operator expressions."""
  2. import warnings
  3. from sympy.core.add import Add
  4. from sympy.core.mul import Mul
  5. from sympy.core.numbers import Integer
  6. from sympy.core.power import Pow
  7. from sympy.physics.quantum import Commutator, AntiCommutator
  8. from sympy.physics.quantum.boson import BosonOp
  9. from sympy.physics.quantum.fermion import FermionOp
  10. __all__ = [
  11. 'normal_order',
  12. 'normal_ordered_form'
  13. ]
  14. def _expand_powers(factors):
  15. """
  16. Helper function for normal_ordered_form and normal_order: Expand a
  17. power expression to a multiplication expression so that that the
  18. expression can be handled by the normal ordering functions.
  19. """
  20. new_factors = []
  21. for factor in factors.args:
  22. if (isinstance(factor, Pow)
  23. and isinstance(factor.args[1], Integer)
  24. and factor.args[1] > 0):
  25. for n in range(factor.args[1]):
  26. new_factors.append(factor.args[0])
  27. else:
  28. new_factors.append(factor)
  29. return new_factors
  30. def _normal_ordered_form_factor(product, independent=False, recursive_limit=10,
  31. _recursive_depth=0):
  32. """
  33. Helper function for normal_ordered_form_factor: Write multiplication
  34. expression with bosonic or fermionic operators on normally ordered form,
  35. using the bosonic and fermionic commutation relations. The resulting
  36. operator expression is equivalent to the argument, but will in general be
  37. a sum of operator products instead of a simple product.
  38. """
  39. factors = _expand_powers(product)
  40. new_factors = []
  41. n = 0
  42. while n < len(factors) - 1:
  43. current, next = factors[n], factors[n + 1]
  44. if any(not isinstance(f, (FermionOp, BosonOp)) for f in (current, next)):
  45. new_factors.append(current)
  46. n += 1
  47. continue
  48. key_1 = (current.is_annihilation, str(current.name))
  49. key_2 = (next.is_annihilation, str(next.name))
  50. if key_1 <= key_2:
  51. new_factors.append(current)
  52. n += 1
  53. continue
  54. n += 2
  55. if current.is_annihilation and not next.is_annihilation:
  56. if isinstance(current, BosonOp) and isinstance(next, BosonOp):
  57. if current.args[0] != next.args[0]:
  58. if independent:
  59. c = 0
  60. else:
  61. c = Commutator(current, next)
  62. new_factors.append(next * current + c)
  63. else:
  64. new_factors.append(next * current + 1)
  65. elif isinstance(current, FermionOp) and isinstance(next, FermionOp):
  66. if current.args[0] != next.args[0]:
  67. if independent:
  68. c = 0
  69. else:
  70. c = AntiCommutator(current, next)
  71. new_factors.append(-next * current + c)
  72. else:
  73. new_factors.append(-next * current + 1)
  74. elif (current.is_annihilation == next.is_annihilation and
  75. isinstance(current, FermionOp) and isinstance(next, FermionOp)):
  76. new_factors.append(-next * current)
  77. else:
  78. new_factors.append(next * current)
  79. if n == len(factors) - 1:
  80. new_factors.append(factors[-1])
  81. if new_factors == factors:
  82. return product
  83. else:
  84. expr = Mul(*new_factors).expand()
  85. return normal_ordered_form(expr,
  86. recursive_limit=recursive_limit,
  87. _recursive_depth=_recursive_depth + 1,
  88. independent=independent)
  89. def _normal_ordered_form_terms(expr, independent=False, recursive_limit=10,
  90. _recursive_depth=0):
  91. """
  92. Helper function for normal_ordered_form: loop through each term in an
  93. addition expression and call _normal_ordered_form_factor to perform the
  94. factor to an normally ordered expression.
  95. """
  96. new_terms = []
  97. for term in expr.args:
  98. if isinstance(term, Mul):
  99. new_term = _normal_ordered_form_factor(
  100. term, recursive_limit=recursive_limit,
  101. _recursive_depth=_recursive_depth, independent=independent)
  102. new_terms.append(new_term)
  103. else:
  104. new_terms.append(term)
  105. return Add(*new_terms)
  106. def normal_ordered_form(expr, independent=False, recursive_limit=10,
  107. _recursive_depth=0):
  108. """Write an expression with bosonic or fermionic operators on normal
  109. ordered form, where each term is normally ordered. Note that this
  110. normal ordered form is equivalent to the original expression.
  111. Parameters
  112. ==========
  113. expr : expression
  114. The expression write on normal ordered form.
  115. independent : bool (default False)
  116. Whether to consider operator with different names as operating in
  117. different Hilbert spaces. If False, the (anti-)commutation is left
  118. explicit.
  119. recursive_limit : int (default 10)
  120. The number of allowed recursive applications of the function.
  121. Examples
  122. ========
  123. >>> from sympy.physics.quantum import Dagger
  124. >>> from sympy.physics.quantum.boson import BosonOp
  125. >>> from sympy.physics.quantum.operatorordering import normal_ordered_form
  126. >>> a = BosonOp("a")
  127. >>> normal_ordered_form(a * Dagger(a))
  128. 1 + Dagger(a)*a
  129. """
  130. if _recursive_depth > recursive_limit:
  131. warnings.warn("Too many recursions, aborting")
  132. return expr
  133. if isinstance(expr, Add):
  134. return _normal_ordered_form_terms(expr,
  135. recursive_limit=recursive_limit,
  136. _recursive_depth=_recursive_depth,
  137. independent=independent)
  138. elif isinstance(expr, Mul):
  139. return _normal_ordered_form_factor(expr,
  140. recursive_limit=recursive_limit,
  141. _recursive_depth=_recursive_depth,
  142. independent=independent)
  143. else:
  144. return expr
  145. def _normal_order_factor(product, recursive_limit=10, _recursive_depth=0):
  146. """
  147. Helper function for normal_order: Normal order a multiplication expression
  148. with bosonic or fermionic operators. In general the resulting operator
  149. expression will not be equivalent to original product.
  150. """
  151. factors = _expand_powers(product)
  152. n = 0
  153. new_factors = []
  154. while n < len(factors) - 1:
  155. if (isinstance(factors[n], BosonOp) and
  156. factors[n].is_annihilation):
  157. # boson
  158. if not isinstance(factors[n + 1], BosonOp):
  159. new_factors.append(factors[n])
  160. else:
  161. if factors[n + 1].is_annihilation:
  162. new_factors.append(factors[n])
  163. else:
  164. if factors[n].args[0] != factors[n + 1].args[0]:
  165. new_factors.append(factors[n + 1] * factors[n])
  166. else:
  167. new_factors.append(factors[n + 1] * factors[n])
  168. n += 1
  169. elif (isinstance(factors[n], FermionOp) and
  170. factors[n].is_annihilation):
  171. # fermion
  172. if not isinstance(factors[n + 1], FermionOp):
  173. new_factors.append(factors[n])
  174. else:
  175. if factors[n + 1].is_annihilation:
  176. new_factors.append(factors[n])
  177. else:
  178. if factors[n].args[0] != factors[n + 1].args[0]:
  179. new_factors.append(-factors[n + 1] * factors[n])
  180. else:
  181. new_factors.append(-factors[n + 1] * factors[n])
  182. n += 1
  183. else:
  184. new_factors.append(factors[n])
  185. n += 1
  186. if n == len(factors) - 1:
  187. new_factors.append(factors[-1])
  188. if new_factors == factors:
  189. return product
  190. else:
  191. expr = Mul(*new_factors).expand()
  192. return normal_order(expr,
  193. recursive_limit=recursive_limit,
  194. _recursive_depth=_recursive_depth + 1)
  195. def _normal_order_terms(expr, recursive_limit=10, _recursive_depth=0):
  196. """
  197. Helper function for normal_order: look through each term in an addition
  198. expression and call _normal_order_factor to perform the normal ordering
  199. on the factors.
  200. """
  201. new_terms = []
  202. for term in expr.args:
  203. if isinstance(term, Mul):
  204. new_term = _normal_order_factor(term,
  205. recursive_limit=recursive_limit,
  206. _recursive_depth=_recursive_depth)
  207. new_terms.append(new_term)
  208. else:
  209. new_terms.append(term)
  210. return Add(*new_terms)
  211. def normal_order(expr, recursive_limit=10, _recursive_depth=0):
  212. """Normal order an expression with bosonic or fermionic operators. Note
  213. that this normal order is not equivalent to the original expression, but
  214. the creation and annihilation operators in each term in expr is reordered
  215. so that the expression becomes normal ordered.
  216. Parameters
  217. ==========
  218. expr : expression
  219. The expression to normal order.
  220. recursive_limit : int (default 10)
  221. The number of allowed recursive applications of the function.
  222. Examples
  223. ========
  224. >>> from sympy.physics.quantum import Dagger
  225. >>> from sympy.physics.quantum.boson import BosonOp
  226. >>> from sympy.physics.quantum.operatorordering import normal_order
  227. >>> a = BosonOp("a")
  228. >>> normal_order(a * Dagger(a))
  229. Dagger(a)*a
  230. """
  231. if _recursive_depth > recursive_limit:
  232. warnings.warn("Too many recursions, aborting")
  233. return expr
  234. if isinstance(expr, Add):
  235. return _normal_order_terms(expr, recursive_limit=recursive_limit,
  236. _recursive_depth=_recursive_depth)
  237. elif isinstance(expr, Mul):
  238. return _normal_order_factor(expr, recursive_limit=recursive_limit,
  239. _recursive_depth=_recursive_depth)
  240. else:
  241. return expr