Navigation

Source code for sympy.physics.secondquant

"""Second quantization operators and states for bosons.This follow the formulation of Fetter and Welecka, "Quantum Theoryof Many-Particle Systems.""""fromcollectionsimportdefaultdictfromsympyimport(Add,Basic,cacheit,Dummy,Expr,Function,I,KroneckerDelta,Mul,Pow,S,sqrt,Symbol,sympify,Tuple,zeros)fromsympy.core.compatibilityimportreducefromsympy.printing.strimportStrPrinterfromsympy.physics.quantum.qexprimportsplit_commutative_partsfromsympy.core.compatibilityimportreducefromsympy.utilities.iterablesimporthas_dupsfromfunctoolsimportreduce__all__=['Dagger','KroneckerDelta','BosonicOperator','AnnihilateBoson','CreateBoson','AnnihilateFermion','CreateFermion','FockState','FockStateBra','FockStateKet','FockStateBosonKet','FockStateBosonBra','BBra','BKet','FBra','FKet','F','Fd','B','Bd','apply_operators','InnerProduct','BosonicBasis','VarBosonicBasis','FixedBosonicBasis','Commutator','matrix_rep','contraction','wicks','NO','evaluate_deltas','AntiSymmetricTensor','substitute_dummies','PermutationOperator','simplify_index_permutations',]classSecondQuantizationError(Exception):passclassAppliesOnlyToSymbolicIndex(SecondQuantizationError):passclassContractionAppliesOnlyToFermions(SecondQuantizationError):passclassViolationOfPauliPrinciple(SecondQuantizationError):passclassSubstitutionOfAmbigousOperatorFailed(SecondQuantizationError):passclassWicksTheoremDoesNotApply(SecondQuantizationError):pass

[docs]classFockState(Expr):""" Many particle Fock state with a sequence of occupation numbers. Anywhere you can have a FockState, you can also have S.Zero. All code must check for this! Base class to represent FockStates. """is_commutative=Falsedef__new__(cls,occupations):""" occupations is a list with two possible meanings: - For bosons it is a list of occupation numbers. Element i is the number of particles in state i. - For fermions it is a list of occupied orbits. Element 0 is the state that was occupied first, element i is the i'th occupied state. """occupations=list(map(sympify,occupations))obj=Basic.__new__(cls,Tuple(*occupations))returnobjdef__getitem__(self,i):i=int(i)returnself.args[0][i]def__repr__(self):return("FockState(%r)")%(self.args)def__str__(self):return"%s%r%s"%(self.lbracket,self._labels(),self.rbracket)def_labels(self):returnself.args[0]def__len__(self):returnlen(self.args[0])

