
---------------------------
(1) MAGNETIC SUSCEPTIBILITY
---------------------------
Input files needed:
- sig.inp1r  (self-energy on real axis)
- sig.inp1i  (self-energy on imaginary axis)
- dmftgkr.in
- dmftgki.in
- suscept.in
------------------
If you want to unfold structure and compute magnetic susceptibility in the extended Brilouin zone,
you must do the following:

1) create case.klist by "klist_generate.py" in the extended BZ. Use all k-points.
   create case.qlist for the q-path where you want your structure factor. You can use
   "qlist_generate.py"

2) correct structure file to have desired atoms in the calculation.
    This is achieved by rearanging grop oprations.

3) edit case.indmfl file to create cluster of atoms, and add "correct" phase 
    factors for unfolding. These phase factors need to be known.


--------------------------------------------------------------------------

4) run "lapw1" on klist generated by "klist_chi.py":
  > x lapw1 -f <case>

5) Print projector:
  > x_dmft.py dmftu -g

---------------------
| Chi on real axis  |
---------------------

6) Create mesh for real axis self-energy, on which chi0 will be computed.
 > LinLogMesh.py -i <input> -o <output>
 will produce "sig.inpr"
 to create sig.inp1r, execute
 > ssplit.py -i sig.inpr -o sig.inp1r

7) Print momentum dependent Green's function on real axis with "dmftgk"
   <prepare dmftgk.in> and execute
   > dmftgk  dmftgr.in
  
  example of dmftgkr.in
"""
g                   # mode [g/e] : we use mode to compute the green's function
BasicArrays.dat     # filename for projector
0                   # matsubara
case.energy         # LDA-energy-file, case.energy(so)(updn)
case.klist          # k-list
case.rotlm          # for reciprocal vectors
Udmft.0             # filename for projector
0.025               # gamma for non-correlated
0.025               # gammac
sig.inp1r           # self-energy name, sig.inp(x)
G_k1r               # name of the output gk
G_local1r           # name of the output glocal
g_local1r           # name of the output glocal-diagonal
"""
8) Unfold Gk on real axis with "unfold.py"
   > unfold.py  dmftgkr.in

  You will obtain G_k1r_ and G_local1r_

9) Compute bubble on the real axis. 

(Optinal:) set number of threads to >= 3 
"export OMP_NUM_THREADS=3"

Prepare input file "suscept.in" and execute:
> dmft_real_bubble suscept.in

example of suscept.in file:
"""
g              # mode [g/e]: we use mode to compute the green's function
case.klist     # filename with k-list
case.qlist     # filename with Qlist
rmesh.dat      # real axis mesh
G_k1r_         # file with real axis k-dependent Gren's function
G_local1r_     # file with real axis local Gren's function
chi0_real.     # name of the Bubble on real axis
G_k1i_         # imaginary axis k-dependent Green's function
G_local1i_     # imaginary axis local Green's function
tvertex.dat    # ctqmc local vertex function
100            # inverse temperature for bose function in Sq(omega)
"""

--------------------------
| Chi on imaginary axis  |
--------------------------

10) Print momentum dependent Green's function on imaginary axis with "dmftgk" :
  sig.inp1 should contain converged imaginary self-energy.
  You can reduce the number of rows (frequencies), however, you should keep twice the number of 
  frequencies kept in vertex calculation. This is because of convolution.
  Copy resulting self-energy into "sig.inp1i"
  Next, prepare <prepare dmftgk.in>  and execute
   > dmftgk dmftgki.in
  
  example of dmftgki.in
"""
g                   # mode [g/e]: we use mode to compute the green's function
BasicArrays.dat     # filename for projector
1                   # matsubara
case.energy         # LDA-energy-file, case.energy(so)(updn)
case.klist          # k-list
case.rotlm          # for reciprocal vectors
Udmft.0             # filename for projector
0.025               # gamma for non-correlated
0.025               # gammac
sig.inp1i           # self-energy name, sig.inp(x)
G_k1i               # name of the output gk
G_local1i           # name of the output glocal
g_local1i           # name of the output glocal-diagonal
"""

