amended previous commit

Author: evaleev

SeQuant.m

@@ -2654,23 +2654,26 @@ degenerately ordered (e.g. there are duplicates in the lists so that
 (* zero out integrals that don't conserve spin quantum numbers *)
 (* 1-body *)
 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/.SQS[particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_]]->If[Ms[bra1spin]==Ms[ket1spin],SQS[particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest]],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];
+result=result/.SQS[particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_],particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin],SQS[particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest]],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=result/.SQS[particleIndex[ket1_String,particleSpace[ket1space__,ket1spin_particleSpin],ket1rest_],
+particleIndex[ket2_String,particleSpace[ket2space__,ket2spin_particleSpin],ket2rest_],
+particleIndex[ket3_String,particleSpace[ket3space__,ket3spin_particleSpin],ket3rest_],
+particleIndex[bra3_String,particleSpace[bra3space__,bra3spin_particleSpin],bra3rest_],
+particleIndex[bra2_String,particleSpace[bra2space__,bra2spin_particleSpin],bra2rest_],
+particleIndex[bra1_String,particleSpace[bra1space__,bra1spin_particleSpin],bra1rest_]]->If[Ms[bra1spin]==Ms[ket1spin]&&Ms[bra2spin]==Ms[ket2spin]&&Ms[bra3spin]==Ms[ket3spin],SQS[particleIndex[ket1,particleSpace[ket1space,ket1spin],ket1rest],particleIndex[ket2,particleSpace[ket2space,ket2spin],ket2rest],
+particleIndex[ket3,particleSpace[ket3space,ket3spin],ket3rest],
+particleIndex[bra3,particleSpace[bra3space,bra3spin],bra3rest],particleIndex[bra2,particleSpace[bra2space,bra2spin],bra2rest],particleIndex[bra1,particleSpace[bra1space,bra1spin],bra1rest]],0];
 result=Expand[result];
 
 If[ SeQuantDebugLevel>=2,Print["in spintrace: after zeroing out spin-nonconserving integrals = ",TraditionalForm[result]];