classBosonState(FockState):""" Base class for FockStateBoson(Ket/Bra). """defup(self,i):""" Performs the action of a creation operator. Examples ======== >>> from sympy.physics.secondquant import BBra >>> b = BBra([1, 2]) >>> b FockStateBosonBra((1, 2)) >>> b.up(1) FockStateBosonBra((1, 3)) """i=int(i)new_occs=list(self.args[0])new_occs[i]=new_occs[i]+S.Onereturnself.__class__(new_occs)defdown(self,i):""" Performs the action of an annihilation operator. Examples ======== >>> from sympy.physics.secondquant import BBra >>> b = BBra([1, 2]) >>> b FockStateBosonBra((1, 2)) >>> b.down(1) FockStateBosonBra((1, 1)) """i=int(i)new_occs=list(self.args[0])ifnew_occs[i]==S.Zero:returnS.Zeroelse:new_occs[i]=new_occs[i]-S.Onereturnself.__class__(new_occs)classFermionState(FockState):""" Base class for FockStateFermion(Ket/Bra). """fermi_level=0def__new__(cls,occupations,fermi_level=0):occupations=list(map(sympify,occupations))iflen(occupations)>1:try:(occupations,sign)=_sort_anticommuting_fermions(occupations,key=hash)exceptViolationOfPauliPrinciple:returnS.Zeroelse:sign=0cls.fermi_level=fermi_levelifcls._count_holes(occupations)>fermi_level:returnS.Zeroifsign%2:returnS.NegativeOne*FockState.__new__(cls,occupations)else:returnFockState.__new__(cls,occupations)defup(self,i):""" Performs the action of a creation operator. If below fermi we try to remove a hole, if above fermi we try to create a particle. if general index p we return Kronecker(p,i)*self where i is a new symbol with restriction above or below. >>> from sympy import Symbol >>> from sympy.physics.secondquant import FKet >>> a = Symbol('a', above_fermi=True) >>> i = Symbol('i', below_fermi=True) >>> p = Symbol('p') >>> FKet([]).up(a) FockStateFermionKet((a,)) A creator acting on vacuum below fermi vanishes >>> FKet([]).up(i) 0 """present=iinself.args[0]ifself._only_above_fermi(i):ifpresent:returnS.Zeroelse:returnself._add_orbit(i)elifself._only_below_fermi(i):ifpresent:returnself._remove_orbit(i)else:returnS.Zeroelse:ifpresent:hole=Dummy("i",below_fermi=True)returnKroneckerDelta(i,hole)*self._remove_orbit(i)else:particle=Dummy("a",above_fermi=True)returnKroneckerDelta(i,particle)*self._add_orbit(i)defdown(self,i):""" Performs the action of an annihilation operator. If below fermi we try to create a hole, if above fermi we try to remove a particle. if general index p we return Kronecker(p,i)*self where i is a new symbol with restriction above or below. >>> from sympy import Symbol >>> from sympy.physics.secondquant import FKet >>> a = Symbol('a', above_fermi=True) >>> i = Symbol('i', below_fermi=True) >>> p = Symbol('p') An annihilator acting on vacuum above fermi vanishes >>> FKet([]).down(a) 0 Also below fermi, it vanishes, unless we specify a fermi level > 0 >>> FKet([]).down(i) 0 >>> FKet([],4).down(i) FockStateFermionKet((i,)) """present=iinself.args[0]ifself._only_above_fermi(i):ifpresent:returnself._remove_orbit(i)else:returnS.Zeroelifself._only_below_fermi(i):ifpresent:returnS.Zeroelse:returnself._add_orbit(i)else:ifpresent:hole=Dummy("i",below_fermi=True)returnKroneckerDelta(i,hole)*self._add_orbit(i)else:particle=Dummy("a",above_fermi=True)returnKroneckerDelta(i,particle)*self._remove_orbit(i)@classmethoddef_only_below_fermi(cls,i):""" Tests if given orbit is only below fermi surface. If nothing can be concluded we return a conservative False. """ifi.is_number:returni<=cls.fermi_levelifi.assumptions0.get('below_fermi'):returnTruereturnFalse@classmethoddef_only_above_fermi(cls,i):""" Tests if given orbit is only above fermi surface. If fermi level has not been set we return True. If nothing can be concluded we return a conservative False. """ifi.is_number:returni>cls.fermi_levelifi.assumptions0.get('above_fermi'):returnTruereturnnotcls.fermi_leveldef_remove_orbit(self,i):""" Removes particle/fills hole in orbit i. No input tests performed here. """new_occs=list(self.args[0])pos=new_occs.index(i)delnew_occs[pos]if(pos)%2:returnS.NegativeOne*self.__class__(new_occs,self.fermi_level)else:returnself.__class__(new_occs,self.fermi_level)def_add_orbit(self,i):""" Adds particle/creates hole in orbit i. No input tests performed here. """returnself.__class__((i,)+self.args[0],self.fermi_level)@classmethoddef_count_holes(cls,list):""" returns number of identified hole states in list. """returnlen([iforiinlistifcls._only_below_fermi(i)])def_negate_holes(self,list):returntuple([-iifi<=self.fermi_levelelseiforiinlist])def__repr__(self):ifself.fermi_level:return"FockStateKet(%r, fermi_level=%s)"%(self.args[0],self.fermi_level)else:return"FockStateKet(%r)"%(self.args[0],)def_labels(self):returnself._negate_holes(self.args[0])

[docs]classFockStateKet(FockState):""" Representation of a ket. """lbracket='|'rbracket='>'

[docs]classFockStateBra(FockState):""" Representation of a bra. """lbracket='<'rbracket='|'def__mul__(self,other):ifisinstance(other,FockStateKet):returnInnerProduct(self,other)else:returnExpr.__mul__(self,other)

classFockStateFermionKet(FermionState,FockStateKet):""" Many-particle Fock state with a sequence of occupied orbits. Each state can only have one particle, so we choose to store a list of occupied orbits rather than a tuple with occupation numbers (zeros and ones). states below fermi level are holes, and are represented by negative labels in the occupation list. For symbolic state labels, the fermi_level caps the number of allowed hole- states. Examples ======== >>> from sympy.physics.secondquant import FKet >>> FKet([1, 2]) #doctest: +SKIP FockStateFermionKet((1, 2)) """def_dagger_(self):returnFockStateFermionBra(*self.args)classFockStateFermionBra(FermionState,FockStateBra):""" See Also ======== FockStateFermionKet Examples ======== >>> from sympy.physics.secondquant import FBra >>> FBra([1, 2]) #doctest: +SKIP FockStateFermionBra((1, 2)) """def_dagger_(self):returnFockStateFermionKet(*self.args)BBra=FockStateBosonBraBKet=FockStateBosonKetFBra=FockStateFermionBraFKet=FockStateFermionKetdef_apply_Mul(m):""" Take a Mul instance with operators and apply them to states. This method applies all operators with integer state labels to the actual states. For symbolic state labels, nothing is done. When inner products of FockStates are encountered (like <a|b>), they are converted to instances of InnerProduct. This does not currently work on double inner products like, <a|b><c|d>. If the argument is not a Mul, it is simply returned as is. """ifnotisinstance(m,Mul):returnmc_part,nc_part=m.args_cnc()n_nc=len(nc_part)ifn_nc==0orn_nc==1:returnmelse:last=nc_part[-1]next_to_last=nc_part[-2]ifisinstance(last,FockStateKet):ifisinstance(next_to_last,SqOperator):ifnext_to_last.is_symbolic:returnmelse:result=next_to_last.apply_operator(last)ifresult==0:returnS.Zeroelse:return_apply_Mul(Mul(*(c_part+nc_part[:-2]+[result])))elifisinstance(next_to_last,Pow):ifisinstance(next_to_last.base,SqOperator)and \
next_to_last.exp.is_Integer:ifnext_to_last.base.is_symbolic:returnmelse:result=lastforiinrange(next_to_last.exp):result=next_to_last.base.apply_operator(result)ifresult==0:breakifresult==0:returnS.Zeroelse:return_apply_Mul(Mul(*(c_part+nc_part[:-2]+[result])))else:returnmelifisinstance(next_to_last,FockStateBra):result=InnerProduct(next_to_last,last)ifresult==0:returnS.Zeroelse:return_apply_Mul(Mul(*(c_part+nc_part[:-2]+[result])))else:returnmelse:returnm