11) Unfold Gk on imaginary axis with "unfold.py"
>  unfold.py  dmftgki.in

  You will obtain G_k1i_ and G_local1i_

12) obtain "tvertex.dat" by ctqmc

13) use the above prepared suscept.in to execute 
 > Suscept.py suscept.in

example of suscept.in
"""
case.klist     # filename with k-list
case.qlist     # filename with Qlist
rmesh.dat      # real axis mesh
G_k1r_         # file with real axis k-dependent Gren's function
G_local1r_     # file with real axis local Gren's function
chi0_real.     # name of the Bubble on real axis
G_k1i_         # imaginary axis k-dependent Green's function
G_local1i_     # imaginary axis local Green's function
tvertex.dat    # ctqmc local vertex function
100            # inverse temperature for bose function in Sq(omega)
"""


14) plot Sqw.xx with PlotSq.py

-------------------------------------------------
-------------------------------------------------

---------------------
(2) SUPERCONDUCTIVITY
---------------------
Input files needed:

- sig.inp1i  (self-energy on imaginary axis)
- dmftgki.in
- BCS.in
- case.klist  (with all k-points)
- case.qlist  (with k-points for solving Elliashberg)
- tvertex.dat (ctqmc vertex function)

---------------------------------------------------------------------------
1) create case.klist by "klist_generate.py" in the extended BZ. Use all k-points.
   create case.qlist, which contains a subset of k-points, on which Elliashberg equations
   will be solved. Typically this is Gamma plane or Gamma+Z plane. Use "klist_generate.py" but 
   reduce the number of points iz z-direction.

2) correct structure file to have desired atoms in the calculation.
    This is achieved by rearanging grop oprations.

3) edit case.indmfl file to create cluster of atoms, and add "correct" phase 
    factors for unfolding. These phase factors need to be known.


--------------------------------------------------------------------------

4) run "lapw1" on klist generated by "klist_chi.py":
  > x lapw1 -f <case>

5) Print projector:
  > x_dmft.py dmftu -g


--------------------------
| Chi on imaginary axis  |
--------------------------
6) Print momentum dependent Green's function on imaginary axis with "dmftgk" :
  sig.inp1 should contain converged imaginary self-energy.
  You can reduce the number of rows (frequencies) and keep only frequencies obtain by 
  vertex calculation. Copy resulting self-energy into "sig.inp1i"
  Next, prepare <prepare dmftgk.in>  and execute

   > dmftgk dmftgki.in
  
  example of dmftgki.in
"""
g                   # mode [g/e]: we use mode to compute the green's function
BasicArrays.dat     # filename for projector
1                   # matsubara
case.energy         # LDA-energy-file, case.energy(so)(updn)
case.klist          # k-list
case.rotlm          # for reciprocal vectors
Udmft.0             # filename for projector
0.025               # gamma for non-correlated
0.025               # gammac
sig.inp1i           # self-energy name, sig.inp(x)
G_k1i               # name of the output gk
G_local1i           # name of the output glocal
g_local1i           # name of the output glocal-diagonal
"""

7) Unfold Gk on imaginary axis with "unfold.py"

>  unfold.py  dmftgki.in

  You will obtain G_k1i_ and G_local1i_

8) obtain "tvertex.dat" by ctqmc

9) To get the particle-hole susceptibility-vertex and particle-particle bubble
  you should execute:

 > Elliashberg.py BCS.in

