Part 1: Constrained random-phase approximation ¶

Content¶

 1 Band structure
 2 Wannier projection
 3 Constrained random-phase approximation

1 Band structure

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the ground state wavefunction.
• Generate the unoccupied states in preparation for a constrained random-phase approximation (cRPA) calculation.
• Project the density of states (including fat bands).

1.1 Task¶

Calculate the ground state wavefunction and unoccupied states of primitive NiO, projecting the DOS onto the p- and d- states using a 12x12x12 k-mesh.

Strongly correlated materials are systems in which electron-electron interactions play a dominant role and cannot be adequately described by independent-particle approximations such as standard DFT. These systems typically include elements with partially filled $d$ and $f$ electron orbitals which are localized. Correlation effects lead to phenomena such as metal-insulator transitions, magnetism, and unconventional superconductivity. To model such systems, several extensions of DFT have been developed: constrained random-phase approximation (cRPA), DFT+U, and dynamical mean-field theory (DMFT). We will cover cRPA in part 1 of the tutorials and DFT+U and DMFT in part 2.

cRPA is a method that enables calculating the effective interaction parameter $U$, $J$, and $U'$ for model Hamiltonians [Phys. Rev. B 77 (2008) 085122, Phys. Rev. B 112 (2025) 245102]. The main idea is to neglect screening effects of specific target states in the screened Coulomb interaction $W$ of the $GW$ method. The resulting partially screened Coulomb interaction is usually evaluated in a localized basis that spans the target space and is described by the model Hamiltonian. Usually, the target space is low-dimensional (up to 5 states) and therefore allows for the application of a higher-level theory, such as DMFT.

A cRPA calculation is performed in several steps:

1. A DFT ground state calculation (e01)
2. Unoccupied Kohn-Sham (KS) eigenstates are calculated (e01)
3. Projection of density of states (including fat bands) (e01)
4. Wannier projection of states (Ni d and O p in e02)
5. cRPA calculation (e03)

In this exercise, you will calculate the ground state of NiO and generate unoccupied eigenstates in preparation for analysing the band structure.

1.2 Determine the DFT ground state¶

1.2.1 Input files¶

The input files to run the DFT ground state calculation are prepared at $TUTORIALS/strongly_correlated/e01_bands. Check them out!

POSCAR:
NiO

NiO
 4.171
  0.5    0.5    0.0
  0.0    0.5    0.5
  0.5    0.0    0.5
Ni O
1   1
Direct
 0.00  0.00  0.00
 0.50  0.50  0.50

INCAR:
A DFT run using the default generalized-gradient-approximation functional of Perdew, Burke, and Ernzerhof (GGA-PBE).

The INCAR tags are standard DFT ones: SYSTEM provides the calculation name, NBANDS specifies the total number of Kohn-Sham orbitals in the calculation, ISMEAR = -1 ensures that Fermi smearing is used (necessary for the LFINITE_TEMPERATURE later), SIGMA is the smearing width, and EDIFF is the global break condition for the electronic self-consistency-loop.

SYSTEM = NiO
ISMEAR = -1 ; SIGMA = 0.1

EDIFF = 1E-8
NBANDS = 48
ALGO = NORMAL

KPOINTS:
An equally spaced (12$\times$12$\times$12) Γ-centered grid of $\bf{k}$ points.

Automatically generated mesh
  0
Gamma
 12 12 12

POTCAR:
For cRPA calculations, one should use the POTCAR files with _GW appended to their name. These potentials were constructed to reproduce the atomic scattering properties over a larger energy range than the conventional potentials. This is needed when one needs to compute a large number of unoccupied states. You can check the TITEL field in your POTCAR file.

   PAW_PBE Ni_GW 31Mar2010
   PAW_PBE O_GW_new 19Mar2012           

1.2.2 Run the calculation¶

Open a terminal, navigate to this example's directory and run VASP by entering the following:

cd $TUTORIALS/strongly_correlated/e01_*
mpirun -np 4 vasp_std