[docs]defapply_operators(e):""" Take a sympy expression with operators and states and apply the operators. Examples ======== >>> from sympy.physics.secondquant import apply_operators >>> from sympy import sympify >>> apply_operators(sympify(3)+4) 7 """e=e.expand()muls=e.atoms(Mul)subs_list=[(m,_apply_Mul(m))forminiter(muls)]returne.subs(subs_list)

[docs]classInnerProduct(Basic):""" An unevaluated inner product between a bra and ket. Currently this class just reduces things to a product of Kronecker Deltas. In the future, we could introduce abstract states like ``|a>`` and ``|b>``, and leave the inner product unevaluated as ``<a|b>``. """is_commutative=Truedef__new__(cls,bra,ket):assertisinstance(bra,FockStateBra),'must be a bra'assertisinstance(ket,FockStateKet),'must be a key'returncls.eval(bra,ket)@classmethoddefeval(cls,bra,ket):result=S.Onefori,jinzip(bra.args[0],ket.args[0]):result*=KroneckerDelta(i,j)ifresult==0:breakreturnresult@property

[docs]defbra(self):"""Returns the bra part of the state"""returnself.args[0]

@property

[docs]defket(self):"""Returns the ket part of the state"""returnself.args[1]

[docs]classNO(Expr):""" This Object is used to represent normal ordering brackets. i.e. {abcd} sometimes written :abcd: Applying the function NO(arg) to an argument means that all operators in the argument will be assumed to anticommute, and have vanishing contractions. This allows an immediate reordering to canonical form upon object creation. >>> from sympy import symbols >>> from sympy.physics.secondquant import NO, F, Fd >>> p,q = symbols('p,q') >>> NO(Fd(p)*F(q)) NO(CreateFermion(p)*AnnihilateFermion(q)) >>> NO(F(q)*Fd(p)) -NO(CreateFermion(p)*AnnihilateFermion(q)) Note: If you want to generate a normal ordered equivalent of an expression, you should use the function wicks(). This class only indicates that all operators inside the brackets anticommute, and have vanishing contractions. Nothing more, nothing less. """nargs=1is_commutative=Falsedef__new__(cls,arg):""" Use anticommutation to get canonical form of operators. Employ associativity of normal ordered product: {ab{cd}} = {abcd} but note that {ab}{cd} /= {abcd}. We also employ distributivity: {ab + cd} = {ab} + {cd}. Canonical form also implies expand() {ab(c+d)} = {abc} + {abd}. """# {ab + cd} = {ab} + {cd}arg=sympify(arg)arg=arg.expand()ifarg.is_Add:returnAdd(*[cls(term)forterminarg.args])ifarg.is_Mul:# take coefficient outside of normal ordering bracketsc_part,seq=arg.args_cnc()ifc_part:coeff=Mul(*c_part)ifnotseq:returncoeffelse:coeff=S.One# {ab{cd}} = {abcd}newseq=[]foundit=Falseforfacinseq:ifisinstance(fac,NO):newseq.extend(fac.args)foundit=Trueelse:newseq.append(fac)iffoundit:returncoeff*cls(Mul(*newseq))# We assume that the user don't mix B and F operatorsifisinstance(seq[0],BosonicOperator):raiseNotImplementedErrortry:newseq,sign=_sort_anticommuting_fermions(seq)exceptViolationOfPauliPrinciple:returnS.Zeroifsign%2:return(S.NegativeOne*coeff)*cls(Mul(*newseq))elifsign:returncoeff*cls(Mul(*newseq))else:pass#since sign==0, no permutations was necessary# if we couldn't do anything with Mul object, we just# mark it as normal orderedifcoeff!=S.One:returncoeff*cls(Mul(*newseq))returnExpr.__new__(cls,Mul(*newseq))ifisinstance(arg,NO):returnarg# if object was not Mul or Add, normal ordering does not applyreturnarg@property

