Skip to content

Commit

Permalink
- added sqOrder[] function defining canonical order of indices in SeQ…
Browse files Browse the repository at this point in the history 
…uant (compare by spin first, then by space, then by index) and use this function in orderedForm[]

- spintrace can be handle composite particle spaces and expressions involving SQS
  • Loading branch information
@evaleev
evaleev committed Jul 6, 2021
1 parent 0ffb6a1 commit 26b4834
Showing 1 changed file with 45 additions and 29 deletions.
74 changes: 45 additions & 29 deletions SeQuant.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@
(* this converts spin tag into m_s *)
Ms[particleSpin[A]]=1/2;
Ms[particleSpin[B]]=-1/2;
Ms[particleSpin[none]]=0;

stringAppendSpin[str_,s_particleSpin]:=If[s===particleSpin[none],str, "\!\(\*SubscriptBox[\("<>str<>"\), \("<>If[s===particleSpin[A],"\[Alpha]","\[Beta]"]<>"\)]\)"];

Expand Down Expand Up @@ -180,6 +181,10 @@
(* return particleSpace of particleIndex *)
indexSpace[a_particleIndex] :=
Cases[a,_particleSpace][[1]];

(* return spin-free particleSpace of particleIndex *)
indexSpaceSpinFree[a_particleIndex] :=
DeleteCases[indexSpace[a],_particleSpin];

(* return the particleType of a particleIndex *)
indexParticle[a_particleIndex] :=
Expand Down Expand Up @@ -218,14 +223,16 @@
(a[[1]]===b[[1]])&&(a[[2]]===b[[2]]);
Protect[Equal];

(* canonical ordering of particleIndex objects is defined as follows: order by particleSpace first, then alphabetically by the symbol *)
Unprotect[OrderedQ];
OrderedQ[A_particleSpace,B_particleSpace] :=
If[ OrderedQ[A[[2]],B[[2]]],
OrderedQ[A[[1]],B[[1]]],
True
];
Protect[OrderedQ];
(* canonical ordering of particleIndex objects is defined as follows: order by spin, then the rest of particleSpace, then alphabetically by the symbol.
particleSpace objects ordered by *)
sqOrder[A_particleIndex,B_particleIndex] :=
If[ Ms[indexSpin[A]]!=Ms[indexSpin[B]],
If[Ms[indexSpin[A]]>Ms[indexSpin[B]],1,-1],
If[ Order[indexSpaceSpinFree[A],indexSpaceSpinFree[B]]==0,
Order[A[[1]],B[[1]]],
Order[indexSpaceSpinFree[A],indexSpaceSpinFree[B]]
]];
sqOrder[a_,b_]:=Order[a,b];