example of BCS.in
"""
case.klist              # filename with k-list
case.qlist              # filename with Qlist
G_k1i_                  # imaginary axis k-dependent Green's function
G_local1i_              # imaginary axis local Green's function
tvertex.dat             # ctqmc local vertex function
0.0                     # Small shift: if magnetic susceptibility is diverging at ordering wave vector Q, we can add a small constant to avoid the existance of the pole
CheckC4Symmetry=True    # Should we check C4 symmetry
EnforceC4Symmetry=True  # Should we enforce C4 symmetry
symm=[2,3]              # which orbitals to symmetrize with C4
rest=[0,1,4]            # which are the rest of the orbitals
imk =[0,1,1]            # we are enforcing G(k)_{symm[0],rest[k]}=G_{R(k)}_{symm[1],rest[k]}. Here imk=0 (imk=1) if R  is clockwise (counter-clockwise) rotation.
"""

You will obtain "Chi0pp.dat", which contains particle-particle bubble, and
"Gpp.dat", "Gpm.dat", "Gmm.dat". The latter contain the 
particle-hole susceptibility-vertex (Gamma*Chi*Gamma) on the case.qlist and at zero energy, 
and at the following brach of the complex plane:
Gpp=(Gamma*Chi*Gamma)_Q(imu=0+,imu'=0+,iOm=0), 
Gpm=(Gamma*Chi*Gamma)_Q(imu=0+,imu'=0-,iOm=0), 
Gmm=(Gamma*Chi*Gamma)_Q(imu=0-,imu'=0-,iOm=0)

------------------
| BCS Equations  |
------------------
10) Use the particle-hole susceptibility vertex Gamma*Chi*Gamma to compute the 
particle-particle vertex, and set-up BCS equations, and solve them with direct diagonalization:

 > BCS.py BCS.in

You will obtain gs_symmetryr.1, gs_symmetryr.2,..., which stand for the pairing field in the leading channel, 
and the subleading channels, respectively.


11) Display the symmetry of pairing field in orbital space with 

> PlotBCS.py gs_symmetryr.1


----------------------
| Fermi surface plot |
----------------------
12) Create a very dense mesh for the Fermi surface plot in the same plane in which BCS equations were solved. 
    It is very imortant to use the same region in k-space for Fermi surface plot, because simple interpolation is used 
    in the code.  You can use "klist2_generate.py".  
   Obtain vector file on this mesh 

  > klist2_generate.py > case.klist
  > x lapw1 -f <case>

13) Print projector:

  > x_dmft.py dmftu -g


14) Obtain eigenvectors and eigenvalues on this k-mesh at the Fermi level. First create a self-energy file with a single row 
  with sigma at zero frequency (sigma(omega=0) here named sig.inp1_zero. The execute:

  > dmftgk dmftgke.in

example of "dmftgke.in":
"""
e                   # mode [g/e]: we use mode to compute eigenvalues and eigenvectors
BasicArrays.dat     # filename for projector
0                   # matsubara
case.energy         # LDA-energy-file, case.energy(so)(updn)
case.klist          # k-list
case.rotlm          # for reciprocal vectors
Udmft.0             # filename for projector
0.025               # gamma for non-correlated
0.025               # gammac
sig.inp1_zero       # self-energy name, sig.inp(x)
eigenvalues.dat     # eigenvalues
UR.dat              # right eigenvector in orbital basis
UL.dat              # left eigenvector in orbital basis
-1.0                # emin for printed eigenvalues
 1.0                # emax for printed eigenvalues
"""

You will obtain eigenvalues in "eigenvalues.dat" and projected left and right eigenvectors "UL.dat", "UR.dat".
The latter are obtained by  UL==Al*U and UR==U^+*Ar , where the local green's function is obtained by
G_local = \sum_k  U^+*Ar * 1/(omega+mu-eigenvalues) * Al*U == \sum_k UR * 1/(omega+mu-eigenvalues) * UL

15) If necessary, unfold the eigenvectors such that they are expressed in the same unit cell as the BCS equations were solved:

  >  unfold.py  dmftgke.in

You will obtain "UR.dat_" and "UL.dat_"

16) Plot fermi surface with orbital character:
 
 >  fsplot.py eigenvalues.dat

17) Plot gap function on the Fermi surface with:
 
  > BCSplot.py gs_symmetryr.1