[docs]defhas_q_creators(self):""" Return 0 if the leftmost argument of the first argument is a not a q_creator, else 1 if it is above fermi or -1 if it is below fermi. Examples ======== >>> from sympy import symbols >>> from sympy.physics.secondquant import NO, F, Fd >>> a = symbols('a', above_fermi=True) >>> i = symbols('i', below_fermi=True) >>> NO(Fd(a)*Fd(i)).has_q_creators 1 >>> NO(F(i)*F(a)).has_q_creators -1 >>> NO(Fd(i)*F(a)).has_q_creators #doctest: +SKIP 0 """returnself.args[0].args[0].is_q_creator

@property

[docs]defhas_q_annihilators(self):""" Return 0 if the rightmost argument of the first argument is a not a q_annihilator, else 1 if it is above fermi or -1 if it is below fermi. Examples ======== >>> from sympy import symbols >>> from sympy.physics.secondquant import NO, F, Fd >>> a = symbols('a', above_fermi=True) >>> i = symbols('i', below_fermi=True) >>> NO(Fd(a)*Fd(i)).has_q_annihilators -1 >>> NO(F(i)*F(a)).has_q_annihilators 1 >>> NO(Fd(a)*F(i)).has_q_annihilators 0 """returnself.args[0].args[-1].is_q_annihilator

def_remove_brackets(self):""" Returns the sorted string without normal order brackets. The returned string have the property that no nonzero contractions exist. """# check if any creator is also an annihilatorsubslist=[]foriinself.iter_q_creators():ifself[i].is_q_annihilator:assume=self[i].state.assumptions0# only operators with a dummy index can be split in two termsifisinstance(self[i].state,Dummy):# create indices with fermi restrictionassume.pop("above_fermi",None)assume["below_fermi"]=Truebelow=Dummy('i',**assume)assume.pop("below_fermi",None)assume["above_fermi"]=Trueabove=Dummy('a',**assume)cls=type(self[i])split=(self[i].__new__(cls,below)*KroneckerDelta(below,self[i].state)+self[i].__new__(cls,above)*KroneckerDelta(above,self[i].state))subslist.append((self[i],split))else:raiseSubstitutionOfAmbigousOperatorFailed(self[i])ifsubslist:result=NO(self.subs(subslist))ifisinstance(result,Add):returnAdd(*[term.doit()forterminresult.args])else:returnself.args[0]def_expand_operators(self):""" Returns a sum of NO objects that contain no ambiguous q-operators. If an index q has range both above and below fermi, the operator F(q) is ambiguous in the sense that it can be both a q-creator and a q-annihilator. If q is dummy, it is assumed to be a summation variable and this method rewrites it into a sum of NO terms with unambiguous operators: {Fd(p)*F(q)} = {Fd(a)*F(b)} + {Fd(a)*F(i)} + {Fd(j)*F(b)} -{F(i)*Fd(j)} where a,b are above and i,j are below fermi level. """returnNO(self._remove_brackets)def__getitem__(self,i):ifisinstance(i,slice):indices=i.indices(len(self))return[self.args[0].args[i]foriinrange(*indices)]else:returnself.args[0].args[i]def__len__(self):returnlen(self.args[0].args)

def_sqkey(sq_operator):"""Generates key for canonical sorting of SQ operators."""returnsq_operator._sortkey()def_sort_anticommuting_fermions(string1,key=_sqkey):"""Sort fermionic operators to canonical order, assuming all pairs anticommute. Uses a bidirectional bubble sort. Items in string1 are not referenced so in principle they may be any comparable objects. The sorting depends on the operators '>' and '=='. If the Pauli principle is violated, an exception is raised. Returns ======= tuple (sorted_str, sign) sorted_str: list containing the sorted operators sign: int telling how many times the sign should be changed (if sign==0 the string was already sorted) """verified=Falsesign=0rng=list(range(len(string1)-1))rev=list(range(len(string1)-3,-1,-1))keys=list(map(key,string1))key_val=dict(list(zip(keys,string1)))whilenotverified:verified=Trueforiinrng:left=keys[i]right=keys[i+1]ifleft==right:raiseViolationOfPauliPrinciple([left,right])ifleft>right:verified=Falsekeys[i:i+2]=[right,left]sign=sign+1ifverified:breakforiinrev:left=keys[i]right=keys[i+1]ifleft==right:raiseViolationOfPauliPrinciple([left,right])ifleft>right:verified=Falsekeys[i:i+2]=[right,left]sign=sign+1string1=[key_val[k]forkinkeys]return(string1,sign)