This calculation computes the DFT ground state. What is the size of the band gap?
Hint: Execute the following cell to obtain the difference between the highest-occupied-molecular-orbital (HOMO) and lowest-unoccupied-molecular-orbital (LUMO) eigenenergies and to plot the DOS using py4vasp!

The small bandgap indicates that NiO might be a metal; however, due to strongly correlated electrons, it is in fact a Mott insulator.

With this calculation, you have successfully obtained the WAVECAR file necessary for generating the virtual orbitals.

1.3 Compute additional unoccupied Kohn-Sham orbitals¶

In step two of this exercise, we will increase the number of virtual states and determine the long-wave limit of the polarizability (stored in WAVEDER). The input files are prepared in $TUTORIALS/strongly_correlated/e01_bands/unoccupied/. Check them out!

1.3.1 Input files¶

unoccupied/INCAR:
Increase NBANDS with respect to the default and perform a single exact diagonalization of the Hamiltonian (ALGO=Exact). The number in NBANDS has been increased so that we compute not only the occupied bands, but also the unoccupied, which are required for the cRPA calculation. The WAVECAR and WAVEDER files of this calculation will be the starting point for the actual cRPA calculation. LOPTICS=T calculates the frequency-dependent dielectric matrix, and LFINITE_TEMPERATURE=T is required to compute all optical matrix elements.

SYSTEM = NiO
ISMEAR = -1 ; SIGMA = 0.1

EDIFF = 1E-8
# increase number of unoccupied states
# 96 is usually not enough, but will work for this tutorial
NBANDS = 96
# exact diagonalization is required
# to obtain good unoccupied orbitals and energies
ALGO = EXACT

#compute optical matrix elements
LOPTICS = T
# also include all un-occupied optical matrix elements
# since we plan to use the output files of this step for cRPA later on
LFINITE_TEMPERATURE = T

1.3.2 Run the calculation¶

To calculate the unoccupied states, go to the unoccupied directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e01_*/unoccupied
cp ../WAVECAR .
mpirun -np 4 vasp_std

The number in NBANDS has been increased so that it is not only the occupied bands that we calculate but also the unoccupied.

The unoccupied states are necessary for performing the subsequent cRPA calculations. First, you will choose the bands of interest to define the localized basis set in terms of Wannier functions.

1.4 Projection of density of states (including fat bands)¶

In step three of this exercise, we will calculate the density of states and find the p- and d- bands of NiO. The input files are prepared in $TUTORIALS/strongly_correlated/e01_bands/fatbands/. Check them out!

1.4.1 Input files¶

fatbands/INCAR:
EMIN and EMAX define the minimum and maximum limits for the pDOS, with NEDOS determining the number of grid points.

SYSTEM = NiO
ISMEAR = -1 ; SIGMA = 0.1

EDIFF = 1E-8
NBANDS = 48
ALGO = NORMAL
# do not edit wavefunction and charge density
LWAVE = F
# freeze density to interpolate the band structure
ICHARG = 11
# create projections on d- and p- states
LORBIT = 11
# higher resolution for pDOS
EMIN = -5 ; EMAX = 15 ; NEDOS = 2000

1.4.2 Run the calculation¶

To calculate the fatbands, go to the fatbands directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e01_*/fatbands
cp ../unoccupied/{CHG,WAVE}CAR .
mpirun -np 4 vasp_std

Check which bands are most relevant to capture the d- and p- character of Ni and O, respectively. You should see that there are 8 bands around the Fermi energy that capture most of these orbital characters. These bands are bands 2-9, between roughly 2 and -8 eV.

From these states, the target space is selected. You can see that these 8 Bloch bands (-8 to 2 eV) have mixed s-, p-, and d- character (cf. the PROCAR, which contains the spd- and site-projected wave function character of each orbital). We could select only the 5 d-states localized on the Ni atom. However, this leads to less localized basis functions overall. Instead, we include the O p-states (delocalized onto the O atom), which allows mixing and improves the description of the Wannier orbitals around the Fermi energy. We could also include the s-states; however, these are entangled with the next p-states (see the crossing at ~9 eV) so it would be difficult to select the s-states without also including these p-states. Even if these s-states were to be included, their contribution to the 5 bands around the Fermi energy (0 eV) is so small as to be negligible (toggle off the s- and p-states to check this).