(* this functions throws out all multiple occurences of the same index from inds *)
uniqueIndexList[inds_List] :=
Expand Down Expand Up @@ -801,7 +808,7 @@ CR is the head for a general result of a contraction (unless it's a zero)



(* ::Subsection:: *)
(* ::Subsection::Closed:: *)
(* Contract SQS *)


Expand Down Expand Up @@ -2058,11 +2065,11 @@ The canonical order of indices of each type (that belong to the same space) is a
permfac = 1;
creInds = Cases[str,a_particleIndex/;a[[3]]==indexType[cre]];
permfac*=Signature[creInds];
creInds = Sort[creInds];
creInds = Sort[creInds,sqOrder];
annInds = Cases[str,a_particleIndex/;a[[3]]==indexType[ann]];
annInds = Reverse[annInds];
permfac*=Signature[annInds];
annInds = Sort[annInds];
annInds = Sort[annInds,sqOrder];
result = permfac*FlattenAt[SQS[FlattenAt[{creInds,Reverse[annInds]},{{1},{2}}]],{1}];
Return[result];
];
Expand All @@ -2087,7 +2094,7 @@ The canonical order of indices of each type (that belong to the same space) is a
result = tagIndices[oper,intInds];
braInds = Cases[result,a_particleIndex/;a[[3]]==indexType[bra]];
(* canonically reorder bra indices first *)
braIndsOrd = Ordering[braInds];
braIndsOrd = Ordering[braInds,Length[braInds],sqOrder];
braInds = braInds[[braIndsOrd]];
If[ symfac==-1,
permfac *= Signature[braIndsOrd];
Expand All @@ -2100,7 +2107,7 @@ degenerately ordered (e.g. there are duplicates in the lists so that
multiple permutations of the original list re ordered *)
If[ symfac==0,
ketInds=ketInds[[braIndsOrd]],
( ketIndsOrd = Ordering[ketInds];
( ketIndsOrd = Ordering[ketInds,Length[ketInds],sqOrder];
If[ symfac==-1,
permfac*=Signature[ketIndsOrd]
];
Expand Down Expand Up @@ -2332,7 +2339,7 @@ degenerately ordered (e.g. there are duplicates in the lists so that


(*
Remove diconnected terms. Connectedness is determined with respect to the Hamiltonian
Remove disconnected terms. Connectedness is determined with respect to the Hamiltonian
*)

removeDisconnectedTerms[expr_,externalIndices_List:{},hamiltonianOpers_List:defaultHamiltonianOpers] :=
Expand Down Expand Up @@ -2585,7 +2592,7 @@ degenerately ordered (e.g. there are duplicates in the lists so that
Protect[MatchQ];


(* ::Section::Closed:: *)
(* ::Section:: *)
(* Spin*)


Expand Down Expand Up @@ -2620,7 +2627,7 @@ degenerately ordered (e.g. there are duplicates in the lists so that
Do[ (* for every combination of spins *)
(* make a list of index replacements to apply (groupings are omitted) *)
repls=
Table[particleIndex[indexSymbol[indGroups[[g,i]]],particleSpace[space_],rest_]->particleIndex[indexSymbol[indGroups[[g,i]]],particleSpace[space,particleSpin[mstuples[[t,g]]]],rest],{g,Length[indGroups]},{i,Length[indGroups[[g]]]}];
Table[particleIndex[indexSymbol[indGroups[[g,i]]],particleSpace[space__],rest__]->particleIndex[indexSymbol[indGroups[[g,i]]],particleSpace[space,particleSpin[mstuples[[t,g]]]],rest],{g,Length[indGroups]},{i,Length[indGroups[[g]]]}];
repls=Flatten[repls];
If[ SeQuantDebugLevel>=3,Print["in spintrace: t=",t," repls=",repls//TraditionalForm," expr/.repls=",TraditionalForm[expr/.repls]];
];
Expand All @@ -2633,26 +2640,35 @@ degenerately ordered (e.g. there are duplicates in the lists so that

(* convert to nonsymmetric n-body operators, e.g. Overscript[g, _] \[Rule] g *)
(* 2-body *)
result=result/.SQM[OHead[label_String,indexSymm[-1]],particleIndex[bra1__],particleIndex[bra2__],particleIndex[ket1__],particleIndex[ket2__]]->rhoIndex[particleIndex[bra1],particleIndex[bra2]]*rhoIndex[particleIndex[ket1],particleIndex[ket2]]*(SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[ket1],particleIndex[ket2]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[ket2],particleIndex[ket1]]);
result=result/.SQM[OHead[label_String?(RefLink[StringContainsQ,paclet:ref/StringContainsQ][#,StartOfLine~~"Overscript["~~__] &),indexSymm[-1]],particleIndex[bra1__],particleIndex[bra2__],particleIndex[ket1__],particleIndex[ket2__]]->rhoIndex[particleIndex[bra1],particleIndex[bra2]]*rhoIndex[particleIndex[ket1],particleIndex[ket2]]*(SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[ket1],particleIndex[ket2]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[ket2],particleIndex[ket1]]);
(* 3-body *)
result=result/.SQM[OHead[label_String,indexSymm[-1]],particleIndex[bra1__],particleIndex[bra2__],particleIndex[bra3__],particleIndex[ket1__],particleIndex[ket2__],particleIndex[ket3__]]->rhoIndex[particleIndex[bra1],particleIndex[bra2]]*rhoIndex[particleIndex[bra1],particleIndex[bra3]]*rhoIndex[particleIndex[bra2],particleIndex[bra3]]*rhoIndex[particleIndex[ket1],particleIndex[ket2]]**rhoIndex[particleIndex[ket1],particleIndex[ket3]]**rhoIndex[particleIndex[ket2],particleIndex[ket3]](SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket1],particleIndex[ket2],particleIndex[ket3]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket2],particleIndex[ket1],particleIndex[ket3]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket1],particleIndex[ket3],particleIndex[ket2]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket3],particleIndex[ket2],particleIndex[ket1]]+SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket2],particleIndex[ket3],particleIndex[ket1]]+SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket3],particleIndex[ket1],particleIndex[ket2]]);
(* eliminate rho's if EPV terms will be avoided implcitly, but summing multiple terms *)
result=result/.SQM[OHead[label_String?(RefLink[StringContainsQ,paclet:ref/StringContainsQ][#,StartOfLine~~"Overscript["~~__] &),indexSymm[-1]],particleIndex[bra1__],particleIndex[bra2__],particleIndex[bra3__],particleIndex[ket1__],particleIndex[ket2__],particleIndex[ket3__]]->rhoIndex[particleIndex[bra1],particleIndex[bra2]]*rhoIndex[particleIndex[bra1],particleIndex[bra3]]*rhoIndex[particleIndex[bra2],particleIndex[bra3]]*rhoIndex[particleIndex[ket1],particleIndex[ket2]]**rhoIndex[particleIndex[ket1],particleIndex[ket3]]**rhoIndex[particleIndex[ket2],particleIndex[ket3]](SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket1],particleIndex[ket2],particleIndex[ket3]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket2],particleIndex[ket1],particleIndex[ket3]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket1],particleIndex[ket3],particleIndex[ket2]]-SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket3],particleIndex[ket2],particleIndex[ket1]]+SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket2],particleIndex[ket3],particleIndex[ket1]]+SQM[OHead[antisymm2nonsymmSQMlabel[label],indexSymm[0]],particleIndex[bra1],particleIndex[bra2],particleIndex[bra3],particleIndex[ket3],particleIndex[ket1],particleIndex[ket2]]);
(* eliminate rho's if EPV terms will be avoided implicitly, but summing multiple terms *)
result=If[explicitEPV,result,(result/.rhoIndex[__]->1)];
result=Expand[result];
(* canonicalize to make spin tracing easier*)
result=Map[maporder,result,2];
If[ SeQuantDebugLevel>=2,Print["in spintrace: after expanding antisymmetric operators = ",TraditionalForm[result]];
];

(* zero out integrals that don't conserve spin quantum numbers *)
(* 1-body *)
result=result/.SQM[OHead[op__],particleIndex[bra1_String,particleSpace[bra1space_,bra1spin_],bra1rest_],particleIndex[ket1_String,particleSpace[ket1space_,ket1spin_],ket1rest_]]->If[Ms[bra1spin]==Ms[ket1spin],SQM[OHead[op],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],0];
result=result/.rhoIndex[particleIndex[bra1_String,particleSpace[bra1space_,bra1spin_],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space_,ket1spin_],ket1rest___]]->If[Ms[bra1spin]==Ms[ket1spin],rhoIndex[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],1];
result=result/.deltaIndex[particleIndex[bra1_String,particleSpace[bra1space_,bra1spin_],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space_,ket1spin_],ket1rest___]]->If[Ms[bra1spin]==Ms[ket1spin],rhoIndex[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],0];
(* 2-body *)
result=result/.SQM[OHead[label_String,indexSymm[0]],particleIndex[bra1_String,particleSpace[bra1space_,bra1spin_],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space_,bra2spin_],bra2rest_],particleIndex[ket1_String,particleSpace[ket1space_,ket1spin_],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space_,ket2spin_],ket2rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin],SQM[OHead[label,indexSymm[0]],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest]],0];
(* 3-body *)
result=result/.SQM[OHead[label_String,indexSymm[0]],particleIndex[bra1_String,particleSpace[bra1space_,bra1spin_],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space_,bra2spin_],bra2rest_],
particleIndex[bra3_String,particleSpace[bra3space_,bra3spin_],bra3rest_],particleIndex[ket1_String,particleSpace[ket1space_,ket1spin_],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space_,ket2spin_],ket2rest_],
particleIndex[ket3_String,particleSpace[ket3space_,ket3spin_],ket3rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin]&&Ms[bra3spin]==Ms[ket3spin],SQM[OHead[label,indexSymm[0]],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],
result=result/.SQM[OHead[op__],particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_]]->If[Ms[bra1spin]==Ms[ket1spin],SQM[OHead[op],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],0];
result=result/.SQS[particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_]]->If[Ms[bra1spin]==Ms[ket1spin],SQS[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],0];
result=result/.rhoIndex[particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest___]]->If[Ms[bra1spin]==Ms[ket1spin],rhoIndex[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],1];
result=result/.deltaIndex[particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest___]]->If[Ms[bra1spin]==Ms[ket1spin],rhoIndex[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest]],0];
(* 2-body ... N.B. SQM and SQS are canonically ordered, hence can just match spins within columns *)
result=result/.SQM[OHead[label_String,indexSymm[0]],particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin],SQM[OHead[label,indexSymm[0]],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest]],0];
result=result/.SQS[particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin],SQS[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest]],0];
(* 3-body ... N.B. see above *)
result=result/.SQM[OHead[label_String,indexSymm[0]],particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],
particleIndex[bra3_String,particleSpace[bra3space__,bra3spin_particleSpin],bra3rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_],
particleIndex[ket3_String,particleSpace[ket3space__,ket3spin_particleSpin],ket3rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin]&&Ms[bra3spin]==Ms[ket3spin],SQM[OHead[label,indexSymm[0]],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],
particleIndex[bra3,particleSpace[bra3space,bra3spin],bra3rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest],
particleIndex[ket3,particleSpace[ket3space,ket3spin],ket3rest]],0];
result=result/.SQS[particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_],particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],
particleIndex[bra3_String,particleSpace[bra3space__,bra3spin_particleSpin],bra3rest_],particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_],
particleIndex[ket3_String,particleSpace[ket3space__,ket3spin_particleSpin],ket3rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin]&&Ms[bra3spin]==Ms[ket3spin],SQS[particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],
particleIndex[bra3,particleSpace[bra3space,bra3spin],bra3rest],particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest],
particleIndex[ket3,particleSpace[ket3space,ket3spin],ket3rest]],0];
result=Expand[result];
Expand All @@ -2667,7 +2683,7 @@ degenerately ordered (e.g. there are duplicates in the lists so that
];

(* optionally replace rhoIndex with deltaIndex ... and reduce the deltas *)
result=If[!rhoToDelta,result,result/.rhoIndex[particleIndex[bra1_String,particleSpace[bra1space_],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space_],ket1rest___]]->(1-deltaIndex[particleIndex[bra1,particleSpace[bra1space],bra1rest],particleIndex[ket1,particleSpace[ket1space],ket1rest]])];
result=If[!rhoToDelta,result,result/.rhoIndex[particleIndex[bra1_String,particleSpace[bra1space__],bra1rest___],particleIndex[ket1_String,particleSpace[ket1space__],ket1rest___]]->(1-deltaIndex[particleIndex[bra1,particleSpace[bra1space],bra1rest],particleIndex[ket1,particleSpace[ket1space],ket1rest]])];
If[ rhoToDelta&&SeQuantDebugLevel>=2,Print["in spintrace: after replacing rhos with deltas = ",TraditionalForm[result]];
];
result=If[!rhoToDelta,result,reduceWick[Expand[result],extInds,intInds]];
Expand Down

0 comments on commit 26b4834

Please sign in to comment.