[docs]defevaluate_deltas(e):""" We evaluate KroneckerDelta symbols in the expression assuming Einstein summation. If one index is repeated it is summed over and in effect substituted with the other one. If both indices are repeated we substitute according to what is the preferred index. this is determined by KroneckerDelta.preferred_index and KroneckerDelta.killable_index. In case there are no possible substitutions or if a substitution would imply a loss of information, nothing is done. In case an index appears in more than one KroneckerDelta, the resulting substitution depends on the order of the factors. Since the ordering is platform dependent, the literal expression resulting from this function may be hard to predict. Examples ======== We assume the following: >>> from sympy import symbols, Function, Dummy, KroneckerDelta >>> from sympy.physics.secondquant import evaluate_deltas >>> i,j = symbols('i j', below_fermi=True, cls=Dummy) >>> a,b = symbols('a b', above_fermi=True, cls=Dummy) >>> p,q = symbols('p q', cls=Dummy) >>> f = Function('f') >>> t = Function('t') The order of preference for these indices according to KroneckerDelta is (a, b, i, j, p, q). Trivial cases: >>> evaluate_deltas(KroneckerDelta(i,j)*f(i)) # d_ij f(i) -> f(j) f(_j) >>> evaluate_deltas(KroneckerDelta(i,j)*f(j)) # d_ij f(j) -> f(i) f(_i) >>> evaluate_deltas(KroneckerDelta(i,p)*f(p)) # d_ip f(p) -> f(i) f(_i) >>> evaluate_deltas(KroneckerDelta(q,p)*f(p)) # d_qp f(p) -> f(q) f(_q) >>> evaluate_deltas(KroneckerDelta(q,p)*f(q)) # d_qp f(q) -> f(p) f(_p) More interesting cases: >>> evaluate_deltas(KroneckerDelta(i,p)*t(a,i)*f(p,q)) f(_i, _q)*t(_a, _i) >>> evaluate_deltas(KroneckerDelta(a,p)*t(a,i)*f(p,q)) f(_a, _q)*t(_a, _i) >>> evaluate_deltas(KroneckerDelta(p,q)*f(p,q)) f(_p, _p) Finally, here are some cases where nothing is done, because that would imply a loss of information: >>> evaluate_deltas(KroneckerDelta(i,p)*f(q)) f(_q)*KroneckerDelta(_i, _p) >>> evaluate_deltas(KroneckerDelta(i,p)*f(i)) f(_i)*KroneckerDelta(_i, _p) """# We treat Deltas only in mul objects# for general function objects we don't evaluate KroneckerDeltas in arguments,# but here we hard code exceptions to this ruleaccepted_functions=(Add,)ifisinstance(e,accepted_functions):returne.func(*[evaluate_deltas(arg)forargine.args])elifisinstance(e,Mul):# find all occurences of delta function and count each index present in# expression.deltas=[]indices={}foriine.args:forsini.atoms():ifsinindices:indices[s]+=1else:indices[s]=0# geek counting simplifies logic belowifisinstance(i,KroneckerDelta):deltas.append(i)fordindeltas:# If we do something, and there are more deltas, we should recurse# to treat the resulting expression properlyifindices[d.killable_index]:e=e.subs(d.killable_index,d.preferred_index)iflen(deltas)>1:returnevaluate_deltas(e)elifindices[d.preferred_index]andd.indices_contain_equal_information:e=e.subs(d.preferred_index,d.killable_index)iflen(deltas)>1:returnevaluate_deltas(e)else:passreturne# nothing to do, maybe we hit a Symbol or a numberelse:returne

