nmr_q 0.001
sets the 
(in cartesian coordinates), which determines the
periodicity of the applied external field.
Default: nmr_q 0.01
Set the output flag nmrshift to get a plot of the susceptibility (section 6.9).
If a line
checkpoint nmr
is entered, the code will checkpoint the NMR shift calculation after each kpoint. You can recover from a crash by giving an input_flags of chkptchi.
IMPORTANT: Right now, the NMR shift can only be computed for INSULATORS. You have to set number_bands to contain only the VALENCE BANDS, or else the result will be wrong.
An entry pw_job band_structure computes the band structure. You should set screening_type to previous, and provide a converged CD file for the starting density. THERE WILL BE NO SELFCONSISTENT ITERATION, i.e. the bandstructure is plotted corresponding to the CD file, or to the atomic valence charge if no intact CD file is found. Furthermore, YOU HAVE TO USE THE FERMILEVEL OF A PREVIOUS SCF RUN. You do so by giving a line like:
fermi_level fermi_energy_in_Ryd
Default: fermi_level 0
If you want to use a different number of bands from that used in an earlier SCF loop, you may specify this number using:
number_alt_bands <number of bands for the bandstructure>
You also have to determine the lines in the BZ along which you wish to plot the bandstructure. This is done by specifying a sequence of lines in the following way:
begin bandstructure label left_middle_right kpoint startx starty startz endx endy endz number_of_bins label left_middle_right kpoint startx starty startz endx endy endz number_of_bins .. label left_middle_right kpoint startx starty startz endx endy endz number_of_bins end bandstructure
NOTE: The syntax of this structure has changed in version 5.0. Previously this ended with just end, and this has changed, so that now you need to end with end bandstructure.
All coordinates are given with respect to the reciprocal lattice vectors. The labels are mandatory; they have to be triples, separated by an underscore. Notice that xmgr allows greek characters using escape sequences (see example below).
The kpoints will be generated into the file KPOINTS. The output of the calculation can be found in the file BANDSTRUC. Energies are in eV, and are plotted against an x-coordinate that corresponds to the physical distances in cartesian k-space, i.e. 1/a.u. The file can be displayed with the popular graphics program ``xmgr''.
Example:
begin bandstructure label \8G\4_\8\D\4_X kpoint 0.0 0.0 0.0 0.0 0.0 0.5 10 label _V_W kpoint 0.0 0.0 0.5 0.0 0.5 0.5 10 end bandstructure
The density of states is output in units of states/eV/unit cell. It does not count both spin states. 2.1).
begin pw_jobs pw_job scf pw_job band_structure end pw_jobs output_flags momdens
You also need to specify the Fermi level according to a previous scf run: fermi_level fermi_energy_in_Ryd
If you don't specify this correctly, you will get garbage!
GW-interfaced Version of Paratec with f-states was implemented by Eric Chang. Further development was done by Xavier Blase, Gian-Marco Riganese and others.
To calculate dielectric function, we need 2 wavefunctions from paratec. They are called CWFE and CWFEq. CWFE requires lots of conductance bands, but CWFEq only requires wavefunction of bands (valence bands + 2).
Procedure:
(1) Generate a converged charge density CD, save this CD to a file CD.conv.
You must use the same charge density each time when you perform the following jobs.
(2) Compute CWFE. In the input file of paratec, use the options
number_bands 90 # this number means the number of conductance bands used in output_flags gwr # the summation of dielectric function calculationIf the system does not have inversion symmetry, you need to run the complex version:
number_bands 90 # this is a large number compared to the number of valence bands output_flags gwcAfter run paratec, rename the file GWR (or GWC for complex version) to CWFE. You may need to use fewer or more bands in the option number_bands.
If one wants to choose only selected wavefunctions for the GW calculation then the following procedure may be implemented. For selected bands around the Fermi energy one uses
num_gwout 4in order to choose 4 bands on either side of the Fermi energy. If one wants a different energy than the Fermi energy to obtain wavefunctions near then one uses
gwout_mid_energy -0.45If instead of defining a number of states in the above manner, one wants to obtain all wavefunctions within an energy window then
gwout_low_energy -1.0 gwout_high_energy 2.0is employed.
(3) Compute CWFEq. In the input file of paratec, use the options
number_bands 5 # this number is (number of valence bands + 2) output_flags gwr # use gwc if you want to run complex version of the GW code gwshift 0.005 0.005 0.005You should use number_bands slightly more than the number of valence bands in the system. gwshift gives a wavefuction at slightly shifted grid. Rename the file GWR to CWFEq after paratec computation.
GWC or GWR 
CWFE for the shifted grid (need many bands)
GWC or GWR CWFEq for slightly shifted+q grid
(need only valence bands,need to read k-points from a file)
(4) Compute VXC and CD95. VXC is the information of exchange-correlation functional;
CD95 is the charge density. In the input file of paratec, use the option
output_flags gwscreening(5) Compute wfn0. In the input file of paratec, use the option
k_shift 0 0 0 number_bands 90 # test the convergence of this number output_flags gwr # use gwc if you want to run complex version of the GW codeRename GWR or GWC to wfn0.
(6) Now you are ready to use CWFE and CWFEq to compute eps0mat and epsmat. eps0mat, epsmat, wfn0, VXC, CD95 are required to do GW calculation.
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' - [8pi/(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) [8pi/(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 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->0) should be the macroscopic epsilon,
local field effects
1/(1-8*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: delta_GG' + Omega^2/(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
looking at sig.inp:
screened_coulomb_cutoff 6.25 ! same as in xi0.inp
bare_coulomb_cutoff 12.0 ! never less than wavefunction_cutoff
wavefunction_cutoff 12.0 ! same as in paratec
number_bands 44 ! These two lines have the same meaning as
band_occupation 4*1 40*0 ! in xi0.inp
emergency_cutoff 6
#Define the bands for which we want to compute sigma
band_index_min 1 ! band_index_min and band_index_max fix the range
band_index_max 10 ! of calculation in terms of bands
#finite_difference_spacing 0.1
number_kpoints 1
begin kpoints
0.0000 0.0000 0.0000 1.0 ! Gamma
end
number_offdiag 3 ! choose off-diagonal entries of Epsilon to look at
begin offdiag
1 3 <1k|Sigma(E)|3k>
1 2
2 1
end
evs 0.0000 ! scissor shift (energy unit is eV)
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 = x + 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))
econd --> econd+(ecs+ecdel*(econd-ec0))
Updating Xi : epsilon -> 1 + wp^2/E_av^2
doesn't matter. E_av = average gap (which doesn't change much)
tensor_bra 0 1 0 tensor_ket 0 0 1to select the yz component (in cartesian coordinates) of a tensor. The ``tensor'' input is also used for specifying the alignment of the external B-field for plotting purposes.
tensor_bra 1 0 0 tensor_ket 1 0 0
Example:
begin line_plot line 0.0 0.0 0.0 1.0 1.0 1.0 100 plot data along (111) on 100 mesh points end line_plot
The start and end coordinates of line_plot are relative to the realspace lattice vectors.
Paratec can produce Data Explorer and Khoros files that visualize wavefunctions at particular kpoints, bands. The square of the wavefunctions is plotted in realspace. Wave function plotting is switched on by setting pw_job scf and output_flags waveplot.
The wave functions to be plotted are specified as pairs of the irreducible kpoint number and band number.
Example:
begin plot_wave_function wavefunction 1 12 irreducible kpoint 1, band 12 wavefunction 2 10 irreducible kpoint 2, band 10 end plot_wave_function
To produce a simple three-dimensional charge density plot, one simply adds the flag cdplot to the list of output flags, and runs an scf calculation. If an energy window is specified with the energy window structure:
begin energy_window window (start_energy in eV) (end_energy in eV) end energy_windowthen the charge density for the plot is generated only by using wave functions with eigenvalues within the energy window (specified in eV). To find suitable values for start_energy and end_energy, look at the eigenvalues printed out in previous runs. A gaussian smearing is applied, and the occupation numbers are ignored. This feature only works if the wave functions have been computed, i.e. it does not work with the pot_plot plane wave job.
Warning: the energy window specification does not include an adjustment for the fermi level. Also, unoccupied bands will be included in the charge density if your energy window includes energies above the fermi level.
Both line plotting and three-dimensional plotting of the effective Kohn-Sham potential (without the nonlocal pseudopotential part) are available via the pot_plot plane-wave job. For more see section 5.1.1.
A ball-and-stick file BALLS.dx or BALLS.kh can be generated and viewed with the IBM Data Explorer or Khoros. This is done by adding the output flag ballnstick to the output_flags field (section 6.9). To use the data explorer, please refer to the README file in /usr/local/codes/paratec/dx. To use Khoros instead, set output_flags khoros.
number_of_spins 1
(for non-spin polarized LDA)
or
number_of_spins 2
(for spin-polarized LSDA)
Default: number_of_spins 1
To get a spin-polarized starting configuration for the self-consistent charge, use the spin polarization factors where you input the crystal structure.
A magnetic field
 |