Below you can see the contribution of the p- and d-states to the DOS using py4vasp. Try changing NiO_bands.dos.plot("d") or NiO_bands.dos.plot("p") to NiO_bands.dos.plot("Ni") or NiO_bands.dos.plot("O") to see where the states are localized.

Looking at the DOS, you can see that there is some contribution from the p-states to the d-dominated states around the Fermi energy (0 eV), and some from the d-states to the p-dominated states around -5 eV, showing the importance of including the p states in our basis. Additionally, you can see that at around 7 eV, there is almost no contribution from the d-states, and almost no contribution from the s-states around 0 eV, justifying our exclusion of the s-states. From these 8 Bloch bands, we can select the 8 Wannier orbitals. In the next exercise, we will do exactly this, generating the Wannier basis for the cRPA calculation.

1.5 Questions¶

1. What sort of smearing must you use for zero-bandgap or entangled systems? (hint: those that will use LFINITE_TEMPERATURE later)
2. What spd-character do the bands around the Fermi energy (0 eV) consist of?

2 Wannier projection

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• generate Wannier functions using Wannier90.
• plot their projection using gnuplot.

2.1 Task¶

Generate Wannier functions using Wannier90 for primitive NiO using a 12x12x12 k-mesh.

Wannier functions are an important tool for creating a small local basis [Phys. Rev. B 56 (1997) 12847, Phys. Rev. B 75 (2007) 195121]. The Bloch orbital $\psi_{n\bf k}^\sigma ({\bf r})$ is related to the cell-periodic function $u^\sigma_{n\bf k}({r})$ by:

$$\tag{1} \psi_{n\bf k}^\sigma ({\bf r}) = e^{i{\bf k r}} u^\sigma_{n\bf k}({r}) $$

where $n$ is the band index, $\textbf{k}$ is the crystal momentum, and $\sigma$ is the spin. A Wannier function $|w_{i\bf R}^\sigma \rangle$ is obtained from a linear combination of these Bloch orbitals via

$$\tag{2} |w_{i\bf R}^\sigma \rangle = \frac{1}{N_k}\sum_{n\bf k} e^{-i {\bf k R}} T_{i n}^{\sigma({\bf k})} | \psi_{n\bf k}^\sigma \rangle $$

where $N_k$ is the number of $\textbf{k}$-points, and $T_{i n}^{\sigma({\bf k})}$ is the transformation matrix. We use T here due to the prevalence of U in DFT+U calculations (cf. the widely used U-notation: $U_{mn\textbf{k}}$).

Usually, the basis set is localized such that the interaction between periodic images can be neglected. This allows us to work with the Wannier functions in the unit cell at $\bf R=0$:

$$\tag{3} |w_{i}^\sigma \rangle = \frac{1}{N_k}\sum_{n\bf k} T_{i n}^{\sigma({\bf k})} | \psi_{n\bf k}^\sigma \rangle $$

For a model Hamiltonian, e.g. cRPA, only a small subset of Bloch functions are used. These target states are typically centered around the chemical potential (or Fermi energy) and are strongly localized around ions. The model Hamiltonian can be solved successfully only if the target states are well represented by the Wannier basis.

In this exercise, you will project Wannier orbitals using Wannier90 and plot them to visualize their components.

2.2 Input files¶

The input files to run this example are prepared at $TUTORIALS/strongly_correlated/e02_wannier_projection.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

POSCAR
NiO

NiO
 4.171
  0.5    0.5    0.0
  0.0    0.5    0.5
  0.5    0.0    0.5
Ni O
1   1
Direct
 0.00  0.00  0.00
 0.50  0.50  0.50