[docs]defsubstitute_dummies(expr,new_indices=False,pretty_indices={}):""" Collect terms by substitution of dummy variables. This routine allows simplification of Add expressions containing terms which differ only due to dummy variables. The idea is to substitute all dummy variables consistently depending on the structure of the term. For each term, we obtain a sequence of all dummy variables, where the order is determined by the index range, what factors the index belongs to and its position in each factor. See _get_ordered_dummies() for more inforation about the sorting of dummies. The index sequence is then substituted consistently in each term. Examples ======== >>> from sympy import symbols, Function, Dummy >>> from sympy.physics.secondquant import substitute_dummies >>> a,b,c,d = symbols('a b c d', above_fermi=True, cls=Dummy) >>> i,j = symbols('i j', below_fermi=True, cls=Dummy) >>> f = Function('f') >>> expr = f(a,b) + f(c,d); expr f(_a, _b) + f(_c, _d) Since a, b, c and d are equivalent summation indices, the expression can be simplified to a single term (for which the dummy indices are still summed over) >>> substitute_dummies(expr) 2*f(_a, _b) Controlling output: By default the dummy symbols that are already present in the expression will be reused in a different permuation. However, if new_indices=True, new dummies will be generated and inserted. The keyword 'pretty_indices' can be used to control this generation of new symbols. By default the new dummies will be generated on the form i_1, i_2, a_1, etc. If you supply a dictionary with key:value pairs in the form: { index_group: string_of_letters } The letters will be used as labels for the new dummy symbols. The index_groups must be one of 'above', 'below' or 'general'. >>> expr = f(a,b,i,j) >>> my_dummies = { 'above':'st', 'below':'uv' } >>> substitute_dummies(expr, new_indices=True, pretty_indices=my_dummies) f(_s, _t, _u, _v) If we run out of letters, or if there is no keyword for some index_group the default dummy generator will be used as a fallback: >>> p,q = symbols('p q', cls=Dummy) # general indices >>> expr = f(p,q) >>> substitute_dummies(expr, new_indices=True, pretty_indices=my_dummies) f(_p_0, _p_1) """# setup the replacing dummiesifnew_indices:letters_above=pretty_indices.get('above',"")letters_below=pretty_indices.get('below',"")letters_general=pretty_indices.get('general',"")len_above=len(letters_above)len_below=len(letters_below)len_general=len(letters_general)def_i(number):try:returnletters_below[number]exceptIndexError:return'i_'+str(number-len_below)def_a(number):try:returnletters_above[number]exceptIndexError:return'a_'+str(number-len_above)def_p(number):try:returnletters_general[number]exceptIndexError:return'p_'+str(number-len_general)aboves=[]belows=[]generals=[]dummies=expr.atoms(Dummy)ifnotnew_indices:dummies=sorted(dummies)# generate lists with the dummies we will inserta=i=p=0fordindummies:assum=d.assumptions0ifassum.get("above_fermi"):ifnew_indices:sym=_a(a);a+=1l1=aboveselifassum.get("below_fermi"):ifnew_indices:sym=_i(i);i+=1l1=belowselse:ifnew_indices:sym=_p(p);p+=1l1=generalsifnew_indices:l1.append(Dummy(sym,**assum))else:l1.append(d)expr=expr.expand()terms=Add.make_args(expr)new_terms=[]forterminterms:i=iter(belows)a=iter(aboves)p=iter(generals)ordered=_get_ordered_dummies(term)subsdict={}fordinordered:ifd.assumptions0.get('below_fermi'):subsdict[d]=next(i)elifd.assumptions0.get('above_fermi'):subsdict[d]=next(a)else:subsdict[d]=next(p)subslist=[]final_subs=[]fork,vinsubsdict.items():ifk==v:continueifvinsubsdict:# We check if the sequence of substitutions end quickly. In# that case, we can avoid temporary symbols if we ensure the# correct substitution order.ifsubsdict[v]insubsdict:# (x, y) -> (y, x), we need a temporary variablex=Dummy('x')subslist.append((k,x))final_subs.append((x,v))else:# (x, y) -> (y, a), x->y must be done last# but before temporary variables are resolvedfinal_subs.insert(0,(k,v))else:subslist.append((k,v))subslist.extend(final_subs)new_terms.append(term.subs(subslist))returnAdd(*new_terms)

