CRPA of SrVO3: Difference between revisions

From VASP Wiki
VisualWikitext
 
(70 intermediate revisions by 4 users not shown)
Line 1: Line 1:
{{Template:gw}}
{{Template:GW - Tutorial}}
The following tutorial describes how to perform [[Constrained–random-phase–approximation_formalism|cRPA calculations]], which is available as of VASP 6.
== Task ==  
== Task ==  
 
Calculation of the Coulomb matrix elements <math>U_{ijkl}(\omega=0)</math> in the constrained Random Phase Approximation ([[Constrained-random-phase approximation|cRPA]]) of SrVO<sub>3</sub> between the Vanadium t<sub>2g</sub> states.
Calculation of the Coulomb matrix elements U<sub>ijkl</sub> in the constrained Random Phase Approximation ([[Constrained Random Phase Approximation|CRPA]]) of SrVO<sub>3</sub>.
----
----


Performing a CRPA calculation with VASP is a 3-step procedure: a DFT groundstate calculation, a calculation to obtain a number of virtual orbitals, and the actual CRPA calculation itself.
Performing a cRPA calculation with VASP is a 3-step procedure: a DFT groundstate calculation, a calculation to obtain a number of virtual orbitals, and the actual cRPA calculation itself.


'''N.B.:''' This example involves quite a number of individual calculations.
'''N.B.:''' This example involves quite a number of individual calculations.
The easiest way to run this example is to execute:
The easiest way to run this example is to execute:
  ./doall.sh
  ./doall.sh
And compare the output of the different steps (DFT, GW, HSE) by:
./plotall.sh


In any case, one can consider the <tt>doall.sh</tt> script to be an overview of the steps described below.
In any case, one can consider the <tt>doall.sh</tt> script to be an overview of the steps described below.
Line 21: Line 18:
*{{TAG|INCAR}} (see INCAR.DFT)
*{{TAG|INCAR}} (see INCAR.DFT)
   
   
  {{TAGBL|SYSTEM}}  = SrVO3                       # system name
  {{TAGBL|SYSTEM}}  = SrVO3   # system name
  {{TAGBL|NBANDS}} = 36                           # small number of bands
  {{TAGBL|NBANDS}} = 36       # small number of bands
  {{TAGBL|ISMEAR}} = 0                             # Gaussian smearing
  {{TAGBL|ISMEAR}} = -1        # Fermi smearing
  {{TAGBL|EDIFF}} = 1E-8                           # high precision for groundstate calculation
{{TAGBL|SIGMA}} = 0.1        # electronic temperature in eV (1eV ~ 11604K)
  {{TAGBL|KPAR}} = 2                               # parallelization of k-points in two groups
  {{TAGBL|EDIFF}} = 1E-8       # high precision for groundstate calculation
  {{TAGBL|KPAR}} = 2           # parallelization of k-points in two groups


Copy the aforementioned file to {{TAG|INCAR}}:
Copy the aforementioned file to {{TAG|INCAR}}:
Line 52: Line 50:
The {{TAG|KPOINTS}} file describes how the first Brillouin zone is sampled.
The {{TAG|KPOINTS}} file describes how the first Brillouin zone is sampled.
In the first step we use a uniform k-point sampling:
In the first step we use a uniform k-point sampling:
*{{TAG|KPOINTS}} (see KPOINTS.BULK)
*{{TAG|KPOINTS}}
<pre>
<pre>
Automatically generated mesh
Automatically generated mesh
      0
0
Gamma
Gamma
  4 4 4
  4 4 4
Line 62: Line 60:


'''Mind''': this is definitely not dense enough for a high-quality description of SrVO<sub>3</sub>, but in the interest of speed we will live with it.
'''Mind''': this is definitely not dense enough for a high-quality description of SrVO<sub>3</sub>, but in the interest of speed we will live with it.
Copy the aforementioned file to {{TAG|KPOINTS}}:
cp KPOINTS.BULK KPOINTS


and run VASP. If all went well, one should obtain a {{TAG|WAVECAR}} file containing the PBE wavefunction.
Run VASP. If all went well, one should obtain a {{TAG|WAVECAR}} file containing the PBE wavefunction.


== Obtain DFT virtual orbitals and long-wave limit ==
== Obtain DFT virtual orbitals and long-wave limit ==
Use following {{TAG|INCAR}} file to increase the number of virtual states and to determine the long-wave limit of the polarizability (stored in {{TAG|WAVEDER}}):
Use following {{TAG|INCAR}} file to increase the number of virtual states and to determine the long-wave limit of the polarizability (stored in {{TAG|WAVEDER}}):


*{{TAG|INCAR}} (see INCAR.DIAG)
*{{TAG|INCAR}} (see INCAR.PBE)
   
   
  {{TAGBL|SYSTEM}} = SrVO3                         # system name
  {{TAGBL|SYSTEM}} = SrVO3         # system name
  {{TAGBL|ISMEAR}} = 0                             # Gaussian smearing
  {{TAGBL|ISMEAR}} = -1            # Fermi smearing
  {{TAGBL|KPAR}} = 2                               # parallelization of k-points in two groups
{{TAGBL|SIGMA}} = 0.1            # electronic temperature in eV (1eV ~ 11604K)
  {{TAGBL|ALGO}} = Exact                           # exact diagonalization
  {{TAGBL|KPAR}} = 2               # parallelization of k-points in two groups
  {{TAGBL|NELM}} = 1                               # one electronic step suffices, since WAVECAR from previous step is present
  {{TAGBL|ALGO}} = Exact           # exact diagonalization
  {{TAGBL|NBANDS}} = 96                           # need for a lot of bands in GW
  {{TAGBL|NELM}} = 1               # one electronic step suffices, since WAVECAR from previous step is present
  {{TAGBL|LOPTICS}} = .TRUE.                       # we need d phi/ d k  for GW calculations for long-wave limit
  {{TAGBL|NBANDS}} = 96             # need for a lot of bands in GW
  {{TAGBL|LOPTICS}} = .TRUE.       # we need d phi/ d k  for GW calculations for long-wave limit
{{TAGBL|LFINITE_TEMPERATURE}} = T # compute all optical matrix elements (only required for CRPAR)


