Test simplification routines
Part of SNEG package, http://nrgljubljana.ijs.si/sneg
Rok Zitko, rok.zitko@ijs.si, March 2007

<< sneg`sneg`;

snegfermionoperators[c, d, e, o]
snegrealconstants[δ, U, t, U12, J12]

SnegSimplify and related functions try to transform an expression written using the second quantization operators into a simpler expression in terms of  "higher level functions", such as charge, spin, hopping, etc. operators. The operators are wrapped in HoldForm[] forms which can be removed using HoldForm[i_]->i transformation rule.

number[c[]]
x = SnegSimplify[%]
x // FullForm
x /. HoldForm[i_] -> i

Single-impurity Anderson model

expr = δ (number[c[]] - 1) + (U/2) 
            pow[number[c[]] - 1, 2] // Expand
SnegSimplifyNumber[%]
% // FullSimplify

FullSimplify[expr, 
      TransformationFunctions -> {Automatic, SnegSimplifyNumber, 
          SnegSimplifyHubbard}]

      
FullSimplify[
      expr == (% //. HoldForm[i_] :> i)]

number[c[], UP] + number[c[], DO]
SnegSimplifyNumber[%]
2 (number[c[], UP] + number[c[], DO])
SnegSimplifyNumber[%]
Expand[2 (number[c[], UP] + number[c[], DO])]
SnegSimplifyNumber[%]

spinx[c[]]
SnegSimplifySpin[%]
spiny[c[]]
SnegSimplifySpin[%]
spinz[c[]]
SnegSimplifySpin[%]

spinx[c[]] // Expand
SnegSimplifySpin[%]
spiny[c[]] // Expand
SnegSimplifySpin[%]
spinz[c[]] // Expand
SnegSimplifySpin[%]

2  spinx[c[]]
SnegSimplifySpin[%]
2  spiny[c[]]
SnegSimplifySpin[%]
2  spinz[c[]]
SnegSimplifySpin[%]

I spinx[c[]] // Expand
SnegSimplifySpin[%]
I spiny[c[]] // Expand
SnegSimplifySpin[%]
I spinz[c[]] // Expand
SnegSimplifySpin[%]

U spinx[c[]]
SnegSimplifySpin[%]
U spiny[c[]]
SnegSimplifySpin[%]
U spinz[c[]]
SnegSimplifySpin[%]

Expand[U spinx[c[]]]
SnegSimplifySpin[%]
Expand[U spiny[c[]]]
SnegSimplifySpin[%]
Expand[U spinz[c[]]]
SnegSimplifySpin[%]

s1 = spinx[c[]] // Expand
SnegSimplifySpin[s1]
s1 = spinx[c[k]]
SnegSimplifySpin[s1]
s1 = spinx[c[k, l]]
SnegSimplifySpin[s1]

For spin-z operator, there are two possible simplifications.

s3 = spinz[c[]]
SnegSimplifyNumber[s3]
SnegSimplifySpin[s3]

s3 = spinz[c[]] // Expand
SnegSimplifyNumber[s3]
SnegSimplifySpin[s3]

nc[spinx[c[]], spinx[d[]]]
SnegFullSimplify[%]
nc[spiny[c[]], spiny[d[]]]
SnegFullSimplify[%]
nc[spinz[c[]], spinz[d[]]]
SnegSimplifySpin[%]

nc[spinx[c[]], spinx[d[]]] // Expand
SnegFullSimplify[%]
nc[spiny[c[]], spiny[d[]]] // Expand
SnegFullSimplify[%]
nc[spinz[c[]], spinz[d[]]] // Expand
SnegSimplifySpin[%]

On-site spin-spin correlation can be expressed in terms of charge operators

spinspin[c[], c[]]
Simplify[%, 
      TransformationFunctions -> {Automatic, SnegSimplifyNumber, 
          SnegSimplifyHubbard}]

sxz = nc[spinx[c[]], spinz[d[]]]
SnegSimplifySpin[sxz]
FullForm[sxz == (% /. HoldForm[i_] -> i)]

sxx = nc[spinx[c[]], spinx[d[]]]
SnegSimplifySpin[sxx]

s2a = nc[spinx[c[]], spiny[d[]]]
s2b = nc[spiny[c[]], spinx[d[]]]
s2 = Expand[s2a + s2b]
SnegSimplifySpin[s2]

s = spinx[c[]]~nc~spinx[d[]]
SnegSimplifySpin[s]

s = 
      spinx[c[]]~nc~spinx[d[]] + spiny[c[]]~nc~spiny[d[]] // 
        Expand
SnegSimplifySpin[%]

s12 = spinxyz[c[]]~inner~spinxyz[d[]]
SnegSimplifySpin[s12]

s12 = 
      spinxyz[c[]]~inner~spinxyz[d[]] // Expand
SnegSimplifySpin[s12]

hop[c[], d[]]
SnegSimplifyHop[%]
hop[c[3], d[3]]
SnegSimplifyHop[%]
hop[c[k], d[l, m]]
  
SnegSimplifyHop[%]

twohop[c[], d[]]
SnegSimplifyHop[%]

Expand [twohop[e[], o[]]]
SnegSimplifyHop[%]

Expand [2  twohop[e[], o[]]]
SnegSimplifyHop[%]

Expand [hop[e[], o[]]]
SnegSimplifyHop[%]

Expand [2  hop[e[], o[]]]
SnegSimplifyHop[%]

hubbard[c[]]
SnegSimplifyHubbard[%]

expr = -hubbard[c[]]
SnegSimplifyHubbard[expr]

hubbard[c[1]]
SnegSimplifyHubbard[%]

hubbard[c[k, l]]
SnegSimplifyHubbard[%]

Special rule (???)

xx = +2 spinz[e[]]~nc~spinz[o[]] - 1/2  chargecharge[e[], o[]] // 
        Expand
SnegSimplifyMisc[%]

chargecharge[e[], o[]]
SnegSimplifyMisc[%]

-1/2  chargecharge[e[], o[]]
SnegSimplifyMisc[%]

xx1 = -2 spinz[e[]]~nc~spinz[o[]]
xx2 = -1/2  chargecharge[e[], o[]]
Expand[xx1 + xx2]
SnegSimplifyMisc[%]
Expand[2 (xx1 + xx2)]
SnegSimplifyMisc[%]

Expand[2 (U + t) (xx1 + xx2)]
SnegSimplifyMisc[%]
Simplify[%]

Two site Hubbard model

H = δ (number[c[]] + number[d[]]) + 
        U/2 pow[number[c[]] - 1, 2] + U/2 pow[number[d[]] - 1, 2] + 
        t hop[c[], d[]]

Timing @ 
      Simplify[H, 
        TransformationFunctions -> 
          allsnegsimplifyrules]
FullSimplify[
      H == (%[[2]] /. HoldForm[i_] -> i)]

Timing @ SnegSimplify[H]

Expand[H]
temp = Simplify[H]
Timing @
      Simplify[temp, 
        TransformationFunctions -> allsnegsimplifyrules]

We switch to even and odd combination of atomic orbitals

Hprime = 
      H /. {c[i__] -> (1/(Sqrt[2])) (e[i] + o[i]), 
            d[i__] -> (1/(Sqrt[2])) (e[i] - o[i])} // Expand

HJ = J12 spinspin[c[], d[]]
HJ = 
        HJ /. {c[i__] -> (1/(Sqrt[2])) (e[i] + o[i]), 
              d[i__] -> (1/(Sqrt[2])) (e[i] - o[i])} // 
          Expand;
x12 = 
        FullSimplify[%, 
          TransformationFunctions -> 
            allsnegsimplifyrules];
FullSimplify[%]
Expand[%]

Two site extended Hubbard model

H = δ (number[c[]] + number[d[]]) + 
        U/2 pow[number[c[]] - 1, 2] + U/2 pow[number[d[]] - 1, 2] + 
        t hop[c[], d[]] + U12 chargecharge[c[], d[]] + 
        J12 spinspin[c[], d[]]

Hprime = 
      H /. {c[i__] -> (1/(Sqrt[2])) (e[i] + o[i]), 
            d[i__] -> (1/(Sqrt[2])) (e[i] - o[i])} // Expand

\(\(\( (*\ exprin : \ what\ needs\ to\ be\ simplified; \ 
      exprout : \ what\ was\ already\ simplified; \ 
      simplifyfnc : \ 
        what\ function\ simplification\ function\ will\ be\ \(\(applied\)\(.\)\)\ *) \)\(
\)\(\(SnegSimplifyOne::error = "Assertion check failed. Bug detected.";\)

    \(SnegSimplifyOne[{exprin_, exprout_}, simplifyfnc_] := 
        Module[{col, res, ok, remainder}, 
col = 
            Collect[exprin, 
              Union[Cases[%, _nc, {1, ∞}]]]; 
Print["COL=", col]; 
res = 
            simplifyfnc[col]; 

 (*\ 
            Bug\ trap\ *) 
cmp1 = 
            Expand[exprin]; 
cmp2 = 
            Expand[res /. 
                HoldForm[x_] :> x]; 
If[
            Simplify[cmp1 - cmp2] ≠ 
              0, 
Print["cmp1=", 
              cmp1]; 
Print["cmp2=", 
              cmp2]; 
Message[
              SnegSimplifyOne::error]; 
Return[
              Null];
]; 

res = 
            Expand[res]; 
If[\(! AtomQ[
                res]\), 
ok = 
              Select[res, 
                MemberQ[#, 
                    HoldForm[_], {0, ∞}] &]; 
remainder = Expand[res - ok], 
 (*\ 
              else\ *) 
ok = 
              res; 
remainder = 
              0;
]; 
{remainder, 
            exprout + ok}
];\)\)\)\)

SnegSimplifyOne[{Hprime, 0}, 
      SnegSimplifyNumber]
SnegSimplifyOne[%, SnegSimplifyHubbard]
SnegSimplifyOne[%, SnegSimplifyHop]
SnegSimplifyOne[%, SnegSimplifyMisc]
SnegSimplifyOne[%, 
      FullSimplify[#, 
          TransformationFunctions -> {Automatic, 
              SnegSimplifySpin}] &]
SnegSimplifyOne[%, SnegSimplifyMisc]
Plus  @@  %
res = % //. ruleSnegSimplifySpin[[1]]

Collect[res, listrealconstants]

Collect[res, 
      Union[Cases[res, HoldForm[_], {1, ∞}]]]