INCAR:
LWRITE_WANPROJ determines whether the Wannier projection file WANPROJ is written; it contains the Wannier projection matrices $T$ for each $k$-point and band involved in the projection. NUM_WANN controls the number of Wannier functions to be constructed, WANNIER90_WIN sets the content of the wannier90.win file and LWANNIER90_RUN tells VASP to run the wannier90 library. Pay attention to include only 8 bands in the Wannier projection with the exclude_bands instruction. Other Wannier90 options define the number of Wannier functions, the number of bands used in the projection, and the orbital types on which the bands are projected. Finally, we also want to plot the corresponding fatbands using a gnuplot file. ALGO=NONE is set so that the wavefunction cannot be changed; this is a post-processing step and so we do not want to change the wavefunction at all.

ISMEAR = -1 ; SIGMA = 0.1

EDIFF = 1E-8
NBANDS = 96
ALGO = NONE
LWAVE = F
#
# write wanproj dataset
# allows to skip wannier projection in subsequent steps
LWRITE_WANPROJ = .TRUE.

# how many wannier states are there in total?
NUM_WANN = 8
# execute Wannier90 in library mode
LWANNIER90_RUN= T

WANNIER90_WIN = "
num_bands = 8
num_wann  = 8
# include only bands 2-9 in the projection
exclude_bands = 1, 10 - 96
# project d-like states on Ni and p-like ones on O
begin projections
Ni:d
O:p
end projections

# also interpolate band structure to validate
# Wannier projection
bands_plot = true
begin kpoint_path
 L  0.500  0.500  0.500  G  0.000  0.000  0.000
 G  0.000  0.000  0.000  X  0.500  0.000  0.500
 X  0.500  0.000  0.500  W  0.500  0.250  0.750
 W  0.500  0.250  0.750  L  0.500  0.500  0.500
 L  0.500  0.500  0.500  K  0.375  0.375  0.750
 K  0.375  0.375  0.750  G  0.000  0.000  0.000
end kpoint_path
# create data for fatbands too that
# visualize d-character on each band
bands_plot_project : 1 - 5
"

KPOINTS:
An equally spaced (12$\times$12$\times$12) Γ-centered grid of $\bf{k}$ points.

Automatically generated mesh
  0
Gamma
 12 12 12

POTCAR:

   PAW_PBE Ni_GW 31Mar2010
   PAW_PBE O_GW_new 19Mar2012           

You can see the contributions from the d and p states using gnuplot below. Create the figure by entering the following into the terminal:

gnuplot -e "set terminal pngcairo; set output 'bands.png'; load 'w90_2d_band_plot.gnu'"

Then take a look at it:

If the Wannier states have been well chosen, then the band structure using the Wannier basis, should look like that calculated earlier using the Bloch basis (cf. Exercise 1.4.2).

The band structure above provides important context for the cRPA calculation. The color coding reflects the orbital character of each band, with the Ni $d$-states shown in red and the O $p$-states in blue. The 5 Wannier bands around the Fermi energy (~0 eV) have largely $d$-character, with an important contribution from the $p$-states, while the three lower energy Wannier bands (around -7 eV) have largely $p$-character, with a significant contribution from the $d$-states.

The 5 $d$-states visible near the Fermi energy are the target states for the cRPA calculation. They will be used to describe the model Hamiltonian. You can see this by looking in the Wannier90 output file (wannier90.wout):

The first five Wannier functions are centered at the Ni atom (0.0, 0.0, 0.0), and the next three at the O atom (2.0855, 2.0855, 2.0855). These correspond to the 5 $d$-orbitals and the 3 $p$-orbitals. You want to select the orbitals with the smallest spread, i.e., the $d$-orbitals centered on the Ni atom ~0.45, rather than the O atom ~0.87 (cf. Eq. 14 of [Phys. Rev. B 56, 12847 (1997)] for more details of the spread).

2.4 Repeat for $d$-only (optional)¶

