====== MeanFieldGroundstate ======

###
//MeanFieldGroundState($O$, $\rho$, $T$)// calculates the ground-state of operator $O$ within mean-field theory, starting from density matrix $\rho$ at temperature $T$.
$rho$ stores the expectation values of $a^{\dagger}_{\tau}a^{\phantom{\dagger}}_{\tau'}$, a table of dimensions $NFermion$ by $NFermion$. The operator is expected to consist of:
$O^{(0)}$ a constant,
$O^{(1)}$ a one particle operator of the form $\sum_{\tau,\tau'} \alpha_{\tau,\tau'} a^{\dagger}_{\tau}a^{\phantom{\dagger}}_{\tau'}$,
$O^{(2)}$ a two particle operator of the form $\sum_{\tau,\tau',\tau'',\tau'''} U_{\tau,\tau',\tau'',\tau'''} a^{\dagger}_{\tau}a^{\dagger}_{\tau'}a^{\phantom{\dagger}}_{\tau''}a^{\phantom{\dagger}}_{\tau'''}$.

The two particle operator will be replaced in mean-field, using the Hartree-Fock approximation by:
\begin{eqnarray}
a^{\dagger}_{i}a^{\dagger}_{j}a^{\phantom{\dagger}}_{k}a^{\phantom{\dagger}}_{l} &\to&\\
\nonumber &-& a^{\dagger}_{i}a^{\phantom{\dagger}}_{k} \langle a^{\dagger}_{j}a^{\phantom{\dagger}}_{l} \rangle \\
\nonumber &+& a^{\dagger}_{i}a^{\phantom{\dagger}}_{l} \langle a^{\dagger}_{j}a^{\phantom{\dagger}}_{k} \rangle \\
\nonumber &+& a^{\dagger}_{j}a^{\phantom{\dagger}}_{k} \langle a^{\dagger}_{i}a^{\phantom{\dagger}}_{l} \rangle \\
\nonumber &-& a^{\dagger}_{j}a^{\phantom{\dagger}}_{l} \langle a^{\dagger}_{i}a^{\phantom{\dagger}}_{k} \rangle \\
\nonumber &-& \langle a^{\dagger}_{i}a^{\phantom{\dagger}}_{l} \rangle \langle a^{\dagger}_{j}a^{\phantom{\dagger}}_{k} \rangle \\
\nonumber &+& \langle a^{\dagger}_{i}a^{\phantom{\dagger}}_{k} \rangle \langle a^{\dagger}_{j}a^{\phantom{\dagger}}_{l} \rangle 
\end{eqnarray}

the function MeanFieldGroundState searches for a self consistent solution such that the lowest eigen-state of the mean-field approximated operator has the same density matrix as is used to calculate the operator.
###

===== Input =====

  * $O$ : Operator
  * $rho$ : Matrix (Table of Table of length $O.NF$) of doubles
  * $T$ : Real

===== Output =====

  * $rho_{MF}$ : Matrix (Table of Table of length $O.NF$ The self consistent density matrix
  * $E_{MF}$ : The mean-field ground state energy
  * $O_{MF}$ The mean-field approximated operator

===== Example =====

###
description text
###

==== Input ====
<code Quanty MeanFieldGroundstate.Quanty>
dofile("../definitions.Quanty")

-- we define an arbitrary operator
Opp = -(2*OppSy + OppLy)*(2*OppSy + OppLy) + (2*OppSy + OppLy) + 0.0000001 * OppLy
-- and a starting density
rho = DensityMatrix(psi1)
-- as well as a temperature needed to average over degenerate states
T = 0.0001
-- calculate the ground-state in mean-field
rhogrd, E, OppMF = MeanFieldGroundState(Opp, rho, T)

-- print the resulting density
print(Chop(rhogrd))
-- print the ground-state energy
print(E)
-- print the Hamiltonian in mean-field, i.e. a potential optimized for the ground-state density
print(Chop(OppMF))

-- lets compare the eigenstates of Opp and OppMF

Npsi=15

-- In order to make sure we have a filling of 2
-- electrons we need to define some restrictions
StartRestrictions = {Nf, Nb,  {"111111",2,2}}

-- We now can create the lowest Npsi eigenstates:
psiList   = Eigensystem(Opp,   StartRestrictions, Npsi)
psiListMF = Eigensystem(OppMF, StartRestrictions, Npsi)


-- We define a list of some operators to look at their expectation values
oppList={Opp, OppMF, OppSy, OppLy}

-- after we've created the eigen states we can look 
-- at a list of their expectation values
-- on the left we show the full eigen-states, on the right the eigen-states of the mean-field approximated operator
print("  <E>     <E>     <S_y>   <L_y>  MF  <E>     <E>     <S_y>   <L_y>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.3f ",Chop(psiList[i]*oppList[j]*psiList[i])))
  end
  io.write(" | ")
  for j=1,#oppList do
    io.write(string.format("%7.3f ",Chop(psiListMF[i]*oppList[j]*psiListMF[i])))
  end 
  io.write("\n")
end
</code>

==== Result ====
<file Quanty_Output MeanFieldGroundstate.out>
{ { 0.37500001606196 , (-5.3632883498142e-09 + 0.375 I) , (-7.5782580190366e-09 + 0.17677669529663 I) , (-0.17677671801051 + 8.4602858479628e-12 I) , (0.12499998393969 + 1.1964467459925e-11 I) , (5.3539869430386e-09 + 0.125 I) } , 
  { (-5.3632883498142e-09 - 0.375 I) , 0.37499998393804 , (0.17677667258276 + 8.460206504589e-12 I) , (7.5782580958409e-09 + 0.17677669529663 I) , (5.3539868459665e-09 - 0.125 I) , (0.12500001606031 - 1.1964561947813e-11 I) } , 
  { (-7.5782580190366e-09 - 0.17677669529663 I) , (0.17677667258276 - 8.460206504589e-12 I) , 0.24999996787773 , (1.0707423359528e-08 + 0.25 I) , (7.564325232475e-09 + 0.17677669529663 I) , (-0.17677667258276 - 8.4601954064865e-12 I) } , 
  { (-0.17677671801051 - 8.4602858479628e-12 I) , (7.5782580958409e-09 - 0.17677669529663 I) , (1.0707423359528e-08 - 0.25 I) , 0.25000003212227 , (0.17677671801051 - 8.4602042449046e-12 I) , (-7.5643252261757e-09 + 0.17677669529663 I) } , 
  { (0.12499998393969 - 1.1964467459925e-11 I) , (5.3539868459665e-09 + 0.125 I) , (7.564325232475e-09 - 0.17677669529663 I) , (0.17677671801051 + 8.4602042449046e-12 I) , 0.37500001606196 , (-5.3435844501437e-09 + 0.375 I) } , 
  { (5.3539869430386e-09 - 0.125 I) , (0.12500001606031 + 1.1964561947813e-11 I) , (-0.17677667258276 + 8.4601954064865e-12 I) , (-7.5643252261757e-09 - 0.17677669529663 I) , (-5.3435844501437e-09 - 0.375 I) , 0.37499998393804 } }
-12.0000001

Operator: Operator
QComplex         =          1 (Real==0 or Complex==1 or Mixed==2)
MaxLength        =          2 (largest number of product of lader operators)
NFermionic modes =          6 (Number of fermionic modes (site, spin, orbital, ...) in the one particle basis)
NBosonic modes   =          0 (Number of bosonic modes (phonon modes, ...) in the one particle basis)

Operator of Length   2
QComplex      =          1 (Real==0 or Complex==1)
N             =         36 (number of operators of length   2)
C  0 A  0 |  1.999999871509219E+00 -5.054768245166573E-25
C  0 A  1 |  4.286855070714231E-08  5.499999999999948E+00
C  0 A  2 |  6.059819870986585E-08  3.535533976643392E+00
C  0 A  3 |  1.817109998858424E-07  3.384084768552549E-11
C  0 A  4 |  1.284874287521554E-07  4.785818575907231E-11
C  0 A  5 | -4.283024376372977E-08  4.999999999999833E-01
C  1 A  0 |  4.286855070714231E-08 -5.499999999999948E+00
C  1 A  1 |  2.000000128490683E+00  3.217037880363704E-25
C  1 A  2 |  1.817109995805311E-07  3.384080224126484E-11
C  1 A  3 | -6.059819876404298E-08  3.535533976643392E+00
C  1 A  4 | -4.283024384034779E-08 -4.999999999999832E-01
C  1 A  5 | -1.284874284745996E-07 -4.785812149118257E-11
C  2 A  0 |  6.059819870986585E-08 -3.535533976643392E+00
C  2 A  1 |  1.817109995805311E-07 -3.384080224126484E-11
C  2 A  2 |  2.000000256978111E+00  2.732361740195914E-25
C  2 A  3 | -8.565938661503036E-08  4.999999999999965E+00
C  2 A  4 | -6.054246758567044E-08  3.535533976643392E+00
C  2 A  5 | -1.817109991364418E-07 -3.384080224126484E-11
C  3 A  0 |  1.817109998858424E-07 -3.384084768552549E-11
C  3 A  1 | -6.059819876404298E-08 -3.535533976643392E+00
C  3 A  2 | -8.565938661503036E-08 -4.999999999999965E+00
C  3 A  3 |  1.999999743021791E+00 -5.054768245166573E-25
C  3 A  4 | -1.817109997748201E-07 -3.384084768552549E-11
C  3 A  5 |  6.054246753149332E-08  3.535533976643392E+00
C  4 A  0 |  1.284874287521554E-07 -4.785818575907231E-11
C  4 A  1 | -4.283024384034779E-08  4.999999999999832E-01
C  4 A  2 | -6.054246758567044E-08 -3.535533976643392E+00
C  4 A  3 | -1.817109997748201E-07  3.384084768552549E-11
C  4 A  4 |  1.999999871509220E+00 -7.122719776549142E-25
C  4 A  5 |  4.278973491884084E-08  5.499999999999948E+00
C  5 A  0 | -4.283024376372977E-08 -4.999999999999833E-01
C  5 A  1 | -1.284874284745996E-07  4.785812149118257E-11
C  5 A  2 | -1.817109991364418E-07  3.384080224126484E-11
C  5 A  3 |  6.054246753149332E-08 -3.535533976643392E+00
C  5 A  4 |  4.278973491884084E-08 -5.499999999999948E+00
C  5 A  5 |  2.000000128490682E+00 -9.188651824014346E-26


  <E>     <E>     <S_y>   <L_y>  MF  <E>     <E>     <S_y>   <L_y>
-12.000 -12.000  -1.000  -1.000  | -12.000 -12.000  -1.000  -1.000 
 -6.000  -6.000   0.000  -2.000  |  -6.000  -6.000   0.000  -2.000 
 -6.000  -6.000  -1.000   0.000  |  -6.000  -6.000  -1.000   0.000 
 -6.000  20.000   1.000   1.000  |  -2.000  -2.000   0.000  -1.000 
 -2.000  -1.737   0.000  -1.000  |  -2.000  -2.000  -1.000   1.000 
 -2.000  -0.263   0.000  -1.000  |  -2.000  -0.000   0.000  -1.000 
 -2.000  14.000   1.000   0.000  |   0.000   4.000   0.000   0.000 
 -2.000  -2.000  -1.000   1.000  |   0.000   4.000   0.000   0.000 
 -2.000  14.000   0.000   2.000  |   0.000   4.000   0.000   0.000 
 -0.000  10.000   1.000  -1.000  |   0.000   8.000   0.000   1.000 
  0.000   4.000   0.000   0.000  |  -0.000  10.000   1.000  -1.000 
  0.000   4.000   0.000   0.000  |   0.000  10.000   0.000   1.000 
  0.000   4.000   0.000   0.000  |  -2.000  14.000   1.000   0.000 
  0.000   8.831   0.000   1.000  |  -2.000  14.000   0.000   2.000 
  0.000   9.169   0.000   1.000  |  -6.000  20.000   1.000   1.000 
</file>

===== Table of contents =====
{{indexmenu>.#1}}