If you want to run the code for Si, do the following:
(1) mkdir Si
(2) cd Si
(3) mkdir paratec; mkdir xi; mkdir sig
(4) if you type pwd you should have the subdirectories:
paratec sig xi
paratec -> has the LDA/pseudopotential/planewave code
sig -> has self-energy program
xi -> has screening program
(ignore subdirectory epsilon)
(2) in your paratec directory you should have:
the following input files:
input.cd
input.cwfe
input.cwfeq
input.wfn0
the pseudopotential:
Si_POT.DAT
and the executable:
paratec.mpi
--- to make the executable, paratec.mpi, do the following:
cd ~
mkdir source
cp paratec.tar source/
cd source
tar -xvf paratec.tar (untar the tar file)
cd paratec
make paratec
cd ~/Si/paratec
ln -s ~/source/paratec/bin/paratec.mpi . ! paratec.mpi is the
executable
----To convert the ascii pseudopotential to the binary form in which
it is read, do the following:
cd ~/source/paratec ! same directory as above
make tools
cd ~/Si/paratec
ln -s ~/source/paratec/bin/kbascbin .
./kbascbin Si_POT.ASC.DAT Si_POT.DAT
Now we are ready to run paratec for four times, one time for each
of the four input files:
input.cd ----> compute the ground state density, and screening.
use output_flags gwscreening, and save CD.
input.cwfe -----> compute the wavefunctions on the
4x4x4 grid, shifted 0.5 0.5 0.5. In this
case we compute 5-10 times the number
of valence bands. Notice the option in the input:
output_flags gwr. This means, produce real
wavefunctions to be read by the program in directory
Xi. Also, make sure that you converge the wavefunctions
to a high accuracy.
By the way, for a system with inversion symmetry,
you need to say, output_flags gwc.
input.cwfeq -----> compute the wavefunction on the
4x4x4 grid, shifted 0.5 0.5 0.5 plus some small
q-vector, specified by the option in the input,
e.g. gw_shift 0.00 0.00 0.001. Compute only one more
than the number of valence bands in this case. Use
option output_flags gwr.
input.wfn0 ---> computes wavefunctions on unshifted grid, 5-10 number
of valence wavefunctions, to be read by
In summary, we have the following new features in paratec
which enable the interface between paratec and GW:
gw_shift 0.00 0.00 0.001
output_flags gwr,gwc, or gwscreening
gwr -> produces real wavefunctions
gwc -> produces complex wavefunctions
gwscreening -> produces CD95 and VXC (for the sigma code)
So do the following:
cd ~/Si/paratec ! go to the right directory
cp input.cd input ! compute ground state charge density
mpprun -n 8 ./paratec.mpi ! run paratec
cp CD CD.save ! save ground state charge density
mv CD95 ../sig ! CD95 is the charge density read by self-energy program
mv VXC ../sig ! VXC is the exchange corr
e. read by self-energy program
cp input.cwfe input ! computes wavefunctions on 4x4x4 0.5 0.5 0.5 k grid
! computes about 5-10 times number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi ! run paratec
mv GWR ../xi/CWFE ! move wavefunctions (GWR or GWC) to xi
cp input.cwfeq input ! computes wavefunctions on 4x4x4 0.5 0.5 0.5 +
small q k grid for xi
! only need one more than number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi ! run paratec
mv GWR ../xi/CWFEq
cp input.wfn0 input ! computes wavefunctions for self-energy
! 4x4x4 unshifted
! 5-10 times number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi
mv GWR ../sig/wfn0
(2) calculate static screening matrix:
q is an element of 4x4x4 unshifted grid (8 q-points)
k,k' are elements of 4x4x4 shifted grid (10 k-points)
q is {k-k'| k,k' shifted}
RPA screening (xi polarization)
xi(G,G';q) = Sum(cvk) [ M(cvkqG) M*(cvkqG') /(E_ck - E_vk+q) ]
M(cvkqG) = <ck|exp(i(G+q).r)|vk+q>
eps(G,G';q) = delta_GG' - [4pi/(q+G)^2 ]xi(G,G';q)
if G=0 and q=0
then we let
eps(0,G';q) = delta_0G' - lim(q->0) [4pi/(q)^2] xi(0,G';q)
program takes inverse of eps(G'G;q) for every q:
-> epsmat (for q not equal to gamma)
-> eps0mat (for q equal to gamma)
cd xi ! you should have the following files:
CWFE
CWFEq
xi0.inp input file
xi0 (or cxi0 for complex version) executable
xi0.inp looks like
epsilon_cutoff 6.25 ! compute eps matrix up to G^2,G'^2 < 6.25 Ry)
number_bands 44 ! number of bands in summation (cond.+valence)
wavefunction_cutoff 12.0 ! same as in paratec
number_qpoints 8 ! number of q-points to compute eps(GG';q)
band_occupation 4*1 40*0 ! band occupation
begin qpoints
0.00 0.00 0.001 1.0 1 ! q-point,divisor, 1 or 0
0.00 0.00 0.25 1.0 0 ! copy from paratec
0.00 0.00 0.50 1.0 0
0.00 0.25 0.25 1.0 0
0.00 0.25 0.50 1.0 0
0.00 0.25 0.75 1.0 0
0.00 0.50 0.50 1.0 0
0.25 0.50 0.75 1.0 0
end
#Values for the 1st step
evs 0.0000
evdel 0.0000
ev0 0.0000
ecs 0.0000
ecdel 0.0000
ec0 0.0000
mpprun -n 8 ./xi0 ! run xi
mv eps0mat ../sig
mv epsmat ../sig
always check Xi(0,0;q->0) and eps^(-1)(0,0;q->0)
in xi0.log
1/eps^(-1)(0,0;q->) should be the macroscopic epsilon,
local field effects
1/(1-4*pi*Xi(0,0;q->0) should also be (roughly) the macroscopic epsilon,
non local field effects
(3) Look at Hybertsen/Louie, PRB 34,5390.
formulas:
(34a) we calculate Sigma_(SEX)
has two terms: 1 + Omega/(DeltaE^2+wp^2)
first term -> Gv (small v is bare coulomb interaction)
second term -> is part of G(W-v)
(34b) we calculate Sigma_(COH)
other part of G(W-v)
We compute
<nk|Sigma(E)|nk>, Sigma = i G_0 W
cd sig
in your directory you should have:
eps0mat
epsmat
CD95
VXC
wfn0
sig.inp
sigma (or csigma)
looking at sig.inp:
number_kpoints 2 # of k-points in <nk|Sigma(E)|nk>
gmaxs 2.50000 # G-vector sum for screened interaction, G<gmaxs
gmaxb 5.00000 # G-vector sum for bare interaction G<gmaxb
ncband 44 # number of bands in summation (n1 in notation of paper)
nbnmin 1 # min band to calculate Sigma for (n in notation of paper)
nbnmax 8
dw 1.0 # don't change.
epscut 3.0 # don't change
ieps 0 # don't change
ecut 12.00 # wavefunction cutoff
number_offdiag 3 # number of offidagonal elements
begin kpoints
0.0 0.0 0.0 2.0
0.0 0.0 0.5 1.0
end
begin offdiag
1 3 <1k|Sigma(E)|3k>
1 2
2 1
end
evs 0.0000 ! scissor shift
evdel 0.0000
ev0 0.0000
ecs 0.0000
ecdel 0.0000
ec0 0.0000
mpprun -n 8 ./sigma ! run Sigma
output looks like this (sig.log) :
k= 0.000 0.000 0.000
i elda e0 x sx ch sig v0 efso enew (ev) spin=1
1 3.884 3.884 -17.107 13.287 -7.228 -11.048 -10.429 3.265 3.467
2 15.775 15.775 -12.700 9.382 -8.410 -11.727 -11.229 15.277 15.383
3 15.775 15.775 -12.700 9.382 -8.410 -11.727 -11.229 15.277 15.383
4 15.775 15.775 -12.700 9.382 -8.410 -11.727 -11.229 15.277 15.383
5 18.356 18.356 -5.672 3.732 -7.827 -9.767 -10.051 18.640 18.579
6 18.356 18.356 -5.672 3.732 -7.827 -9.767 -10.051 18.640 18.579
7 18.356 18.356 -5.672 3.732 -7.827 -9.767 -10.051 18.640 18.579
8 19.096 19.096 -5.770 4.011 -8.825 -10.584 -10.804 19.316 19.268
1 1 --real- -17.107 13.287 -7.228 -11.048 -10.429 del= -0.620
1 2 --real- -0.075 0.122 -0.013 0.034 -0.000 del= 0.034
2 1 --real- -0.075 0.060 0.002 -0.012 -0.000 del= -0.012
elda -> E_LDA
e0 -> E_LDA scissor shifted (linear interpolation)
x -> <nk|Gv|nk>, fock operator
sx +ch -> <nk|G(W-v)|nk>
sig = sx + ch
v0 -> <nk|v_xc|nk>
efs0 = e0 + sig - v0
enew = e0 + Z(sig - v0), where Z is close to 1
Sigma depends on E
E^QP = E_0 + <nk|Sigma(E^QP) - vxc|nk>
Assuming Sigma(E^QP) = Sigma(E_LDA) + Sigma'(E_LDA)(E^QP-E^LDA)
then we can show that
E^QP = E_0 + Z <nk|Sigma(E^LDA) - vxc|nk>
where Z = 1/(1-Sigma'(E_LDA))
------------------------------------------------------------------
We calculate
Sigma(E) = G({Psi_LDA},{E_0}) * W({Psi_LDA,{E_0'})
{E_0} is given by E_LDA with scissor shift specified in sig.inp
{E_0'} is given by E_LDA with scissor shift specified in xi0.inp
{E_0} --> some approximation to E^QP.
How scissor shift is defined:
eval=eval+(evs+evdel*(eval-ev0))/
((1.0d0-evdel)*ryd)
econd=econd+(ecs+ecdel*(econd-ec0))/
((1.0d0-ecdel)*ryd)
Updating Xi : epsilon -> 1 + wp^2/E_av^2
doesn't matter. E_av = average gap (which doesn't change much)