Restart VASP.
Restart VASP.
At this stage it is a good idea to make a safety copy of the {{TAG|WAVECAR}} and {{TAG|WAVEDER}} files since we will repeatedly need them in the calculations that follow:  
At this stage it is a good idea to make a safety copy of the {{TAG|WAVECAR}} and {{TAG|WAVEDER}} files since we will repeatedly need them in the calculations that follow:  
  cp {{TAGBL|WAVECAR}} WAVECAR.DIAG
  cp {{TAGBL|WAVECAR}} WAVECAR.PBE
  cp {{TAGBL|WAVEDER}} WAVEDER.DIAG
  cp {{TAGBL|WAVEDER}} WAVEDER.PBE


== GW Step ==
== cRPA Calculation ==
The actual GW calculation requires a set of one-electron energies and eigenstates. In this case we use the PBE solution obtained from previous step:
Calculate the cRPA interaction parameters for the t2g states by using the PBE wavefunction as input
  cp WAVECAR.DIAG {{TAGBL|WAVECAR}}
  cp WAVECAR.PBE {{TAGBL|WAVECAR}}
  cp WAVEDER.DIAG {{TAGBL|WAVEDER}}
  cp WAVEDER.PBE {{TAGBL|WAVEDER}}


The following {{TAG|INCAR}} file selects the 'single shot' GW calculation also known as G<sub>0</sub>W<sub>0</sub>:
And use following input file as
*{{TAG|INCAR}} (see INCAR.GW0)
*{{TAG|INCAR}} (see INCAR.CRPA) and run vasp
   
   
  {{TAGBL|SYSTEM}} = SrVO3                         # system name
  {{TAGBL|SYSTEM}} = SrVO3           # system name
  {{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
  {{TAGBL|ISMEAR}} = -1              # Fermi smearing
  {{TAGBL|KPAR}} = 2                              # parallelization of k-points in two groups
  {{TAGBL|SIGMA}} = 0.1              # electronic temperature in eV (1eV ~ 11604K)
  {{TAGBL|ALGO}} = GW0                            # GW with iteration in G, W kept on DFT level
  {{TAGBL|NCSHMEM}} = 1              # switch off shared memory for chi
  {{TAGBL|NELM}} = 1                              # one electronic step suffices, since WAVECAR from previous step is present
  {{TAGBL|ALGO}} = CRPA              # Switch on CRPA
  {{TAGBL|NBANDS}} = 96                           # need for a lot of bands in GW
  {{TAGBL|NBANDS}} = 96               # CRPA needs many empty states
  {{TAGBL|PRECFOCK}} = Fast                       # fast mode for FFTs
  {{TAGBL|PRECFOCK}} = Fast           # fast mode for FFTs
  {{TAGBL|ENCUTGW}} = 100                          # small energy cutoff for response function suffices for this tutorial
  {{TAGBL|NTARGET_STATES}} = 1 2 3    # exclude wannier states 1 - 3 in screening
  {{TAGBL|NOMEGA}} = 200                          # large number of real frequency points for Hilbert transforms of W and self-energy
  {{TAGBL|LWRITE_WANPROJ}} = .TRUE.  # write wannier projection file
 
#{{TAGBL|LFINITE_TEMPERATURE}} = T  # T>0 formalism can be avoided here, because U computed only at omega=0
Restarting VASP will overwrite the present {{TAGBL|WAVECAR}} and {{TAGBL|vasprun.xml}} file. Make a copy them for later.  
{{NB|warning|As of version 6.2.0 run this tutorial successfully by adding following lines to INCAR. For older version copy wannier90.win.CRPA to wannier90.win and omit following lines in INCAR.}}
cp {{TAGBL|WAVECAR}} WAVECAR.GW0
  {{TAGBL|NUM_WANN}} = 3
  cp {{TAGBL|vasprun.xml}} vasprun.GW0.xml
  {{TAGBL|WANNIER90_WIN}} = "
 
  num_bands=  96
=== The dielectric function ===
 
To extract the frequency dependent dielectric constant, both in the independent-particle picture as well as including local field effects (either in DFT or in the RPA) and plot the real and imaginary components using ''gnuplot'', execute
 
  ./plotchi
 
== HSE hybrid functional ==
 
To illustrate the kind of results one would obtain for SrVO<sub>3</sub> using the [[Hartree-Fock_and_HF/DFT_hybrid_functionals#range_separated|DFT/Hartree-Fock hybrid functional HSE]], without actually doing a full selfconsistent calculation, we will recalculate the one-electron energies and DOS ({{TAG|ALGO}}=Eigenval) using the HSE functional with DFT orbitals as input
  cp WAVECAR.DIAG {{TAGBL|WAVECAR}}
 
Use the following {{TAG|INCAR}} file:
*{{TAG|INCAR}} (see INCAR.HSE)
   
   
  {{TAGBL|SYSTEM}} = SrVO3                        # system name
  # because bands 21, 22, 23 do not cross with other bands
  {{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
  # one can exclude all other bands in wannierization
  {{TAGBL|KPAR}} = 2                              # parallelization of k-points in two groups
  # and omit the definition of an energy window like so
  {{TAGBL|ALGO}} = Eigenval                        # calulate eigenvalues
  exclude_bands = 1-20, 24-96
{{TAGBL|NELM}} = 1                               # one electronic step suffices, since WAVECAR from previous step is present
{{TAGBL|NBANDS}} = 48                            # small number of bands suffice
{{TAGBL|PRECFOCK}} = Fast                        # fast mode for FFTs
{{TAGBL|LHFCALC}} = .TRUE.                      # switch on Hartree-Fock routines to calculate exact exchange
{{TAGBL|HFSCREEN}} = 0.2                        # HSE06 screening parameter


Restart VASP and make a copy of the wavefunction for post-processing
  begin projections
  cp {{TAGBL|WAVECAR}} WAVECAR.HSE
  V:dxy;dxz;dyz
 
end projections
== Post-processing: Density of states and Band structure for PBE, GW and HSE ==
  "
 
The cRPA interaction values for <math>\omega=0</math> can be found in the {{TAG|OUTCAR}}:
=== Density of States ===
spin components:  1  1, frequency:    0.0000    0.0000
The DOS of the PBE, GW and HSE solution can be calculated in a post-processing step with
*{{TAG|INCAR}} (see INCAR.DOS)
   
   
  {{TAGBL|SYSTEM}} = SrVO3                        # system name
  screened Coulomb repulsion U_iijj between MLWFs:
{{TAGBL|ISMEAR}} = -5                            # Bloechl's tetrahedron method (requires at least 3x3x3 k-points)
        1         2        3
{{TAGBL|ALGO}} = NONE                            # no electronic changes required
    1   3.3459    2.3455    2.3455
{{TAGBL|NELM}} = 1                               # one electronic step suffices, since WAVECAR from previous step is present
    2    2.3455    3.3459    2.3455
{{TAGBL|NBANDS}} = 48                            # number of bands used
    3    2.3455    2.3455    3.3459
{{TAGBL|EMIN}} = -20 ; {{TAGBL|EMAX}} = 20                # smallest/largest energy included in calculation
{{TAGBL|NEDOS}} = 1000                          # sampling points for DOS
{{TAGBL|LORBIT}} = 11                            # calculate l-m decomposed DOS
{{TAGBL|LWAVE}} = .FALSE.                        # do not overwrite WAVECAR
{{TAGBL|LCHARG}} = .FALSE.                      # do not overwrite CHGCAR
 
and requires the apropriate {{TAG|WAVECAR}} file from one of the previous steps. Copy
cp WAVECAR.DIAG WAVECAR
or
cp WAVECAR.GW0 WAVECAR
or
cp WAVECAR.HSE WAVECAR
and restart VASP. The density of states is written to {{TAG|DOSCAR}}, make a copy of this file
cp {{TAGBL|DOSCAR}} DOSCAR.XXX
where XXX is either PBE, GW0 or HSE. Visualize the projected DOS for the V-t2g, V-eg and O-p states with the scriptfile
./plotdos.sh DOSCAR.*
This requires gnuplot to be installed.
 
===  Band structure with <tt>wannier90</tt>===
The band structure can be calculated via Wannier interpolation using <tt>wannier90</tt> in the library mode
*{{TAG|INCAR}} (see INCAR.WAN)
{{TAGBL|SYSTEM}} = SrVO3                        # system name
{{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
{{TAGBL|ALGO}} = NONE                            # no electronic changes required
{{TAGBL|NELM}} = 1                               # one electronic step suffices, since WAVECAR from previous step is present
{{TAGBL|NBANDS}} = 48                            # number of bands used
{{TAGBL|LWAVE}} = .FALSE.                       # do not overwrite WAVECAR
{{TAGBL|LCHARG}} = .FALSE.                       # do not overwrite CHGCAR
{{TAGBL|LWANNIER90_RUN}} = .TRUE.               # run wannier90 in library mode
Use the corresponding wannier90.win.XXX file as input for <tt>wannier90</tt>
cp wannier90.win.XXX wannier90.win
where XXX=PBE, GW0 or HSE and looks similar to
bands_plot = true
   
   
  begin kpoint_path
  screened Coulomb repulsion U_ijji between MLWFs:
R  0.50000000  0.50000000  0.50000000  G  0.00000000  0.00000000  0.00000000
        1        2        3
G  0.00000000  0.00000000  0.00000000  X  0.50000000  0.00000000  0.00000000
    1    3.3459    0.4281    0.4281
X  0.50000000  0.00000000  0.00000000  M  0.50000000  0.50000000  0.00000000
    2    0.4281    3.3459    0.4281
M  0.50000000  0.50000000  0.00000000  G  0.00000000  0.00000000  0.00000000
    3    0.4281    0.4281    3.3459
end kpoint_path
   
   
  # number of wannier states
  screened Coulomb repulsion U_ijij between MLWFs:
num_wann =   3
        1        2        3
    1    3.3459    0.4281    0.4281
    2    0.4281    3.3459    0.4281
    3    0.4281    0.4281   3.3459
   
   
  # number of bloch bands
  averaged interaction parameter
  num_bands=  96
  screened Hubbard U =   3.3459    0.0000
screened Hubbard u =    2.3455    0.0000
screened Hubbard J =    0.4281  -0.0000
The full interaction matrix is written to {{FILE|UIJKL}}.
{{NB|mind|The frequency point <math>\omega</math> can be set by {{TAG|OMEGAMAX}} in the INCAR.}}
For instance to evaluate the cRPA interaction matrix at <math>\omega=10</math> eV, add
  {{TAGBL|OMEGAMAX}} = 10
to the INCAR and restart VASP. In contrast, adding following two lines to the {{TAG|INCAR}}
  {{TAGBL|OMEGAMAX}} = 10
  {{TAGBL|NOMEGAR}} = 0
tells VASP to calculate the interaction on the imaginary frequency axis at <math>\omega=i 10</math>. This can be used to evaluate <math>U</math> at a specific Matsubara frequency point.
 
In addition, the bare Coulomb interaction matrix is calculated for a high {{TAG|VCUTOFF}} and low energy cutoff {{TAG|ENCUTGW}} and written in that order to the {{TAG|OUTCAR}} file. Look for the lines similar to:
spin components:  1  1
 
bare Coulomb repulsion V_iijj between MLWFs:
        1        2        3
    1  16.3485  15.0984  15.0984
    2  15.0984  16.3485  15.0984
    3  15.0984  15.0984  16.3485
 
bare Coulomb repulsion V_ijji between MLWFs:
        1        2        3
    1  16.3485    0.5351    0.5351
    2    0.5351  16.3485    0.5351
    3    0.5351    0.5351  16.3485
 
bare Coulomb repulsion V_ijij between MLWFs:
        1        2        3
    1  16.3485    0.5351    0.5351
    2    0.5351  16.3485    0.5351
    3    0.5351    0.5351   16.3485
   
   
  # GW energy window for t2g states
  averaged bare interaction
  dis_win_min = 7.4
bare Hubbard U =  16.3485  -0.0000
  dis_win_max = 9.95
bare Hubbard u =  15.0984  -0.0000
bare Hubbard J =    0.5351    0.0000
Similar to the effectively screened interaction the full output is written to {{FILE|VIJKL}}.
 
=== cRPA calculation on Matsubara axis ===
{{NB|mind|Available as of VASP.6.5.2.}}
Note that the same frequency grid is used as for {{TAG|ALGO}}=RPA (RPA correlation energy calculation) and can not be changed directly.
To calculate the cRPA interaction for a set of automatically chosen imaginary frequency points use once again the PBE wavefunction as input
  cp WAVECAR.PBE {{TAGBL|WAVECAR}}
  cp WAVEDER.PBE {{TAGBL|WAVEDER}}
Currently, this step requires the {{TAG|WANPROJ}} file from previous step, no wannier90.win file is necessary.
You can also, delete {{TAG|WANPROJ}} and define a wannier projection as in the previous step in the {{FILE|INCAR}}.
 
Select the space-time cRPA algorithm with following file:
*{{TAG|INCAR}} (see INCAR.CRPAR)
   
   
  begin projections
  {{TAGBL|SYSTEM}} = SrVO3            # system name
  V:dxy;dxz;dyz
  {{TAGBL|LFINITE_TEMPERATURE}} = T    # use finite temperature formalism
  end projections
  {{TAGBL|ISMEAR}} = -1                # required for finite temperature algorithm
{{TAGBL|SIGMA}} = 0.1                # electron temperature in eV (1 eV ~ 11000 K)
{{TAGBL|ALGO}} = CRPAR              # Switch on CRPA on imaginary axis
{{TAGBL|NBANDS}} = 96                # CRPA needs many empty states
{{TAGBL|PRECFOCK}} = Fast            # fast mode for FFTs
{{TAGBL|NTARGET_STATES}} = 1 2 3    # exclude wannier states 1 - 3 in screening
{{TAGBL|NCRPA_BANDS}} = 21 22 23    # remove bands 21-23 in screening, currently required for space-time algo
{{TAGBL|NOMEGA}} = 8                # use 8 imaginary frequency points
{{TAGBL|NOMEGA_DUMP}} = 0            # write WFULLxxxx.tmp files at omega=0, used for off-centre Coulomb integrals in next step


Use the corresponding WAVECAR.XXX file as input
cp WAVECAR.XXX {{TAGBL|WAVECAR}}
and restart VASP. If all went well, the Vanadium t2g band dispersion thus obtained, may conveniently be visualized with gnuplot:
gnuplot -persist ./wannier90_band.gnu


:'''N.B.:''' Most modern versions of <tt>gnuplot</tt> will respond with an error message unless you remove the first line of <tt>wannier90_band.gnu</tt> (some deprecated syntax issue).
Run VASP and make a copy of the output file
cp {{TAGBL|OUTCAR}} OUTCAR.CRPAR


=== Alternative way to calculate the PBE band structure ===
After a successful run, the interaction values at {{TAGBL|NOMEGA}}+1 frequencies are written to the {{TAG|OUTCAR}} file, where the first point
VASP allows to interpolate the PBE band structure from the PBE charge density
is always <math>\omega=0</math>:
  cp CHGCAR.DIAG CHGCAR
   spin components: 1 1, frequency:    0.0000    0.0000
  cp WAVECAR.DIAG WAVECAR
by adapting the {{TAG|KPOINTS}} file as follows:
*{{TAG|KPOINTS}} (see KPOINTS.BSTR)
Auto
15
Linemode
reciprocal
0.50000000  0.50000000  0.50000000   !R
  0.00000000 0.00000000  0.00000000  !G
   
   
  0.00000000  0.00000000  0.00000000  !G
  screened Coulomb repulsion U_iijj between MLWFs:
0.50000000  0.00000000  0.00000000  !X
        1        2        3
    1    3.3450    2.3447    2.3447
    2    2.3447    3.3450    2.3447
    3    2.3447    2.3447    3.3450
   
   
  0.50000000  0.00000000  0.00000000  !X
  ...
0.50000000  0.50000000  0.00000000  !M
   
   
0.50000000  0.50000000  0.00000000  !M
   spin components:  1 1, frequency:   0.0000 109.6955
0.00000000  0.00000000  0.00000000   !G
The following {{TAGBL|INCAR}} file tells VASP to interpolate the band structure:
*{{TAG|INCAR}} (see INCAR.BSTR)
{{TAGBL|SYSTEM}} = SrVO3                        # system name
{{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
{{TAGBL|EDIFF}} = 1E-7                          # tight convergence criterion
{{TAGBL|NBANDS}} = 36                            # 36 bands are sufficient
{{TAGBL|LWAVE}} = .FALSE.                        # do not overwrite WAVECAR
  {{TAGBL|LCHARG}} = .FALSE.                      # do not overwrite CHGCAR
  {{TAGBL|ICHARG}} = 11                            # use CHGCAR file for interpolation
{{TAGBL|LORBIT}} = 11                            # compute lm-decomposed states
{{TAGBL|EMIN}} = -20 ; {{TAGBL|EMAX}} = 20                # smallest/largest energy included in calculation
{{TAGBL|NEDOS}} = 1000                          # sampling points for DOS
This PBE band structure and the Wannier-interpolated structures of the PBE, HSE and GW calculation can be compared via
./plotbands.sh
:'''N.B.:''' Mind that this approach works only for DFT wavefunctions, like PBE or LDA.
 
== CRPA Calculation ==
Calculate the CRPA interaction parameters for the t2g states by using the PBE wavefunction as input
  cp WAVECAR.DIAG {{TAGBL|WAVECAR}}
cp WAVEDER.DIAG {{TAGBL|WAVEDER}}
Use following Wannier projection for the basis:
 
*{{TAG|wannier90.win}} (see wannier90.win.CRPA)
num_wann =    3
num_bands=  96
   
   
  # PBE energy window of t2g states (band 21-23)
  screened Coulomb repulsion U_iijj between MLWFs:
dis_win_min = 6.4
        1        2        3
dis_win_max = 9.0
    1  15.2510  14.0759   14.0759
    2  14.0759  15.2510  14.0759
begin projections
    3  14.0759  14.0759  15.2510
   V:dxy;dxz;dyz
end projections
 
Copy this file to wannier90.win
cp wannier90.win.CRPA wannier90.win
 
And use following input file as
*{{TAG|INCAR}} (see INCAR.CRPA)
{{TAGBL|SYSTEM}} = SrVO3                        # system name
{{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
{{TAGBL|NCSHMEM}} = 1                            # switch off shared memory for chi
{{TAGBL|ALGO}} = CRPA                            # Switch on CRPA
{{TAGBL|NBANDS}} = 96                            # CRPA needs many empty states
{{TAGBL|PRECFOCK}} = Fast                        # fast mode for FFTs
{{TAGBL|NTARGET_STATES}} = 1 2 3                # exclude wannier states 1 - 3 in screening
{{TAGBL|LWRITE_WANPROJ}} = .TRUE.               # write wannier projection file


and run VASP. The CRPA interaction values can be found in the {{TAG|OUTCAR}} file after following lines
The complete matrix at zero frequency is also written to {{FILE|UIJKL}}, while the result at the first frequency point of the minimax grid{{cite|Kaltak:PRB:2020}} is found in {{FILE|UIJKL}}.1 and so on. 
screened Coulomb repulsion U_iijj between MLWFs:
including an averaged value:
screened Hubbard U =    3.3746  -0.0000
Make a copy of the output file
cp {{TAGBL|OUTCAR}} OUTCAR.CRPA


=== CRPA calculation on full imaginary frequency axis (optional) ===
==== Optional: Analytic continuation ====
To calculate the CRPA interaction for a set of imaginary frequency points use once again the PBE wavefunction as input
To obtain the effective interaction on the real frequency axis from the imaginary axis (stored in {{FILE|UIJKL}}.*) following python code can be used in a [https://jupyter.org/ jupyter notebook]:
  cp WAVECAR.DIAG {{TAGBL|WAVECAR}}
<syntaxhighlight lang="python" line>
cp WAVEDER.DIAG {{TAGBL|WAVEDER}}
import numpy as np
This step requires uses the WANPROJ file from previous step, no wannier90.win file is necessary.
from scipy.interpolate import AAA
import matplotlib.pyplot as plt
# results stored in NOMEGA dimensional array
nomega=24
u = np.zeros( nomega, dtype='complex' ) 
# one-site indices
idx=[0,40,80 ]
# store one-site bare interaction
v = np.sum( np.loadtxt('VIJKL').T[4][idx] )/len(idx)
for i in range(nomega):
    raw = np.loadtxt( 'UIJKL.{omega:d}'.format(omega=i+1) )
    u[i] = np.sum( raw.T[4][idx] + 1j * raw.T[5][idx] ) / len(idx)
# extract omega points
!grep "omega =" UIJKL.{?,??} | awk '{print $5}' > omegas.dat
omegas=np.loadtxt( 'omegas.dat' )*1j
# use AAA algorithm for analytic continuation
u_cont = AAA(omegas,u )
# plot real part of U and bare interaction V
z = np.linspace( 0, 200, num=1000)
fig, ax = plt.subplots()
ax.plot( z, u_cont(z).real, '-', color='r', label='U')
ax.axhline(y=v, color='b', linestyle='-', label='V')
ax.legend()
# add low-frequency regime as an inset
zlow=np.linspace( 0, 20, num=1000)
inset_ax = ax.inset_axes([0.35, 0.1, 0.6, 0.6]) # [left, bottom, width, height]
inset_ax.plot(zlow, u_cont(zlow), color='red')
plt.xlabel('$\omega$ [eV]')
plt.show()
</syntaxhighlight>
[[File:SrVO3 U omega acont.png|480px]]
{{NB|tip|Increase {{TAG|NOMEGA}} points to resolve more details on the real frequency axis.}}


Select the space-time CRPA algorithm with following file:
== Off-centre Coulomb integrals ==
*{{TAG|INCAR}} (see INCAR.CRPAR)
Every cRPA job writes the effectively screened Coulomb kernel (in reciprocal space) at zero frequency to {{FILE|WFULLxxxx.tmp}} files.  
These files can be read in and off-centre Coulomb integrals can be evaluated using following input:
  {{TAGBL|SYSTEM}} = SrVO3                        # system name
  {{TAGBL|ALGO}} = 2E4WA                          # Compute off-centre Coulomb integrals
  {{TAGBL|ISMEAR}} = 0                            # Gaussian smearing
  {{TAGBL|ISMEAR}} = -1                            # Fermi smearing            
  {{TAGBL|NCSHMEM}} = 1                            # switch off shared memory for chi
  {{TAGBL|SIGMA}} = 0.1                            # electronic temperature
  {{TAGBL|ALGO}} = CRPAR                          # Switch on CRPA on imaginary axis
  {{TAGBL|NBANDS}} = 96                            # use same number of bands as stored in WAVECAR
  {{TAGBL|NBANDS}} = 96                            # CRPA needs many empty states
  {{TAGBL|ALGO}} = 2E4WA                          # Compute off-centre Coulomb integrals
  {{TAGBL|PRECFOCK}} = Fast                        # fast mode for FFTs
  {{TAGBL|PRECFOCK}} = Fast                        # fast mode for FFTs
  {{TAGBL|NTARGET_STATES}} = 1 2 3                # exclude wannier states 1 - 3 in screening
  {{TAGBL|NTARGET_STATES}} = 1 2 3                # Wannier states for which Coulomb integrals are evaluated
{{TAGBL|NCRPA_BANDS}} = 21 22 23                # remove bands 21-23 in screening, currently required for space-time algo
{{TAGBL|NOMEGA}} = 12                            # use 12 imaginary frequency points
{{TAGBL|NTAUPAR}} = 4                            # distribute 12 time points into 4 groups


Run VASP and make a copy of the output file
The bare off-centre Coulomb integrals are written to {{FILE|VRijkl}}, while the effectively screened ones are found in {{FILE|URijkl}}:
cp {{TAGBL|OUTCAR}} OUTCAR.CRPAR
# File generated by VASP contains Coulomb matrix elements
The resulting interactions are written for every imaginary frequency point to the {{TAG|OUTCAR}} file.
  # U_ijkl = [ij|kl]
For instance, to extract the averaged on-site U interaction for each point enter following command
#  I  J  K  L          RE(U_IJKL)          IM(U_IJKL)
  grep "screened Hubbard U" OUTCAR
  # R:   1  0.000000  0.000000  0.000000
resulting in following output
    1  1  1  1        3.3450226866        0.0000000000
  screened Hubbard U =   3.3798  -0.0000
    2  1  1  1        0.0000058776        0.0000006717
screened Hubbard U =    3.4172  -0.0000
    3   1  1  1        0.0000026927      -0.0000003292
screened Hubbard U =    3.5169  -0.0000
  ...
screened Hubbard U =    3.7418  -0.0000
    3  3  3   3        3.3450230976        0.0000000000
  screened Hubbard U =    4.2069   -0.0000
  # R:   2  0.000000 0.000000 1.000000
  screened Hubbard U =   5.0802  -0.0000
    1  1   1  1        0.7321605022      -0.0001554484
  screened Hubbard U =    6.5456  -0.0000
    2  1   1   1        0.0000021144        0.0000002295
  screened Hubbard U =    8.6426   -0.0000
  ...
screened Hubbard U =  11.0815  -0.0000
{{NB|mind|Available as of VASP.6.5.2.}}
screened Hubbard U =   13.3615   -0.0000
screened Hubbard U =  15.0636  -0.0000
  screened Hubbard U =  16.0412  -0.0000
Here each line corresponds to an (increasing) imaginary frequency point.  
The first line is the CRPA interaction at the lowest frequency point and is roughly the same as the value at 0 calculated in previous step.
The last line (interaction at the highest frequency point) approaches the bare Coulomb interaction in this basis, which is also written to the {{TAG|OUTCAR}}:
bare Hubbard U =  16.3485    0.0000


== Downloads ==
[[Media:CRPA of SrVO3.zip| CRPA_of_SrVO3.zip]]


{{Template:gw}}
{{Template:GW - Tutorial}}


Back to the [[The_VASP_Manual|main page]].
Back to the [[The_VASP_Manual|main page]].


[[Category:Examples]]
[[Category:Examples]][[Category:Constrained-random-phase approximation]]

Latest revision as of 09:10, 28 March 2025

The following tutorial describes how to perform cRPA calculations, which is available as of VASP 6.

Task

Calculation of the Coulomb matrix elements in the constrained Random Phase Approximation (cRPA) of SrVO3 between the Vanadium t2g states.

Performing a cRPA calculation with VASP is a 3-step procedure: a DFT groundstate calculation, a calculation to obtain a number of virtual orbitals, and the actual cRPA calculation itself.

N.B.: This example involves quite a number of individual calculations. The easiest way to run this example is to execute: 

./doall.sh

In any case, one can consider the doall.sh script to be an overview of the steps described below.

DFT groundstate calculation

The first step is a conventional DFT (in this case PBE) groundstate calculation. 

SYSTEM  = SrVO3    # system name
NBANDS = 36        # small number of bands
ISMEAR = -1        # Fermi smearing
SIGMA = 0.1        # electronic temperature in eV (1eV ~ 11604K)
EDIFF = 1E-8       # high precision for groundstate calculation
KPAR = 2           # parallelization of k-points in two groups

Copy the aforementioned file to INCAR: 

cp INCAR.DFT INCAR

The POSCAR file describes the structure of the system: 

SrVO3
3.84652  #cubic fit for 6x6x6 k-points
 +1.0000000000  +0.0000000000  +0.0000000000 
 +0.0000000000  +1.0000000000  +0.0000000000 
 +0.0000000000  +0.0000000000  +1.0000000000 
Sr V O
 1 1 3
Direct
 +0.0000000000  +0.0000000000  +0.0000000000 
 +0.5000000000  +0.5000000000  +0.5000000000 
 +0.5000000000  +0.5000000000  +0.0000000000 
 +0.5000000000  +0.0000000000  +0.5000000000 
 +0.0000000000  +0.5000000000  +0.5000000000

This file remains unchanged in the following.

The KPOINTS file describes how the first Brillouin zone is sampled. In the first step we use a uniform k-point sampling: 

Automatically generated mesh
 0
Gamma
 4 4 4
 0 0 0

Mind: this is definitely not dense enough for a high-quality description of SrVO3, but in the interest of speed we will live with it.

Run VASP. If all went well, one should obtain a WAVECAR file containing the PBE wavefunction.

Obtain DFT virtual orbitals and long-wave limit

Use following INCAR file to increase the number of virtual states and to determine the long-wave limit of the polarizability (stored in WAVEDER): 

SYSTEM = SrVO3          # system name
ISMEAR = -1             # Fermi smearing
SIGMA = 0.1             # electronic temperature in eV (1eV ~ 11604K)
KPAR = 2                # parallelization of k-points in two groups
ALGO = Exact            # exact diagonalization
NELM = 1                # one electronic step suffices, since WAVECAR from previous step is present
NBANDS = 96             # need for a lot of bands in GW
LOPTICS = .TRUE.        # we need d phi/ d k  for GW calculations for long-wave limit
LFINITE_TEMPERATURE = T # compute all optical matrix elements (only required for CRPAR)

Restart VASP. At this stage it is a good idea to make a safety copy of the WAVECAR and WAVEDER files since we will repeatedly need them in the calculations that follow: 

cp WAVECAR WAVECAR.PBE
cp WAVEDER WAVEDER.PBE

cRPA Calculation

Calculate the cRPA interaction parameters for the t2g states by using the PBE wavefunction as input 

cp WAVECAR.PBE WAVECAR
cp WAVEDER.PBE WAVEDER

And use following input file as

  • INCAR (see INCAR.CRPA) and run vasp
SYSTEM = SrVO3            # system name
ISMEAR = -1               # Fermi smearing
SIGMA = 0.1               # electronic temperature in eV (1eV ~ 11604K)
NCSHMEM = 1               # switch off shared memory for chi
ALGO = CRPA               # Switch on CRPA
NBANDS = 96               # CRPA needs many empty states
PRECFOCK = Fast           # fast mode for FFTs
NTARGET_STATES = 1 2 3    # exclude wannier states 1 - 3 in screening
LWRITE_WANPROJ = .TRUE.   # write wannier projection file
#LFINITE_TEMPERATURE = T  # T>0 formalism can be avoided here, because U computed only at omega=0
Warning: As of version 6.2.0 run this tutorial successfully by adding following lines to INCAR. For older version copy wannier90.win.CRPA to wannier90.win and omit following lines in INCAR. 
NUM_WANN = 3 
WANNIER90_WIN = "
num_bands=   96

# because bands 21, 22, 23 do not cross with other bands
# one can exclude all other bands in wannierization
# and omit the definition of an energy window like so
exclude_bands = 1-20, 24-96

begin projections
 V:dxy;dxz;dyz
end projections
"

The cRPA interaction values for can be found in the OUTCAR: 

spin components:  1  1, frequency:    0.0000    0.0000

screened Coulomb repulsion U_iijj between MLWFs:
        1         2         3
   1    3.3459    2.3455    2.3455
   2    2.3455    3.3459    2.3455
   3    2.3455    2.3455    3.3459

screened Coulomb repulsion U_ijji between MLWFs:
        1         2         3
   1    3.3459    0.4281    0.4281
   2    0.4281    3.3459    0.4281
   3    0.4281    0.4281    3.3459

screened Coulomb repulsion U_ijij between MLWFs:
        1         2         3
   1    3.3459    0.4281    0.4281
   2    0.4281    3.3459    0.4281
   3    0.4281    0.4281    3.3459

averaged interaction parameter
screened Hubbard U =    3.3459    0.0000
screened Hubbard u =    2.3455    0.0000
screened Hubbard J =    0.4281   -0.0000

The full interaction matrix is written to UIJKL.

Mind: The frequency point can be set by OMEGAMAX in the INCAR.

For instance to evaluate the cRPA interaction matrix at eV, add 

 OMEGAMAX = 10

to the INCAR and restart VASP. In contrast, adding following two lines to the INCAR 

 OMEGAMAX = 10 
 NOMEGAR = 0

tells VASP to calculate the interaction on the imaginary frequency axis at . This can be used to evaluate at a specific Matsubara frequency point.

In addition, the bare Coulomb interaction matrix is calculated for a high VCUTOFF and low energy cutoff ENCUTGW and written in that order to the OUTCAR file. Look for the lines similar to: 

spin components:  1  1
 
bare Coulomb repulsion V_iijj between MLWFs:
        1         2         3
   1   16.3485   15.0984   15.0984
   2   15.0984   16.3485   15.0984
   3   15.0984   15.0984   16.3485
 
bare Coulomb repulsion V_ijji between MLWFs:
        1         2         3
   1   16.3485    0.5351    0.5351
   2    0.5351   16.3485    0.5351
   3    0.5351    0.5351   16.3485
 
bare Coulomb repulsion V_ijij between MLWFs:
        1         2         3
   1   16.3485    0.5351    0.5351
   2    0.5351   16.3485    0.5351
   3    0.5351    0.5351   16.3485

averaged bare interaction
bare Hubbard U =   16.3485   -0.0000
bare Hubbard u =   15.0984   -0.0000
bare Hubbard J =    0.5351    0.0000

Similar to the effectively screened interaction the full output is written to VIJKL.

cRPA calculation on Matsubara axis

Mind: Available as of VASP.6.5.2.

Note that the same frequency grid is used as for ALGO=RPA (RPA correlation energy calculation) and can not be changed directly. To calculate the cRPA interaction for a set of automatically chosen imaginary frequency points use once again the PBE wavefunction as input 

cp WAVECAR.PBE WAVECAR
cp WAVEDER.PBE WAVEDER

Currently, this step requires the WANPROJ file from previous step, no wannier90.win file is necessary. You can also, delete WANPROJ and define a wannier projection as in the previous step in the INCAR.

Select the space-time cRPA algorithm with following file: 

SYSTEM = SrVO3             # system name
LFINITE_TEMPERATURE = T    # use finite temperature formalism 
ISMEAR = -1                # required for finite temperature algorithm
SIGMA = 0.1                # electron temperature in eV (1 eV ~ 11000 K)
ALGO = CRPAR               # Switch on CRPA on imaginary axis
NBANDS = 96                # CRPA needs many empty states
PRECFOCK = Fast            # fast mode for FFTs
NTARGET_STATES = 1 2 3     # exclude wannier states 1 - 3 in screening
NCRPA_BANDS = 21 22 23     # remove bands 21-23 in screening, currently required for space-time algo
NOMEGA = 8                 # use 8 imaginary frequency points
NOMEGA_DUMP = 0            # write WFULLxxxx.tmp files at omega=0, used for off-centre Coulomb integrals in next step


Run VASP and make a copy of the output file 

cp OUTCAR OUTCAR.CRPAR

After a successful run, the interaction values at NOMEGA+1 frequencies are written to the OUTCAR file, where the first point is always : 

 spin components:  1  1, frequency:    0.0000    0.0000

screened Coulomb repulsion U_iijj between MLWFs:
        1         2         3
   1    3.3450    2.3447    2.3447
   2    2.3447    3.3450    2.3447
   3    2.3447    2.3447    3.3450

...

 spin components:  1  1, frequency:    0.0000  109.6955

screened Coulomb repulsion U_iijj between MLWFs:
        1         2         3
   1   15.2510   14.0759   14.0759
   2   14.0759   15.2510   14.0759
   3   14.0759   14.0759   15.2510

The complete matrix at zero frequency is also written to UIJKL, while the result at the first frequency point of the minimax grid[1] is found in UIJKL.1 and so on.

Optional: Analytic continuation

To obtain the effective interaction on the real frequency axis from the imaginary axis (stored in UIJKL.*) following python code can be used in a jupyter notebook: 

import numpy as np 
from scipy.interpolate import AAA
import matplotlib.pyplot as plt
# results stored in NOMEGA dimensional array 
nomega=24
u = np.zeros( nomega, dtype='complex' )  
# one-site indices 
idx=[0,40,80 ]
# store one-site bare interaction
v = np.sum( np.loadtxt('VIJKL').T[4][idx] )/len(idx)
for i in range(nomega):
    raw = np.loadtxt( 'UIJKL.{omega:d}'.format(omega=i+1) ) 
    u[i] = np.sum( raw.T[4][idx] + 1j * raw.T[5][idx] ) / len(idx)
# extract omega points 
!grep "omega =" UIJKL.{?,??} | awk '{print $5}' > omegas.dat
omegas=np.loadtxt( 'omegas.dat' )*1j
# use AAA algorithm for analytic continuation 
u_cont = AAA(omegas,u )
# plot real part of U and bare interaction V 
z = np.linspace( 0, 200, num=1000)
fig, ax = plt.subplots()
ax.plot( z, u_cont(z).real, '-', color='r', label='U')
ax.axhline(y=v, color='b', linestyle='-', label='V')
ax.legend()
# add low-frequency regime as an inset 
zlow=np.linspace( 0, 20, num=1000)
inset_ax = ax.inset_axes([0.35, 0.1, 0.6, 0.6])  # [left, bottom, width, height]
inset_ax.plot(zlow, u_cont(zlow), color='red')
plt.xlabel('$\omega$ [eV]')
plt.show()

Tip: Increase NOMEGA points to resolve more details on the real frequency axis.

Off-centre Coulomb integrals

Every cRPA job writes the effectively screened Coulomb kernel (in reciprocal space) at zero frequency to WFULLxxxx.tmp files. These files can be read in and off-centre Coulomb integrals can be evaluated using following input: 

ALGO = 2E4WA                           # Compute off-centre Coulomb integrals
ISMEAR = -1                            # Fermi smearing             
SIGMA = 0.1                            # electronic temperature
NBANDS = 96                            # use same number of bands as stored in WAVECAR
ALGO = 2E4WA                           # Compute off-centre Coulomb integrals
PRECFOCK = Fast                        # fast mode for FFTs
NTARGET_STATES = 1 2 3                 # Wannier states for which Coulomb integrals are evaluated

The bare off-centre Coulomb integrals are written to VRijkl, while the effectively screened ones are found in URijkl: 

# File generated by VASP contains Coulomb matrix elements
# U_ijkl = [ij|kl] 
#  I   J   K   L          RE(U_IJKL)          IM(U_IJKL)
# R:    1  0.000000  0.000000  0.000000
   1   1   1   1        3.3450226866        0.0000000000
   2   1   1   1        0.0000058776        0.0000006717
   3   1   1   1        0.0000026927       -0.0000003292
...
   3   3   3   3        3.3450230976        0.0000000000
# R:    2  0.000000  0.000000  1.000000
   1   1   1   1        0.7321605022       -0.0001554484
   2   1   1   1        0.0000021144        0.0000002295
...
Mind: Available as of VASP.6.5.2.

Downloads

CRPA_of_SrVO3.zip

Back to the main page.

  1. ↑ M. Kaltak and G. Kresse, Phys. Rev. B. 101, 205145 (2020).
Retrieved from "https://www.vasp.at/wiki/index.php?title=CRPA_of_SrVO3&oldid=30278"