(6.1) |
can be applied to stabilize for example an antiferromagnetic phase. The magnetic field acts only on the spin, AND ALTERS THE ENERGY AND WAVEFUNCTIONS!! This is done with an entry like:
magnetic_field_kvec 1 1 0 1
The last number switches the field on (=1) or off (=0), the first 3 INTEGER numbers specify the wave number in reciprocal lattice coordinates.
The amplitude (normalized such that 
Ryd) is entered as:
magnetic_field_strength 0.001
The option input_flags enables you to read additional input from various files. The flags for the various inputs you want to read are concatenated by underscore.
input_flags flag1_flag2_flag3 (sets 3 input flags)
Implemented flags:
For the direct energy minimization (optimize insulator), it is often useful to checkpoint the wave functions and charge density during the iteration. This can be done with
checkpoint_emin frequencyExample:
checkpoint_emin 50 writes the wave functions to disk every 50 iterations. By default, no checkpointing is done.
Checkpointing can also be controlled for the SCF method. By default the charge density and potential at each SCF cycle are not written to disk as this is too costly. They are written at convergence or when the maximum number of cycles is reached. If one want them checkpointed then use
io_scf true
The option output_flags enables you to get additional numbers out of the code. The flags for the various outputs you want to obtain are concatenated by underscore.
output_flags flag1_flag2_flag3 (sets 3 output flags)
Implemented flags:
in
realspace. If requested, line plots are done also.
WARNING: paratec names the wavefn files using only one digit for the k points, so if you plot the same band at multiple k points (which have an index greater than 9), the wavefunctions will overwrite all but the first. This is a bug, but you can work around it, and I (David Roundy) don't feel like fixing it right now. Sorry.