DFT and projections

We will perform DFT+DMFT calculations for the charge-transfer insulator NiO. We start from scratch and provide all necessary input files to do the calculations: First for doing a single-shot calculation (and then for charge-self consistency).

VASP setup

We start by running a simple VASP calculation to converge the charge density initially. Use the INCAR, POSCAR, and KPOINTS provided and use your own POTCAR file.

Let us take a look in the INCAR, where we have to specify local orbitals as basis for our many-body calculation.

System  = NiO

# for PLOT DOS
# ISMEAR = -5
# better convergence for small kpt grids
ISMEAR = 2
SIGMA = 0.05

# converge wave functions
EDIFF = 1.E-7
NELMIN = 35

# the energy window to optimize projector channels (absolute)
EMIN = -3
EMAX = 10

LMAXMIX = 6

# switch off all symmetries
ISYM = -1

# project to Ni d and O p states
LORBIT = 14
LOCPROJ = "1 : d : Pr
           2 : p : Pr"

LORBIT = 14 takes care of optimizing the projectors in the energy window defined by EMIN and EMAX. We switch off all symmetries with ISYM=-1 since symmetries are not implemented in the later DMFT scripts. Finally, we select the relevant orbitals for atom 1 (Ni, d-orbitals) and 2 (O, p-orbitals) by the two LOCPROJ lines. For details refer to the VASP wiki on the LOCPROJ flag. The projectors are stored in the file LOCPROJ.

PLOVASP

Next, we post-process the projectors, which VASP stored in the file LOCPROJ. We do this by invoking plovasp plo.cfg which is configured by an input file, e.g., named plo.cfg.

[General]
# DOSMESH = -21 55 400
HK = True 
COMPLEMENT = True

[Group 1]
SHELLS = 1 2
EWINDOW = -9 2
NORMION = False
NORMALIZE = True
BANDS = 2 9

[Shell 1] # Ni d shell
LSHELL = 2
IONS = 1
CORR = True
TRANSFORM = 1.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0
            0.0  0.0  0.0  1.0  0.0
            0.0  0.0  0.0  0.0  1.0

[Shell 2] # O p shell
LSHELL = 1
IONS = 2
CORR = False
TRANSFORM = 1.0  0.0  0.0
            0.0  1.0  0.0
            0.0  0.0  1.0

Here, in [General] we set the basename and the grid for calculating the density of states. In [Group 1] we define a group of two shells which are orthonormalized with respect to states in an energy window from -9 to 2 for all ions simultaneously (NORMION = False). We define the two shells, which correspond to the Ni d states and the O p states. Only the Ni shell is treated as correlated (CORR = True), i.e., is supplemented with a Coulomb interaction later in the DMFT calculation. Here, we chose to use the Hamiltonian mode of the vasp converter by setting COMPLEMENT=TRUE, and specifying to use explicitly only bands with indices 2 to 9 (BANDS). This is optional but later used in the post-processing.

Converting to hdf5 file

We run the whole conversion to a dft_tools readable h5 archive by running the converter script provided python converter.py

from triqs_dft_tools.converters.vasp import *
import triqs_dft_tools.converters.plovasp.converter as plo_converter

# Generate and store PLOs
plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')

# run the converter
Converter = VaspConverter(filename = 'vasp', proj_or_hk = 'hk')
Converter.convert_dft_input()

A h5 archive should be created with the name vasp.h5 Now we are all set to perform a dmft calculation.

DMFT

dmft script

Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on SrVO3, we will not go into detail here but simply provide the script nio.py. Following Kunes et al. in PRB 75 165115 (2007) we use \(U=8\) and \(J=1\). We select \(\beta=5\) instead of \(\beta=10\) to ease the problem slightly. For simplicity we fix the double-counting potential to \(\mu_{DC}=59\) eV by:

DC_value = 59.0
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_value=DC_value)

For sensible results run this script in parallel on at least 20 cores. As a quick check of the results, we can compare the orbital occupation from the paper cited above (\(n_{eg} = 0.54\) and \(n_{t2g}=1.0\)) and those from the cthyb output (check lines Orbital up_0 density: for a t2g and Orbital up_2 density: for an eg orbital). They should coincide well.

Local lattice Green’s function for all projected orbitals

We calculate the local lattice Green’s function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script NiO_local_lattice_GF.py. The result is saved in the h5 file as G_latt_orb_it<n_it>, where <n_it> is the number of the last DMFT iteration.

Spectral function on real axis: MaxEnt

To compare to results from literature we make use of the maxent triqs application and calculate the spectral function on real axis. Use this script to perform a crude but quick calculation: maxent.py using a linear real axis and a line-fit analyzer to determine the optimal \(\alpha\). The output is saved in the h5 file in DMFT_results/Iterations/G_latt_orb_w_o<n_o>_it<n_it>, where <n_o> is the number of the orbital and n_it is again the number of the last iteration. The real axis information is stored in DMFT_results/Iterations/w_it<n_it>.

Charge self-consistent DMFT

In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP INCAR such that VASP reads the updated charge density after each step. We add the lines:

ICHARG = 5
NELM = 1000
NELMIN = 1000
IMIX=1
BMIX=0.5
AMIX=0.02

which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file GAMMA. We change the mixing here to stabilize the updating, which can be problem for charge ordered systems. Vasp is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number.

We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the CHGCAR and nio.h5 together with the other VASP input files and plo.cfg in a new directory. We use a script called vasp_dmft to invoke VASP in the background and start the DMFT calculation together with plovasp and the converter. This script assumes that the dmft sript contains a function dmft_cycle() and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforwardly (for a working example see nio_csc.py). The most important new lines are:

SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
SK.calc_density_correction(dm_type='vasp')
correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()

where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file GAMMA. The correlation energy is calculated via Migdal-Galitski formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.

To help convergence, we keep the density (i.e., the GAMMA file) fixed for a few DFT iterations. We do so since VASP uses an iterative diagonalization. Within these steps we still need to update the projectors and recalculate the GAMMA file to not shuffle eigenvalues around by accident.

We can start the whole machinery by executing:

vasp_dmft -n <n_procs> -i <n_iters> -j <n_iters_dft> -p <vasp_exec>  nio_csc.py