classKeyPrinter(StrPrinter):"""Printer for which only equal objects are equal in print"""def_print_Dummy(self,expr):return"(%s_%i)"%(expr.name,expr.dummy_index)def__kprint(expr):p=KeyPrinter()returnp.doprint(expr)def_get_ordered_dummies(mul,verbose=False):"""Returns all dummies in the mul sorted in canonical order The purpose of the canonical ordering is that dummies can be substituted consistently across terms with the result that equivalent terms can be simplified. It is not possible to determine if two terms are equivalent based solely on the dummy order. However, a consistent substitution guided by the ordered dummies should lead to trivially (non-)equivalent terms, thereby revealing the equivalence. This also means that if two terms have identical sequences of dummies, the (non-)equivalence should already be apparent. Strategy -------- The canoncial order is given by an arbitrary sorting rule. A sort key is determined for each dummy as a tuple that depends on all factors where the index is present. The dummies are thereby sorted according to the contraction structure of the term, instead of sorting based solely on the dummy symbol itself. After all dummies in the term has been assigned a key, we check for identical keys, i.e. unorderable dummies. If any are found, we call a specialized method, _determine_ambiguous(), that will determine a unique order based on recursive calls to _get_ordered_dummies(). Key description --------------- A high level description of the sort key: 1. Range of the dummy index 2. Relation to external (non-dummy) indices 3. Position of the index in the first factor 4. Position of the index in the second factor The sort key is a tuple with the following components: 1. A single character indicating the range of the dummy (above, below or general.) 2. A list of strings with fully masked string representations of all factors where the dummy is present. By masked, we mean that dummies are represented by a symbol to indicate either below fermi, above or general. No other information is displayed about the dummies at this point. The list is sorted stringwise. 3. An integer number indicating the position of the index, in the first factor as sorted in 2. 4. An integer number indicating the position of the index, in the second factor as sorted in 2. If a factor is either of type AntiSymmetricTensor or SqOperator, the index position in items 3 and 4 is indicated as 'upper' or 'lower' only. (Creation operators are considered upper and annihilation operators lower.) If the masked factors are identical, the two factors cannot be ordered unambiguously in item 2. In this case, items 3, 4 are left out. If several indices are contracted between the unorderable factors, it will be handled by _determine_ambiguous() """# setup dicts to avoid repeated calculations in key()args=Mul.make_args(mul)fac_dum=dict([(fac,fac.atoms(Dummy))forfacinargs])fac_repr=dict([(fac,__kprint(fac))forfacinargs])all_dums=reduce(set.union,list(fac_dum.values()),set())mask={}fordinall_dums:ifd.assumptions0.get('below_fermi'):mask[d]='0'elifd.assumptions0.get('above_fermi'):mask[d]='1'else:mask[d]='2'dum_repr=dict([(d,__kprint(d))fordinall_dums])def_key(d):dumstruct=[facforfacinfac_dumifdinfac_dum[fac]]other_dums=reduce(set.union,[fac_dum[fac]forfacindumstruct],set())fac=dumstruct[-1]ifother_dumsisfac_dum[fac]:other_dums=fac_dum[fac].copy()other_dums.remove(d)masked_facs=[fac_repr[fac]forfacindumstruct]ford2inother_dums:masked_facs=[fac.replace(dum_repr[d2],mask[d2])forfacinmasked_facs]all_masked=[fac.replace(dum_repr[d],mask[d])forfacinmasked_facs]masked_facs=dict(list(zip(dumstruct,masked_facs)))# dummies for which the ordering cannot be determinedifhas_dups(all_masked):all_masked.sort()returnmask[d],tuple(all_masked)# positions are ambiguous# sort factors according to fully masked stringskeydict=dict(list(zip(dumstruct,all_masked)))dumstruct.sort(key=lambdax:keydict[x])all_masked.sort()pos_val=[]forfacindumstruct:ifisinstance(fac,AntiSymmetricTensor):ifdinfac.upper:pos_val.append('u')ifdinfac.lower:pos_val.append('l')elifisinstance(fac,Creator):pos_val.append('u')elifisinstance(fac,Annihilator):pos_val.append('l')elifisinstance(fac,NO):ops=[opforopinfacifop.has(d)]foropinops:ifisinstance(op,Creator):pos_val.append('u')else:pos_val.append('l')else:# fallback to position in string representationfacpos=-1while1:facpos=masked_facs[fac].find(dum_repr[d],facpos+1)iffacpos==-1:breakpos_val.append(facpos)return(mask[d],tuple(all_masked),pos_val[0],pos_val[-1])dumkey=dict(list(zip(all_dums,list(map(_key,all_dums)))))result=sorted(all_dums,key=lambdax:dumkey[x])ifhas_dups(iter(dumkey.values())):# We have ambiguitiesunordered=defaultdict(set)ford,kindumkey.items():unordered[k].add(d)forkin[kforkinunorderediflen(unordered[k])<2]:delunordered[k]unordered=[unordered[k]forkinsorted(unordered)]result=_determine_ambiguous(mul,result,unordered)returnresultdef_determine_ambiguous(term,ordered,ambiguous_groups):# We encountered a term for which the dummy substitution is ambiguous.# This happens for terms with 2 or more contractions between factors that# cannot be uniquely ordered independent of summation indices. For# example:## Sum(p, q) v^{p, .}_{q, .}v^{q, .}_{p, .}## Assuming that the indices represented by . are dummies with the# same range, the factors cannot be ordered, and there is no# way to determine a consistent ordering of p and q.## The strategy employed here, is to relabel all unambiguous dummies with# non-dummy symbols and call _get_ordered_dummies again. This procedure is# applied to the entire term so there is a possibility that# _determine_ambiguous() is called again from a deeper recursion level.# break recursion if there are no ordered dummiesall_ambiguous=set()fordummiesinambiguous_groups:all_ambiguous|=dummiesall_ordered=set(ordered)-all_ambiguousifnotall_ordered:# FIXME: If we arrive here, there are no ordered dummies. A method to# handle this needs to be implemented. In order to return something# useful nevertheless, we choose arbitrarily the first dummy and# determine the rest from this one. This method is dependent on the# actual dummy labels which violates an assumption for the canonization# procedure. A better implementation is needed.group=[dfordinorderedifdinambiguous_groups[0]]d=group[0]all_ordered.add(d)ambiguous_groups[0].remove(d)stored_counter=_symbol_factory._countersubslist=[]fordin[dfordinorderedifdinall_ordered]:nondum=_symbol_factory._next()subslist.append((d,nondum))newterm=term.subs(subslist)neworder=_get_ordered_dummies(newterm)_symbol_factory._set_counter(stored_counter)# update ordered list with new informationforgroupinambiguous_groups:ordered_group=[dfordinneworderifdingroup]ordered_group.reverse()result=[]fordinordered:ifdingroup:result.append(ordered_group.pop())else:result.append(d)ordered=resultreturnorderedclass_SymbolFactory(object):def__init__(self,label):self._counterVar=0self._label=labeldef_set_counter(self,value):""" Sets counter to value. """self._counterVar=value@propertydef_counter(self):""" What counter is currently at. """returnself._counterVardef_next(self):""" Generates the next symbols and increments counter by 1. """s=Symbol("%s%i"%(self._label,self._counterVar))self._counterVar+=1returns_symbol_factory=_SymbolFactory('_]"]_')# most certainly a unique label@cacheitdef_get_contractions(string1,keep_only_fully_contracted=False):""" Returns Add-object with contracted terms. Uses recursion to find all contractions. -- Internal helper function -- Will find nonzero contractions in string1 between indices given in leftrange and rightrange. """# Should we store current level of contraction?ifkeep_only_fully_contractedandstring1:result=[]else:result=[NO(Mul(*string1))]foriinrange(len(string1)-1):forjinrange(i+1,len(string1)):c=contraction(string1[i],string1[j])ifc:# print "found contraction",csign=(j-i+1)%2ifsign:coeff=S.NegativeOne*celse:coeff=c## Call next level of recursion# ============================## We now need to find more contractions among operators## oplist = string1[:i]+ string1[i+1:j] + string1[j+1:]## To prevent overcounting, we don't allow contractions# we have already encountered. i.e. contractions between# string1[:i] <---> string1[i+1:j]# and string1[:i] <---> string1[j+1:].## This leaves the case:oplist=string1[i+1:j]+string1[j+1:]ifoplist:result.append(coeff*NO(Mul(*string1[:i])*_get_contractions(oplist,keep_only_fully_contracted=keep_only_fully_contracted)))else:result.append(coeff*NO(Mul(*string1[:i])))ifkeep_only_fully_contracted:break# next iteration over i leaves leftmost operator string1[0] uncontractedreturnAdd(*result)# @cacheit