If you are interested in seeing the effect that using only the d-orbitals would have, go into the d-only directory and run the calculations there.

cd $TUTORIALS/strongly_correlated/e02_wannier_projection/d-only
cp ../e01_*/unoccupied/WAVECAR .
mpirun -np 4 vasp_std

cd $TUTORIALS/strongly_correlated/e02_wannier_projection/d-only
vimdiff ../INCAR INCAR

grep -A 5 "Final State" e02_*/wannier90.wout | tail -5

  WF centre and spread    1  ( -0.000000, -0.000000, -0.000000 )     1.69142625
  WF centre and spread    2  ( -0.000000, -0.000000, -0.000000 )     0.95788648
  WF centre and spread    3  ( -0.000000, -0.000000, -0.000000 )     0.95788572
  WF centre and spread    4  ( -0.000000, -0.000000, -0.000000 )     1.69142209
  WF centre and spread    5  ( -0.000000, -0.000000, -0.000000 )     0.95788435

2.5 Questions¶

1. What does the WANPROJ file contain?
2. What orbital characters do the target states have?
3. Which ions are the target states centered around?

3 Constrained random-phase approximation

$\uparrow$

¶

By the end of this tutorial, you will be able to:

• Calculate the cRPA Hubbard U, U' and J parameters using a coarse and fine k-mesh.
• Calculate the cRPA off-center terms.
• Refit the Ohno potential to the Hubbard U potential and investigate the frequency dependence of plasmons.

3.1 Task¶

Calculate the Hubbard U, U', and J parameters for NiO (primitive) using cRPA with 8 Wannier orbitals, comparing the effects of k-mesh (4x4x4 and 6x6x6) and the number of frequency points (12 and 24) on the Coulomb potential.

With carefully chosen Wannier orbitals $w_{a}^{*\sigma}({\bf r})$, the constrained random-phase approximation (cRPA) interaction can be calculated. This will provide the effective interaction matrix $U_{ijkl}$ in the Wannier basis stored in the WANPROJ file:

$$\tag{4} U^{\sigma\sigma'}_{ijkl} = \int {\rm d}{\bf r}\int {\rm d}{\bf r}' w_{i}^{*\sigma}({\bf r}) w_{j}^{\sigma}({\bf r}) U({\bf r},{\bf r}',\omega) w_{k}^{*\sigma'}({\bf r}') w_{l}^{\sigma'}({\bf r}') $$

where $U({\bf r},{\bf r}',\omega)$ is the Coulomb potential in frequency space $\omega$.

In practice, one often looks at the Hubbard-Kanamori interactions as a guideline $$\tag{5} {\cal U }^{\sigma\sigma'} = \frac 1 N \sum_{i\in \cal T} U_{iiii}^{\sigma\sigma'} $$ $$\tag{6} {\cal U' }^{\sigma\sigma'} = \frac{1}{N(N-1)}\sum_{i,j \in{\cal T}, i\neq j}^N U_{ijji}^{\sigma\sigma'} $$ $$\tag{7} {\cal J }^{\sigma\sigma'} = \frac{1}{N(N-1)} \sum_{i,j \in{\cal T}, i\neq j}^N U_{ijij}^{\sigma\sigma'} $$

These are for a single unit cell (at $R=0$; cf. Vijkl and Uijkl). For the extended systems ($R>0$), the corresponding Coulomb integrals (cf. VRijkl and URijkl) are discussed later.

For the actual DFT+U calculations [Phys. Rev. B 44 (1991) 943,Phys. Rev. B 57 (1998) 1505], a spherical averaging is recommended, which will be discussed in the DFT+DMFT tutorial.

3.2 Input files¶

The input files to run this example are prepared at $TUTORIALS/strongly_correlated/e03_CRPA.

VASP looks in the current directory for four main input files, i.e., POSCAR, INCAR, KPOINTS, and POTCAR. Check them out!

POSCAR
NiO

NiO
 4.171
  0.5    0.5    0.0
  0.0    0.5    0.5
  0.5    0.0    0.5
Ni O
1   1
Direct
 0.00  0.00  0.00
 0.50  0.50  0.50

./INCAR and fine_grid/INCAR:
To perform a cRPA calculation, specify ALGO=CRPA ("constrained RPA"). This step requires a localized basis that can be set with LOCALIZED_BASIS. If a file WANPROJ is present in the working directory, VASP will use it as the corresponding basis. In this tutorial, we reuse the Wannier projection from previous step. NTARGET_STATES controls which Wannier states are excluded from the cRPA. LFINITE_TEMPERATURE switches on the finite-temperature formalism (requires Fermi smearing, ISMEAR=-1) for cRPA and is necessary to perform calculations with zero-bandgap or entangled systems. Do not forget to set NBANDS to the same value that was used to generate the WAVECAR file in the previous step.

SYSTEM = NiO
ISMEAR = -1 ; SIGMA = 0.1
NBANDS = 96
# this options selects the wannier states in the target space
# the screening of these states will be removed from
# the effective interaction in GW
NTARGET_STATES = 1 2 3 4 5

ALGO = CRPA
LSCRPA = T
# adding following option to also compute multi-center terms
LTWO_CENTER = T
# use this option to sample a specific point on the
# real frequency axis
# OMEGAMAX = 3.2

# finite temperature can be switched off
# when calculating only one single frequency point
# this makes the code run faster
#LFINITE_TEMPERATURE = T

plasmons/INCAR:

SYSTEM = NiO
ISMEAR = -1 ; SIGMA = 0.1

EDIFF = 1E-8
NBANDS = 96
ALGO = CRPAR
# use more frequency points
NOMEGA = 12
# always use Matsubara formalism in euclidean spacetime
LFINITE_TEMPERATURE = T
NTARGET_STATES = 1 2 3 4 5
LSCRPA = T

# restrict memory
MAXMEM = 64000

./ and ./plasmons/KPOINTS:
An equally spaced (4$\times$4$\times$4) Γ-centered grid of $\bf{k}$ points.

Automatically generated mesh
  0
Gamma
 4 4 4

fine_grid/KPOINTS:

Automatically generated mesh
  0
Gamma
 6 6 6

3.3 cRPA: coarse k-mesh¶

To calculate the unoccupied states, go to the e03_CRPA directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e03_*
cp ../e01_*/unoccupied/WAVE{CAR,DER} .
cp ../e02_*/WANPROJ .
mpirun -np 4 vasp_std

Run vasp and look for lines of Hubbard U in the OUTCAR file.

Note, the complete data is stored in vaspout.h5.

3.4 cRPA: fine k-mesh¶

Now we can try increasing the $\textbf{k}$-point mesh and see what impact it has. Go to the fine_grid directory and take a look at the KPOINTS file. It is a 6x6x6 grid, rather than a 4x4x4 one (NB. the grid must be commensurate with the one used for the DFT step, i.e., 12x12x12; 8x8x8 would not be).

cd $TUTORIALS/strongly_correlated/e03_*/fine_grid
cp ../../e01_*/unoccupied/WAVE{CAR,DER} .
cp ../../e02_*/WANPROJ .
mpirun -np 4 vasp_std

Once the calculation has finished, take a look at the Coulomb interaction as before:

By increasing the $k$-mesh, the bare interaction barely changes, showing that it is already converged. We want the screened Hubbard U ($U$), u ($U'$), and J ($J$) terms for DMFT, so it is these that we should focus on. J is already converged, barely changing with increasing $k$-mesh. U and u are not yet converged (with 0.1 eV), so would require a still denser $k$-mesh.

3.5 cRPA: off-center terms¶

LTWO_CENTER was set in the cRPA calculation. This calculates the off-center Coulomb integrals. Thus the dataset uijkl in vaspout.h5 contains an additional dimension for the Wannier center R.

This dataset stores the following integrals (where the $R$ tells us how the Coulomb potential decays between cells): $$\tag{8} U_{ijkl}(R) = \int {\rm d}r \int {\rm d}r' w_i({\bf r}-{\bf R}) w_j({\bf r}-{\bf R}) U({\bf r}, {\bf r}',\omega=0) w_k({\bf r}') w_l({\bf r}') $$ The spatial decay of the Hubbard Kanamori $U$ is sometimes studied with the Ohno potential $$\tag{9} \frac{U(R)}{U(0)} = \sqrt{\frac{\delta}{R+\delta}} $$ Take a look at the data structure in the vaspout.h5 file:

This model can be plotted with py4vasp, which adjusts $\delta$ to the selected data. Try changing the radius_max parameter to restrict the points considered in fitting.

Each point is the nearest neighbor interaction between Ni sites in the crystal. O has been excluded by only taking the 5 $d$-orbitals as our basis for cRPA.

3.5.1 Filter noisy data and refit Ohno potential¶

The Ohno potential follows roughly the Hubbard U potential but is significantly distorted by noisy data points for longer distances [Theor. Chim. Acta 2 (1964) 219, npj Computational Materials 10 (2024) 286]. This noise comes mostly from insufficient $k$-point sampling and can be reduced by improving $k$-point sampling. Alternatively, the agreement can be improved by removing some of the noisy data points from the fit. We will do this in a few steps:

1. Extract the data from the above plot.
2. Remove outliers (like those above 0.74 eV beyond 10A).
3. Plot the filtered data and fit the Ohno potential.

NB. In this tutorial, we have removed the points that are due to noise in order to show how the fitting should work. In practice, you should increase the $k$-points mesh and the planewave energy cutoff, as the noise comes from using too sparse an FFT grid.

3.6 Plasmons: Frequency dependence¶

The cRPA calculations can be evaluated for large unit cells using a Euclidean spacetime algorithm [Comput. Phys. Commun., 117 (1999) 21100174-X), Phys. Rev. B 94 (2016) 165109]. This approach performs all computations internally on the imaginary time axis. In this regard, the ordinary Minkowski spacetime metric $\eta=(-+++)$ is changed to $\eta=(++++)$ and thus becomes one of Euclidean space $\mathbb{R}^4$, allowing for sparse representations of the effective Coulomb kernel.

This approach allows you to investigate the effective interaction parameters for large systems using more memory than a conventional cRPA calculation in Minkowski spacetime.

This computation is memory-intensive and takes longer for small systems compared to the Minkowski algorithm, but is much faster for larger systems. Using four MPI ranks and a 4x4x4 $\textbf{k}$-point mesh requires about 8GB of RAM.

NB. There is a long pause after estimated memory requirement per rank ****.* MB, per node ****.* MB while nothing is printed; the calculation is still running.

cd $TUTORIALS/strongly_correlated/e03_CRPA/plasmons/
cp ../WAVE{CAR,DER} . 
cp ../WANPROJ . 
mpirun -np 4 vasp_std

Inspect the vaspout.h5 file to check what results are now stored there:

The calculation determines the Hubbard interaction at NOMEGA imaginary frequency points and stores it in the uijkl dataset. This data is of minor practical interest. The actual result we are interested in is the Hubbard interaction as a function of real-frequency. The real-frequency dependence can be obtained from analytic continuation by fitting the imaginary frequency data to a Padé expression of the form:

$$\tag{10} U(z) = \frac{ (z-u_1)(z-u_2)\cdots(z-u_N) }{ (z-v_1)(z-v_2)\cdots(z-v_M) }, \quad u_1, v_1, \cdots, u_N, v_M \in \mathbb{C} $$

In py4vasp, the AAA algorithm is used to determine the roots $u_1,v_1,\cdots,u_N,v_M$ of the polynomial in the numerator and denominator, respectively [SIAM J. Sci. Comput. 20 (2018) A1494]. However, it can yield almost identical roots for the numerator and denominator, which should cancel out in an exact calculation. This phenomenon is known as Froissart doublets. Due to numerical noise, they often remain in the expression and deteriorate the analytical continuation. The clean_up_tol parameter can be used as a clean up of these doublets.

You can see several peaks and a wavy structure. The first peak (~ 5.5 eV) is a plasmon peak and corresponds roughly to the distance between the center of mass of the $p$- and $d$- states around the Fermi energy (cf. the DOS and Wannier functions). This is the most common transition in the system at low energy. Changing the clean_up_tol parameter does not change the position of the peak significantly, indicating that it is a real plasmon peak that does not stem from a Froissart doublet. Before this peak (<5 eV), there is almost constant behavior of the effective Coulomb interaction, indicating that the static approximation to the Coulomb potential in Hubbard-like models for NiO is appropriate.

At high frequencies, the screened Coulomb potential U begins to approach the bare Coulomb potential V. To resolve higher frequencies, more points are required.

Let's try increasing NOMEGA to see how this improves the resolution at higher frequencies by increasing the number of integration points from $12$ to $24$. For the sake of conserving time, we will not run this calculation now. If you have the time, you can run it yourself.

cd $TUTORIALS/strongly_correlated/e03_CRPA/plasmons/nomega_24
cp ../../../e01_*/unoccupied/WAVE{CAR,DER} .
cp ../../../e02_*/WANPROJ .
cp ../{INCAR,POSCAR,POTCAR,KPOINTS} .
sed -i 's/NOMEGA = 12/NOMEGA = 24/g' INCAR
mpirun -np 4 vasp_std

With $24$ imaginary frequency points, the AAA fit shows more structure at higher frequencies. The wavy structure beyond $50$ eV is much more sensitive to the clean_up_tol parameter at lower frequencies, indicating the presence of Froissart doublets. This numerical noise is usually reduced by increasing the number of data points. Since, the Minimax isometry method selects the imaginary frequency points already very well [Phys. Rev. B 101 (2020) 205145], increasing NOMEGA will not improve the result here. In contrast, increasing the number of $k$-points stabilizes the high frequency (Note the $24$ imaginary frequency points, rather than $12$) structure and a lower clean_up_tol setting can be chosen. This computation takes relatively long with 4 MPI ranks and uses about 12 GBs of memory.

You may skip this computation and focus on the analysis, since the result is already present in the corresponding path.

This makes the two plasmon peaks more distinct and the Coulomb potential more stable at higher frequencies. The peaks correspond to two different transitions. The first peak at 5.5 eV, corresponds roughly to the distance between the $d$ and $p$ centers, while the main plasmon peak at 50 eV corresponds to collective excitations, typical for GW-type approximations, such as cRPA. At higher frequencies, we can see how the screened Coulomb potential $U$ begins to approach the bare Coulomb potential $V$. It is only at higher density $k$-point meshes and more imaginary frequency points that the intermediate structure is better resolved.

The key takeaway from this is that the analytic continuation is not a black box. You can see that the first plasmon peak does not change much as we change the $k$-mesh and the number of frequency points. This tells us that this peak is reliable. This low-frequency regime is where the most interesting physics happens for model Hamiltonians. The peak at higher frequencies (~50 eV) is much harder to converge. This peak corresponds to roughly the energy of Gamma-radiation, and so is not likely to have much meaning for us. These higher frequencies are much harder to converge and so the results not as reliable. When you do a calculation, do not worry too much about these; focus on converging the low-frequency peaks that have a clear physical meaning.

3.7 Questions¶

1. What does the first plasmon peak at 5.5 eV correspond to?
2. What does the main plasmon peak at 50 eV correspond to?
3. What happens to the screened Coulomb potential at high frequencies?
4. What tag do you need to include for zero-bandgap or entangled systems in the cRPA?

Well done! You have finished Part 1 of the strongly correlated tutorial!¶

Keep going and get started with Part 2. In the file browser, navigate to $TUTORIALS/strongly_correlated and open part2.ipynb!

Go to Top $\uparrow$