[docs]defwicks(e,**kw_args):""" Returns the normal ordered equivalent of an expression using Wicks Theorem. Examples ======== >>> from sympy import symbols, Function, Dummy >>> from sympy.physics.secondquant import wicks, F, Fd, NO >>> p,q,r = symbols('p,q,r') >>> wicks(Fd(p)*F(q)) # doctest: +SKIP d(p, q)*d(q, _i) + NO(CreateFermion(p)*AnnihilateFermion(q)) By default, the expression is expanded: >>> wicks(F(p)*(F(q)+F(r))) # doctest: +SKIP NO(AnnihilateFermion(p)*AnnihilateFermion(q)) + NO(AnnihilateFermion(p)*AnnihilateFermion(r)) With the keyword 'keep_only_fully_contracted=True', only fully contracted terms are returned. By request, the result can be simplified in the following order: -- KroneckerDelta functions are evaluated -- Dummy variables are substituted consistently across terms >>> p, q, r = symbols('p q r', cls=Dummy) >>> wicks(Fd(p)*(F(q)+F(r)), keep_only_fully_contracted=True) # doctest: +SKIP KroneckerDelta(_i, _q)*KroneckerDelta(_p, _q) + KroneckerDelta(_i, _r)*KroneckerDelta(_p, _r) """ifnote:returnS.Zeroopts={'simplify_kronecker_deltas':False,'expand':True,'simplify_dummies':False,'keep_only_fully_contracted':False}opts.update(kw_args)# check if we are already normally orderedifisinstance(e,NO):ifopts['keep_only_fully_contracted']:returnS.Zeroelse:returneelifisinstance(e,FermionicOperator):ifopts['keep_only_fully_contracted']:returnS.Zeroelse:returne# break up any NO-objects, and evaluate commutatorse=e.doit(wicks=True)# make sure we have only one term to considere=e.expand()ifisinstance(e,Add):ifopts['simplify_dummies']:returnsubstitute_dummies(Add(*[wicks(term,**kw_args)fortermine.args]))else:returnAdd(*[wicks(term,**kw_args)fortermine.args])# For Mul-objects we can actually do somethingifisinstance(e,Mul):# we dont want to mess around with commuting part of Mul# so we factorize it out before starting recursionc_part=[]string1=[]forfactorine.args:iffactor.is_commutative:c_part.append(factor)else:string1.append(factor)n=len(string1)# catch trivial casesifn==0:result=eelifn==1:ifopts['keep_only_fully_contracted']:returnS.Zeroelse:result=eelse:# non-trivialifisinstance(string1[0],BosonicOperator):raiseNotImplementedErrorstring1=tuple(string1)# recursion over higher order contractionsresult=_get_contractions(string1,keep_only_fully_contracted=opts['keep_only_fully_contracted'])result=Mul(*c_part)*resultifopts['expand']:result=result.expand()ifopts['simplify_kronecker_deltas']:result=evaluate_deltas(result)returnresult# there was nothing to doreturne