Part 1: Spin-polarized calculations ¶
Content¶
1 Density of states and band structure for ferromagnetic hexagonal-close-packed Co
2 Magnetic unit cell of antiferromagnetic body-centered-cubic Cr
3 Density of states of $e_g$ and $t_{2g}$ orbitals in antiferromagnetic NiO
4 Including on-site Coulomb interaction and pseudopotential selection for NiO relaxation
By the end of this tutorial, you will be able to:
- perform a spin-polarized calculation
- obtain the total magnetization and local magnetic moments
- plot the projected density of states (DOS) for a spin-polarized system
- compute the band structure in a second step after self-consistency has been reached
1.1 Task¶
Perform a spin-polarized DFT calculations for hexagonal-close-packed (HCP) cobalt with initial ferromagnetic moments and obtain the local magnetic moment as well as the DOS. Then, perform a band-structure calculation along the high-symmetry path Γ-M-K-Γ.
In a spin-polarized calculation with ISPIN = 2, the Kohn–Sham (KS) orbitals are two component vectors in spin space. Here, the two components correspond to spin up and down with a global spin-quantization axis. The spin-up and spin-down KS orbitals are allowed to have different eigenvalues. This allows describing collinear magnetic structures also for itinerant electrons. The interaction between the two spin components stems solely from the exchange-correlation functional and potential that depends on the spin-up and spin-down charge densities.
Mind that without spin-orbit coupling the direction of the spin in real space is arbitrary. Only the direction relative to neighboring magnetic moments is relevant and thus ferromagnetic and antiferromagnetic order can be distinguished. In practice, you can perform spin-polarized calculations with the executable vasp_std. Only when spin-orbit coupling is included or noncollinear magnetic structures (with a local spin-quantization axis) are considered, do you need to use the executable vasp_ncl.
1.2 Input¶
The input files to run this example are prepared at $TUTORIALS/magnetism/e01_hcp-Co. Check them out!
Co2
1.0
2.4719999032000000 0.0000000000000000 0.0000000000000000
-1.2359999516013529 2.1408147143262166 0.0000000000000000
0.0000000000000000 0.0000000000000000 4.0211440800000000
Co
2
direct
0.3333333333333333 0.6666666666666667 0.2500000000000000 Co
0.6666666666666667 0.3333333333333333 0.7500000000000000 Co
POSCAR.equiv
Co2
2.4719999032
1.0000000000000000 0.0000000000000000 0.0000000000000000
-0.5000000000000000 0.8660254037853866 0.0000000000000000
0.0000000000000000 0.0000000000000000 1.6266764714653246
Co
2
direct
0.3333333333333333 0.6666666666666667 0.2500000000000000 Co
0.6666666666666667 0.3333333333333333 0.7500000000000000 Co
INCAR.scf
SYSTEM = FM Co
ALGO = Normal
PREC = Accurate ! accurate FFTs
EDIFF = 1e-5 ! convergence criterion
ENCUT = 350 ! energy cutoff
ISPIN = 2
MAGMOM = 2 2 ! initial magnetic moments
LMAXMIX = 4
LASPH = T
ISMEAR = 2 ! Methfessel-Paxton method
SIGMA = 0.1 ! broadening
INCAR.dos
SYSTEM = FM Co
ALGO = Normal
PREC = Accurate ! accurate FFTs
EDIFF = 1e-5 ! convergence criterion
ENCUT = 350 ! energy cutoff
ISPIN = 2
LMAXMIX = 4
LASPH = T
# DOS
ICHARG = 11
ISMEAR = -5
NEDOS = 5001
# on-site spin moments and partial DOS
LORBIT = 11
band/INCAR
SYSTEM = FM Co
ALGO = Normal
PREC = Accurate ! accurate FFTs
EDIFF = 1e-5 ! convergence criterion
ENCUT = 350 ! energy cutoff
ISPIN = 2
LMAXMIX = 4
LASPH = T
# bands
ICHARG = 11
LORBIT = 11
Regular k-point mesh
0 ! 0 -> determine number of k points automatically
Gamma ! generate a Gamma centered mesh
6 6 4 ! subdivisions N_1, N_2 and N_3 along the reciprocal lattice vectors
0 0 0 ! optional shift of the mesh (s_1, s_2, s_3)
band/KPOINTS
k points for band structure in HCP structure
10
line-mode
reciprocal
0.0000000000 0.0000000000 0.0000000000 Γ
0.5000000000 0.0000000000 0.0000000000 M
0.5000000000 0.0000000000 0.0000000000 M
0.3333333333 0.3333333333 0.0000000000 K
0.3333333333 0.3333333333 0.0000000000 K
0.0000000000 0.0000000000 0.0000000000 Γ
Pseudopotential of Co.
The POSCAR file contains two Cobalt atoms in a hexagonal structure. Can you recall what the keyword direct in line 8 means?
Click to see the answer!
The positions are provided in direct coordinates, also called fractional coordinates. The alternative is specifying the positions in Cartesian coordinates.
There are three INCAR files provided: First, INCAR.scf to obtain a self-consistent solution for the spin-up and spin-down charge densities. Here,
- ISPIN turns on spin-polarized calculations.
- The magnetic moments are initialized on both Co sites to 2$\mu_B$ with MAGMOM.
- As we wish to restart the calculation, it is important to set LMAXMIX to 4 for $d$ electrons (and 6 for $f$ electrons). Otherwise the charge density written to the CHGCAR file is incomplete. But in any case, it might speed up convergence as the charge-density mixer consideres states up to LMAXMIX.
- LASPH should be considered a standard setting, especially for magnetic calculations. It includes aspherical contribution to the gradient corrections inside the PAW spheres.
- ISMEAR= 2 selects the Methfessel-Paxton method of order 2, which is good for metallic systems, but should be avoided for semi-conductors and insulators.
Second, INCAR.dos includes settings to obtain a smooth DOS:
- ICHARG = 11 reads the charge density from the CHGCAR file and fixes it throughout the calculation. Note that for a spin-polarized calculation with ISPIN=2, the CHGCAR file contains the spin-up and spin-down charge densities.
- ISMEAR = -5 sets the tetrahedron method with Blöchl corrections, which yields smoother DOS but is less suited for the electronic convergence.
- NEDOS sets the number of grid points on which the DOS is evaluated, and
- LORBIT = 11 must be set to obtain the projected DOS, i.e., in terms of angular-momentum and magnetic quantum numbers $l$ and $m$. Additionally, the partial magnetization can be found in the OUTCAR file.
Third, INCAR.band is provided. Check if you understand all tags and search the VASP Wiki if need be.
The KPOINTS file sets a uniform $\mathbf{k}$ mesh in the first Brilloun zone for the self-consistent field (SCF) and the DOS calculation. Why do we choose the same number of subdivisions for the first two directions and fewer in the third direction? Consider the structure defined in the POSCAR file to find the answer.
Click to see the answer!
The first two lattice vectors, $\mathbf{a}_1$ and $\mathbf{a}_2$, have the same length: $|\mathbf{a}_1|=|\mathbf{a}_2|$. Note that $\mathbf{a}_1$ points in x direction, while $\mathbf{a}_2$ points in $(-1/2, \sqrt(3)/2,0)$. The third lattice vector, $\mathbf{a}_3$, points in z direction and is longer by roughly a factor of $|\mathbf{a}_3|/|\mathbf{a}_1|\approx 1.6$. To emphasize these properties already in the input, you could also rewrite the POSCAR file as in POSCAR.equiv.
With this in mind, recall that a long lattice vector in real space, leads to a short lattice vector in reciprocal space. So, what ever devision you choose for the third direction, you should choose roughly a factor of $1.6$ more subdevisions to obtain the same resolution along the other reciprocal lattice vectors.
The band/KPOINTS file contains a high-symmetry path for the band-structure calculation.
1.3 Calculation¶
Open a terminal and navigate to this example's directory, copy the INCAR file for the self-consistent (scf) calculation and run it by entering the following:
cd $TUTORIALS/magnetism/e01_*
cp INCAR.scf INCAR
mpirun -np 2 vasp_std
In the standard output (to the terminal) you find a summary stating the total energy and the total magnetization of the system. How does the total magnetization (mag=) compare to the experimental value of $1.7 \mu_B$/Co [T. WAKIYAMA, J. Phys. Colloques 32 (1971) C1-340-C1-341]? Plot the structure to get a hint!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_hcp_Co")
mycalc.structure.plot()
Click to see the answer!
There are two Co atoms in the system. So, to compare to the experimental value, the total spin magnetization should be divided by two. The computed value is thus a bit smaller but in relatively good agreement.
While self-consistency is reached much more reliably using the Methfessel-Paxton method, it is better to use the tetrahedron method to obtain a smooth DOS. The partial DOS is only computed when LORBIT>=10. In a second step, we thus fix the charge density and compute the partial DOS with the tetrahedron method:
cp INCAR.dos INCAR
mpirun -np 2 vasp_std
Plot the partial DOS using py4vasp! What are the main features?
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_hcp_Co")
mycalc.dos.plot( "up, down" )
mycalc.dos.plot( "d" )
Click to see the answer!
Most of the magnetization emerges from the d orbitals. The spin-up and spin-down DOS have similar shape but are shifted in energy, such that there is a majority and a minority spin. The relative energy shift is the so-called exchange splitting in the context of the Stoner model.
To perform further analysis, you can get the data as a dictionary using py4vasp:
mycalc.dos.to_dict("d")
Alternatively, check out the PROCAR file.
mycalc.local_moment.to_dict()
sum(mycalc.local_moment.to_dict()['total'][0])
Click to see the answer!
These are the local magnetic moments in the $s$, $p$ and $d$ orbital on each site. To obtain the total magnetic moment per site the values in the python dictionary have to be summed. Mind that in this calculation only the magnetic moments due to the electronic spin, but not due to the orbital angular momentum, is considered. In $d$ orbitals the orbital moments are small compared to the spin moments, but for $f$ orbitals the contribution might be sizable. It is possible to consider the orbital moments in VASP, but this is covered in part 2.
Do the on-site magnetic moments correspond to the so-called magnetic moments as defined in equation 1?
$$\tag{1} m_{\textrm{atom}} = g m \mu_B $$
Here, $g$ is the $g$ factor), $m$ is the magnetic quantum number and $\mu_B$ is the Bohr magneton. The magnetic moment in Eq. (1) is, e.g., measured in neutron scattering spectroscopy. Or does the computed value correspond to the effective magnetic moment?
$$\tag{2} m^{\textrm{eff}}_{\textrm{atom}} = g \sqrt{(s(s+1))} \mu_B $$
$s$ is the spin quantum number. The effective magnetic moment can be measured in muon scattering or indirectly obtained through Curie's law by measuring the specific heat.
Plot the local magnetic moments using py4vasp! But before, think about which direction the local moments have. Why?
Click to see the answer!
Without spin-orbit coupling the spin-quantization axis is not coupled to real space. Thus, the calculation contains no information about the direction of the magnetic moments. We only know their relative orientation to each other. As the partial magnetization has the same sign for both Co sites, we know the on-site moments are ferromagnetically coupled, as expected when we set the initial magnetic moments using the MAGMOM tag. The plot below hence shows a tentative global direction.
mycalc.local_moment.plot()
Next, run the band-structure calculation in the subdirectory band starting from the converged charge density:
cp CHGCAR POSCAR POTCAR band/.
cd band
mpirun -np 2 vasp_std
Now, plot the band structure using py4vasp! What effect does the exchange splitting we saw in the DOS have on the bands?
import py4vasp
my_band_calc = py4vasp.Calculation.from_path( "./e01_hcp_Co/band")
my_band_calc.band.plot()
Click to see the answer!
The bands of the minority spin are moved up in energy by the exchange splitting. See, for instance, the three spin-up bands around $-0.5\,$eV between $M$ and $K$. The same dispersion can be seen around $1.3\,$eV in the spin-down bands.
1.4 Questions¶
- What does the ISPIN tag set?
- Which INCAR tag must be set to obtain local magnetic moments?
- What is exchange splitting and where can it be observed?
By the end of this tutorial, you will be able to:
- perform a spin-polarized calculation for an antiferromagnetic (AFM) material
- restart a spin-polarized calculation
2.1 Task¶
Perform a spin-polarized calculation for body-centered-cubic (BCC) Cr with AFM order with propagation vector $\mathbf{q}=(1/2,1/2,1/2)$. Consider convergence with respect to k points and compute the band structure using a KPOINTS_OPT file. Finally, compute the partial DOS.
The propagation vector $\mathbf{q}$ describes how the magnetic moments vary spatially within a bulk material. It represents the spatial periodicity of the magnetic order, and can be used to define the magnetic unit cell, which might be larger than the crystallographic unit cell.
Let's consider BCC Cr as an example. In its paramagnetic phase, there is only one atom per primitive unit cell. However, in the antiferromagnetic (AFM) phase, the primitive unit cell contains two atoms, one located at fractional coordinates $(0,0,0)$ and the other at $(0.5,0.5,0.5)$. In this AFM configuration, the magnetic moments on the Cr atoms alternate in direction, and the magnetic unit cell is twice the size of the crystallographic unit cell.
To describe this AFM order, we can introduce a propagation vector $\mathbf{q}=(1/2,1/2,1/2)$ in units of the reciprocal lattice vectors. This vector indicates that the magnetic moments oscillate with a periodicity of two lattice constants, which corresponds to the doubling of the unit cell along the body diagonal.
2.2 Input¶
The input files to run this example are prepared at $TUTORIALS/magnetism/e02_bcc_Cr. Check them out!
BCC Cr
2.9688993551160281
1.000000000 0.000000000 0.000000000
0.000000000 1.000000000 0.000000000
0.000000000 0.000000000 1.000000000
Cr
2
Cartesian
0.000000000 0.000000000 0.000000000 Cr
0.500000000 0.500000000 0.500000000 Cr
INCAR.scf
SYSTEM = AFM Cr
ALGO = Normal
PREC = Accurate ! accurate FFTs
ENCUT = 350 ! energy cutoff
EDIFF = 1e-5 ! convergence criterion
ISPIN = 2
MAGMOM = 2 -2 ! initial magnetic moments
LMAXMIX = 4
LASPH = T
ISMEAR = 0 ! Gaussian smearing
SIGMA = 0.1 ! broadening
LORBIT = 11
LKPOINTS_OPT = F
INCAR.fine
SYSTEM = AFM Cr
ALGO = Normal
PREC = Accurate ! accurate FFTs
ENCUT = 350 ! energy cutoff
EDIFF = 1e-5 ! convergence criterion
ISPIN = 2
MAGMOM = -0.1 0.1 ! only to break the symmetry
LMAXMIX = 4
LASPH = T
ISMEAR = 0 ! Gaussian smearing
SIGMA = 0.1 ! broadening
LORBIT = 11
LKPOINTS_OPT = T
ICHARG = 1
INCAR.dos
SYSTEM = AFM Cr
ALGO = None ! only postprocessing
PREC = Accurate ! accurate FFTs
ENCUT = 350 ! energy cutoff
ISPIN = 2
MAGMOM = -0.1 0.1 ! only to break the symmetry
LMAXMIX = 4
LASPH = T
ISMEAR = -5 ! tetrahedron method
LORBIT = 11
LKPOINTS_OPT = F
ICHARG = 11
NEDOS = 5001
LH5 = T
KPOINTS.scf
Regular k-point mesh
0 ! 0 -> determine number of k points automatically
Gamma ! generate a Gamma centered mesh
4 4 4 ! subdivisions N_1, N_2 and N_3 along the reciprocal lattice vectors
0 0 0 ! optional shift of the mesh (s_1, s_2, s_3)
KPOINTS.fine
Regular k-points on finer mesh
0
Gamma
8 8 8
k points for band structure in BCC structure
10 ! intersections
line-mode
reciprocal
0.0000000000 0.0000000000 0.0000000000 Γ
0.0000000000 0.5000000000 0.0000000000 X
0.0000000000 0.5000000000 0.0000000000 X
0.5000000000 0.5000000000 0.0000000000 M
0.5000000000 0.5000000000 0.0000000000 M
0.0000000000 0.0000000000 0.0000000000 Γ
Pseudopotential of Cr.
The POSCAR file shows how the scaling factor in line 2 acts. In particular, not only the lattice matrix but also the Cartesian coordinates are scaled. The defined structure is the magnetic unit cell of BCC Cr. In this case, the conventional unit cell of the crystal structure coincides with the magnetic unit cell.
Concerning the INCAR.scf file, check that you understand the meaning of all tags. LKPOINTS_OPT = F ensures that the KPOINTS_OPT file is ignored for now. The initial magnetic moments are set to AFM order by MAGMOM = 2 -2. This breaks the symmetry of the crystal structure, i.e., the two Cr sites are no longer equivalent.
In the subsequent calculation, we will use INCAR.fine on a finer $\mathbf{k}$ mesh. Also, here, the MAGMOM is set in order to break the symmetry. Note that the magnetic moments are not initialized to the values of MAGMOM because ICHARG is set to 1 and, thus, the spin-up and spin-down charge densities are read from the CHGCAR file. We want to compute the band structure for the finer $\mathbf{k}$ mesh, so we set LKPOINTS_OPT=True, i.e., its default.
The INCAR.dos turns off electronic minimization with ALGO = None. This is convenient for postprocessing but requires a converged WAVECAR or vaspwave.h5 file.
The KPOINTS file and the KPOINTS.fine file define a Γ-centered $\mathbf{k}$ mesh with equally spaced $\mathbf{k}$ points. The only difference is the $\mathbf{k}$-mesh density.
The KPOINTS_OPT file is used in combination with the KPOINTS.fine file in the band-structure calculation.
2.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/magnetism/e02_*
cp INCAR.scf INCAR
cp KPOINTS.scf KPOINTS
mpirun -np 2 vasp_std
You can see in the summary line of the stdout that the total magnetization vanishes as expected for AFM order. Make sure to save the vaspout.h5 file for later:
cp vaspout.h5 bk.h5
Then print the on-site magnetic moments using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_file( "./e02_bcc_Cr/bk.h5")
mycalc.local_moment.print()
# or
[sum(mycalc.local_moment.to_dict()['total'][i]) for i in range(2)]
Check the computational cost by looking for Total CPU time used in the OUTCAR file. Create a backup for comparison and restart the calculation with a finer $\mathbf{k}$ mesh from the charge density obtained with the coarser $\mathbf{k}$ mesh!
cp OUTCAR OUTCAR.bk
cp INCAR.fine INCAR
cp KPOINTS.fine KPOINTS
rm WAVECAR
mpirun -np 2 vasp_std
Make sure to save the vaspout.h5 file for later:
cp vaspout.h5 bk2.h5
Note that we delete the WAVECAR file so there is no attempt to read it. The attempt would fail because the number of $\mathbf{k}$ points has changed, but the calculation is still influenced by the attempt. That is because when no WAVECAR file is present, the KS orbitals are initialized randomly and the first few steps of the electronic minimization the charge density is fixed. After a failed attempt to read a WAVECAR file, the charge density is not fixed for the first few steps and we do not take full advantage of having already converged on a coarser $\mathbf{k}$ mesh.
mycalc = py4vasp.Calculation.from_file( "./e02_bcc_Cr/bk2.h5")
mycalc.local_moment.print()
Compare the OUTCAR files! Search for LOOP to compare the computational effort. Then, go to the end of the files and compare the total charge and magnetization (x). What do you see?
Click to get a hint!
You may want to use one or all of these commands:
grep LOOP OUTCAR
grep LOOP OUTCAR.bk
vimdiff OUTCAR OUTCAR.bk
Click to see the answer!
The calculation including more $\mathbf{k}$ points is clearly more expensive. We note that the on-site magnetic moments were not converged on the coarser $\mathbf{k}$ mesh. The on-site charge, on the other hand, barely changed after increasing the $\mathbf{k}$-mesh density. This emphasizes that one needs to perform a convergence study exactly for the quantity of interest.
Plot the band structure along the $\mathbf{k}$ path specified in the KPOINTS_OPT file using py4vasp! You can click the label in the legend to activate and deactivate plotting a specific line.
mycalc.band.plot("kpoints_opt")
Finally, compute the DOS in a pure postprocessing calculation! Additionally, we set LH5=True in order to plot the magnetization.
cd $TUTORIALS/magnetism/e02_*
cp INCAR.dos INCAR
mpirun -np 2 vasp_std
Plot the DOS at one of the Cr sites! Is there any exchange splitting? How do the spin-up and spin-down DOS differ?
mycalc = py4vasp.Calculation.from_path( "./e02_bcc_Cr/")
mycalc.dos.plot("1")
Click to see the answer!
The majority spin has a larger DOS than the minority spin. There is no exchange splitting like in the ferromagnetic case.
Can you tell what the Cr DOS at site 2 will look like?
Click to see the answer!
Due to the AFM order, the majority and minority spin will be opposite to site 1.
mycalc.dos.plot("2")
Plot the magnetic structure using py4vasp!
mycalc.local_moment.plot()
While the local magnetic moments rely on some projection method (selected by LORBIT), the magnetization on the plane wave grid is a simple byproduct of spin-polarized DFT. It can be read from the vaspwave.h5 file and shown using a quiver plot as follows:
mycalc.density.to_contour("magnetization", a=0.5)+mycalc.density.to_quiver(a=0.5)
Here, the argument for the contour plot selects the magnetization and a plane orthogonal to the first lattice vector at the fraction of 0.5. For the quiver plot, the default density is the magnetization (since a charge density is a scalar in every point and cannot be nicely illustrated as a quiver plot) so only the plane is parsed as an argument. Hence, the graph shows the magnetization density of site 2.
2.4 Questions¶
- Does MAGMOM influence the symmetrization of the charge density?
- What is ALGO=None?
- Is it strictly necessary to converge the total energy with respect to the $\mathbf{k}$-mesh density?
- How can a calculation be restarted when the $\mathbf{k}$ mesh has changed?
By the end of this tutorial, you will be able to:
- set the Fermi energy mid-gap
- adjust the density mixer to converge more quickly
- plot combinations of orbitals in py4vasp
3.1 Task¶
Compute the partial DOS of Ni $e_g$ and $t_{2g}$ orbitals in AFM NiO.
To understand the magnetic properties of a system, it is often useful to consider the orbital splitting in crystal field theory as a first step. NiO is a transition metal oxide, where Ni$^{2+}$ is neighboring six equivalent O$^{2-}$ ions which together forms an NiO$_{6}$ octahedron. The Ni $d$ orbitals split into $t_{2g}$ and $e_g$ orbitals according to their symmetry properties with respect to the operations of the octahedral point group (O$_h$) as summarized in the character table of O$_h$. The octahedral point group is the group of symmetry operations for an octahedron and comprises all rotations on the Ni site that leave the system invariant.
The names $t_{2g}$ and $e_g$ are Mulliken symbols for irreducible representations. The $e_g$ set contains the d$_{z^2}$ and $d_{x^2-y^2}$ orbitals. The 'e' refers to doubly degenerate, i.e., there are two orbitals with the same energy. The 'g' refers to gerade, which is German for 'even' and means symmetric with respect to inversion. The $t_{2g}$ manifold contains the d$_{xy}$, d$_{xz}$, and d$_{yz}$ orbitals. The 't' refers to triply degenerate, and the 'g', as before, refers to gerade. '2' indicates that the manifold is anti-symmetric with respect to a vertical mirror plane perpendicular to the principal axis.
The lobes of the $e_{g}$ orbitals point towards the O$^{2-}$ ligands, while the $t_{2g}$ can maximally avoid the direction of the ligands. The negative point charge of the O$^{2-}$ ligands repels the $d$ electrons. Thus, $e_{g}$ is expected to be higher in energy due to the octahedral crystal field.
In the ionic limit, Ni$^{2+}$ has an electronic configuration of [Ar]3d$^8$. 6 of the electrons are expected to fully occupy the $t_{2g}$ orbitals, while 2 electrons may form a total-spin-$1$ configuration in the $e_{g}$ orbitals. So, we can expect a maximum value of $2\mu_B/$f.u. for the magnetization.
3.2 Input¶
The input files to run this example are prepared at $TUTORIALS/magnetism/e03_NiO. Check them out!
AFM NiO
4.17000000000000
1.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 1.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 1.0000000000000000
Ni O
2 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
SYSTEM = NiO
ISTART = 0
ISPIN = 2
MAGMOM = 2.0 -2.0 2*0
LASPH = T
EDIFF = 1E-8
ISMEAR = -5
EFERMI = MIDGAP
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
LORBIT = 11
EMIN = 3.5
EMAX = 10
NEDOS = 5001
k points
0
gamma
6 6 6
0 0 0
Pseudopotential of Ni and O_s.
Read about density mixing on the VASP Wiki! Consider the options of IMIX. For IMIX=4, the Pulay-mixing method is used. There are two coefficients for both the charge (AMIX, BMIX) and the magnetic degrees of freedom (AMIX_MAG, BMIX_MAG), respectively. For magnetic calculations, it can be advantages to use linear mixing, i.e. setting BMIX and BMIX_MAG to very small values, but not zero for numerical reasons.
The Fermi energy is set midgap using EFERMI = MIDGAP.
3.3 Calculation¶
Open a terminal, navigate to this example's directory, run VASP by entering the following:
cd $TUTORIALS/magnetism/e03_*
mpirun -np 2 vasp_std
Plot the DOS of the Ni $e_g$ and $t_{2g}$ orbitals using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e03_NiO")
dos = mycalc.dos.read("dxy, dyz, dz2, dxz, dx2y2")
(
py4vasp.plot(dos["energies"], dos["dxy_up"] + dos["dyz_up"] + dos["dxz_up"], "t2g") +
py4vasp.plot(dos["energies"], -dos["dz2_up"] -dos["dx2y2_up"], "eg")
)
The orbitals are combined to obtain the linear combination that yields the irreducible representations. The Fermi energy is in the middle of the gap.
Let us compute the size of the gap!
import numpy as np
def gap(calc, select=""):
homo = np.amax((calc.band.to_dict()['bands_up'][calc.band.to_dict()['occupations_up'] > 0.5]))
lumo = np.amin((calc.band.to_dict(select)['bands_up'][calc.band.to_dict(select)['occupations_up'] < 0.5]))
return lumo - homo
print("Band gap (in eV):", gap(mycalc))
The size of the gap is very clearly underestimated compared to the experimentally measured gap of about 4.2eV, see discussion by Dudarev et al., PRB 57, 1505 (1998). This is due to the inherent issue that DFT underestimates the exchange interaction. There are various strategies to improve the computed size of the gap. One is to include on-site Coulomb interaction $U$ which will be discussed in the next exercise.
Inspecting the DOS further, we see that the $t_{2g}$ orbitals are fully occupied. Now, consider the occupation of the $e_{g}$ orbitals on each Ni site! Click the labels in the legend to selectively (de)activate the lines.
dos = mycalc.dos.read("1(dz2, dx2y2) 2(dz2, dx2y2)")
(
py4vasp.plot(dos["energies"], dos["Ni_1_dz2_up"] + dos["Ni_1_dx2y2_up"], "Ni 1 eg up")+
py4vasp.plot(dos["energies"], -dos["Ni_1_dz2_down"] - dos["Ni_1_dx2y2_down"], "Ni 1 eg down")+
py4vasp.plot(dos["energies"], dos["Ni_2_dz2_up"] + dos["Ni_2_dx2y2_up"], "Ni 2 eg up")+
py4vasp.plot(dos["energies"], -dos["Ni_2_dz2_down"] - dos["Ni_2_dx2y2_down"], "Ni 2 eg down")
)
We see that the minority spin alternates between the Ni sites. And we see that, on each site, the majority spin is almost fully occupied and the minority spin almost fully unoccupied.
To see the resulting on-site magnetic moment for one Ni site, execute the following:
def Ni_moment(calc):
return sum(calc.local_moment.to_dict()['total'][0])
print("Ni moment: ", Ni_moment(mycalc))
How does it compare to the experimental value of about 1.9$\mu_B$? For an overview see Kwon and Min, PRB, 62, 73 (2000).
As a final step, you can visualize the magnetic structure!
mycalc.local_moment.plot()
3.4 Questions¶
- How can one select linear density mixing?
- Which tag sets the Fermi energy mid gap?
- Is the DOS of $t_{2g}$ orbitals directly computed in VASP?
4 Including on-site Coulomb interaction and pseudopotential selection for NiO relaxation ¶
By the end of this tutorial, you will be able to:
- perform a volume relaxation using a fit of an equation of state
- select the projector-augmented–wave (PAW) potential and compare calculated values using different PAWs
- run an LSDA+U calculation using the method by Dudarev et al. Phys. Rev. B 57, 1505 (1998)
4.1 Task¶
Fit the equation of state AFM NiO for U=8eV and J=0.95eV considering different PAW pseudopotentials and monitor the equilibrium lattice constant, the electronic bandgap, and the on-site magnetic moment.
The DFT+U method is an extension of density functional theory (DFT) designed to better treat strongly correlated electrons, especially in transition metal oxides and f-electron systems (like rare earths or actinides). The standard local density approximation (LDA) and generalized gradient approximation (GGA) often fails for such systems because it underestimates the on-site Coulomb repulsion (Hubbard U), which leads to incorrect predictions like metallic behavior for known Mott insulators. Traditionally the method is often referred to as LDA+U or LSDA+U although it can be combined with other exchange-correlation functionals, line GGA or even meta GGA.
This exercise is a quick interlude focussing on choosing the PAW potential and the DFT+U method which is also important for nonmagnetic calculations.
4.2 Input¶
The input files to run this example are partially prepared at $TUTORIALS/magnetism/e04_NiO_LSDA+U. Check them out!
POSCAR.tmp
AFM NiO
1
1.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 1.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 1.0000000000000000
Ni O
2 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
INCAR.tmp
SYSTEM = NiO
ISTART = 0
# Magnetism
ISPIN = 2
MAGMOM = 2.0 -2.0 2*0
LASPH = T
# FFT grid
ENCUT = 500
PREC = Acc
# Convergence criterion
EDIFF = 1E-8
# Smearing/k integration
ISMEAR = -5
EFERMI = MIDGAP
# Algorithm and density mixing
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
# compute bandgap
BANDGAP = WEIGHT
# compute on-site moment
LORBIT = 11
k points
0
gamma
6 6 6
0 0 0
POT_GGA_PAW_PBE/
Ni, Ni_GW, Ni_pv, Ni_sv_GW, and O, O_GW, O_GW_new, O_h, O_h_GW, O_s, O_s_GW
Based on this, let us create the input for the individual calculations. Compare different PAW pseudopotentials! Look in the directory POT_GGA_PAW_PBE and see what possible choices there are. Read about available PAW pseudopotentials and try deciding which one to use in this calculation. Consider the value of ENMAX and the atomic configuration underlying each potential.
We will use the following combination:
cd $TUTORIALS/magnetism/e04_*
mkdir PAW_NiO PAW_Ni_sv_GWO_s
cat POT_GGA_PAW_PBE/Ni/POTCAR POT_GGA_PAW_PBE/O/POTCAR > PAW_NiO/POTCAR
cat POT_GGA_PAW_PBE/Ni_sv_GW/POTCAR POT_GGA_PAW_PBE/O_s/POTCAR > PAW_Ni_sv_GWO_s/POTCAR
Why?
Click to see the answer!
There are 7 options for the oxygen potential: All of them correspond to the same atomic configuration. The first aspect is that they are distinct by their hardness. _s is particularly soft with a low energy cutoff, and _h is particularly hard, which is important in molecules with short bond lengths. A second aspect is that the potentials with _GW are particularly suited for calculations involving unoccupied states, e.g., for optical properties. Since here we investigate a bulk material (not very short bond length) where the main property of interest is magnetism, a ground state property emerging from the Ni 3d state, we can choose a computationally cheap O potential without significant loss of accuracy.
For Ni, there are 4 potentials, and some have a different atomic configuration. However, all of them include the 3d shell, and we do not expect the additional inclusion of p or s states to influence the magnetic properties. Again, magnetism is a ground state property and it is not necessary to opt for _GW.
Next, we need subdirectories for each PAW potential and U value. Execute the following in the Terminal:
for d in PAW_Ni_sv_GWO_s PAW_NiO
do
cd $d
mkdir U_0.0 U_3.0 U_8.0
cd ..
done
We write POSCAR files based on POSCAR.tmp centered around an educated guess for the equilibrium lattice constant for each PAW potential and U value. Execute the following:
import ase.io
pot = ["PAW_NiO","PAW_Ni_sv_GWO_s"]
u = [0.0, 3.0, 8.0]
a = {
"PAW_NiO": [4.064, 4.064, 4.035],
"PAW_Ni_sv_GWO_s": [4.067, 4.085, 4.119]
}
scale = [0.950, 0.975, 1.000, 1.025, 1.050]
s=ase.io.read("./e04_NiO_LSDA+U/POSCAR.tmp")
for p in pot:
for i in range(3):
for l in range(5):
# structure
s_out=s.copy()
s_out.set_cell(s.cell*a[p][i]*scale[l],scale_atoms=True)
# to file
filename="./e04_NiO_LSDA+U/"+p+"/U_{:.1f}/POSCAR.".format(u[i])+str(l)
ase.io.write(filename,s_out)
Copy the POTCAR files in the appropriate subdirectories!
cd $TUTORIALS/magnetism/e04_*
for d in PAW_Ni_sv_GWO_s PAW_NiO
do
for u in U_0.0 U_3.0 U_8.0
do
cp $d/POTCAR $d/$u/.
done
done
The INCAR files are all based on INCAR.tmp but differ in the block that defined the exchange-correlation functional. Here, we choose LDA by setting XC = CA. Then, various values for the on-site Coulomb repulsion are added. Finally, these can be copied to the subdirectories as follows:
cd $TUTORIALS/magnetism/e04_*
cp INCAR.tmp INCAR.U_0.0
echo """
# Exchange-correlation functional
XC = CA
""" >> INCAR.U_0.0
cp INCAR.tmp INCAR.U_3.0
echo """
# Exchange-correlation functional
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1 # l quantum number where U is added for each atom; -1 is no U added
LDAUU = 3.00 0.00 # on-site Coulomb interaction (in eV) for each atom
LDAUJ = 0.95 0.00 # on-site exchange interaction (in eV) for each atom
""" >> INCAR.U_3.0
cp INCAR.tmp INCAR.U_8.0
echo """
# Exchange-correlation functional
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1 # l quantum number where U is added for each atom; -1 is no U added
LDAUU = 8.00 0.00 # on-site Coulomb interaction (in eV) for each atom
LDAUJ = 0.95 0.00 # on-site exchange interaction (in eV) for each atom
""" >> INCAR.U_8.0
for d in PAW_Ni_sv_GWO_s PAW_NiO
do
for u in U_0.0 U_3.0 U_8.0
do
cp INCAR.$u $d/$u/INCAR
cp KPOINTS $d/$u/.
cp run $d/$u/.
done
done
To get an overview of the subdirectories in which we will run the calculations, execute the following in the Terminal:
cd $TUTORIALS/magnetism/e04_*
ls PAW_Ni*/U*
4.3 Calculation¶
Since this calculation takes quite some time, we recommend skipping the calculation described in Section 4.3.1 and using the data prepared in Section 4.3.2 instead!
Alternatively, you can continue with Section 4.3.1 and read about DFT+U on the VASP Wiki. You may also check the references therein. Particularly, S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
4.3.1 Performing the calculation¶
Open a terminal, navigate to this example's directory, run VASP for the two different POTCAR files, 3 U values and 4 different lattice constants, respectively. Execute the following:
cd $TUTORIALS/magnetism/e04_*
for d in PAW_Ni_sv_GWO_s PAW_NiO
do
for u in U_0.0 U_3.0 U_8.0
do
cd $d/$u
vasp_rm
bash run
cd ../..
done
done
where run contains the following instructions:
for i in 0 1 2 3 4
do
cp POSCAR.$i POSCAR
mpirun -np 4 vasp_std
cp vaspout.h5 vaspout_$i.h5
cp OUTCAR OUTCAR_$i
vasp_rm
done
Shortly after you started the calculation you can open $TUTORIALS/magnetism/e04_NiO_LSDA+U/NiO/U_0.0/OUTCAR_2 look for fundamental gap. Read the documentation of the BANDGAP tag and print the same data using the Python script below:
import py4vasp as pv
calc = pv.Calculation.from_file("e04_NiO_LSDA+U/PAW_NiO/U_0.0/vaspout_2.h5")
calc.bandgap.print()
The following Python script uses some helper functions to extract the total energy, on-site magnetic moment and fundamental bandgap for the set of parameters.
import py4vasp as pv
import numpy as np
import pandas as pd
pot = ["PAW_NiO","PAW_Ni_sv_GWO_s"]
u = [0.0, 3.0, 8.0]
a = {
"PAW_NiO": [4.064, 4.064, 4.035],
"PAW_Ni_sv_GWO_s": [4.067, 4.085, 4.119]
}
scale = [0.950, 0.975, 1.000, 1.025, 1.050]
def get_energy(i,l,p):
c=pv.Calculation.from_file("./e04_NiO_LSDA+U/"+p+"/U_{:.1f}/vaspout_{:}.h5".format(u[i],l))
return c.energy[-1].to_dict()['energy(sigma->0)']
def get_volume(i,l,p):
b=a[p][i]*scale[l]
print(b, b**3)
return (a[p][i]*scale[l])**3
def get_lattice_constant(i,l,p):
return a[p][i]*scale[l]
def get_U(i,l,p):
return u[i]
def get_fundamental(i,l,p):
c=pv.Calculation.from_file("./e04_NiO_LSDA+U/"+p+"/U_{:.1f}/vaspout_{:}.h5".format(u[i],l))
return c.bandgap.read()['fundamental']
def get_moment(i,l,p):
c=pv.Calculation.from_file("./e04_NiO_LSDA+U/"+p+"/U_{:.1f}/vaspout_{:}.h5".format(u[i],l))
return round(np.sum(c.local_moment.read()['total'][0]),3)
data = {
"lattice constants": [get_lattice_constant(i,l,p) for i in range(3) for l in range(5) for p in pot],
"U (in eV)": [get_U(i,l,p) for i in range(3) for l in range(5) for p in pot],
"total energy (in eV)": [get_energy(i,l,p) for i in range(3) for l in range(5) for p in pot],
"PAW potential": [p for i in range(3) for l in range(5) for p in pot ],
"fundamental bandgap (in eV)": [get_fundamental(i,l,p) for i in range(3) for l in range(5) for p in pot],
"on-site moment (in muB)": [get_moment(i,l,p) for i in range(3) for l in range(5) for p in pot]
}
df = pd.DataFrame(data)
df.to_csv('e04_NiO_LSDA+U/mydata.csv', index=False)
4.3.2 Skipping the calculation¶
Execute the following to read the data from file instead of computing it as described in the previous section.
import pandas as pd
df = pd.read_csv('e04_NiO_LSDA+U/data.csv')
When you prepared the whole setup but skipped Section 4.3.1, you can execute the following in the Terminal in order to run a single calculation
cd $TUTORIALS/magnetism/e04_*
cd PAW_NiO/U_8.0
cp POSCAR.2 POSCAR
mpirun -np 4 vasp_std
cp vaspout.h5 vaspout_2.h5
cp OUTCAR OUTCAR_2
Then, you can open $TUTORIALS/magnetism/e04_NiO_LSDA+U/NiO/U_8.0/OUTCAR_2 and look for fundamental gap. Read the documentation of the BANDGAP tag and print the same data using the Python script below:
import py4vasp as pv
calc = pv.Calculation.from_file("e04_NiO_LSDA+U/PAW_NiO/U_8.0/vaspout_2.h5")
calc.bandgap.print()
4.3.3 Analysis¶
Let us fit the Birch-Murnaghan equation of state! df_extended is a dataframe that will hold both the setup and the analysis.
from ase.eos import EquationOfState
extended_dfs = []
for p, group in df.groupby('PAW potential'):
for u_value, g in group.groupby('U (in eV)'):
volumes = g['lattice constants'].values**3
energies = g['total energy (in eV)'].values
eos = EquationOfState(volumes, energies, eos='birchmurnaghan')
v0, e0, B = eos.fit()
# Add results to a copy of the group
g_extended = g.copy()
g_extended['equilibrium lattice constant (in Å)'] = v0**(1/3)
g_extended['equilibrium volume (in ų)'] = v0
g_extended['minimum energy (in eV)'] = e0
g_extended['bulk modulus in (eV/ų)'] = B
# Append to the list
extended_dfs.append(g_extended)
# Combine all groups back into one DataFrame
df_extended = pd.concat(extended_dfs, ignore_index=True)
df_extended
Using the scripts below, plot the equilibrium lattice constant, magnetic moment and bandgap as a function of the U value for both choices of the PAW potential. Do all quantities depend on the choice of the PAW potential? Is the dependence only present in DFT+U or also for the plain LSDA calculation?
import matplotlib.pyplot as plt
plot_data = {}
for p, group in df_extended.groupby('PAW potential'):
us, mags, gaps, elat = [], [], [], []
for u_value, g in group.groupby('U (in eV)'):
# Get the equilibrium lattice constant for this group
a_eq = g['equilibrium lattice constant (in Å)'].iloc[0]
# Find the row with lattice constant closest to equilibrium
idx = (g['lattice constants'] - a_eq).abs().idxmin()
closest_row = g.loc[idx]
# Record U and magnetic moment
us.append(u_value)
mags.append(closest_row["on-site moment (in muB)"])
gaps.append(closest_row['fundamental bandgap (in eV)'])
elat.append(g['equilibrium lattice constant (in Å)'].values[-1])
# Store in dict by PAW potential
plot_data[p] = (us, mags, gaps, elat)
plt.figure(figsize=(8, 5))
for p, (us, mags, gaps, elat) in plot_data.items():
plt.plot(us, elat, marker='o', label=p)
plt.xlabel("U value (in eV)")
plt.ylabel("Equilibrium lattice constant (in Å)")
plt.title("Equilibrium lattice constant vs U value for each PAW Potential")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
The equilibrium lattice constant shows a dramatically different behaviour as a function of U for the two PAW potentials. For the Ni_sv_GW the Ni p and s orbitals are both in the valence. These hybridize with the Ni d orbitals. For increasing U values the d states are pushed away from the Fermi energy and decrease in size. Thus, the entangled sp orbitals expand and the lattice constant increases.
On the other hand, the simple Ni PAW potential holds the sp orbitals in the frozen core. Hence, they cannot expand and the lattice constant does not increase. Instead the lattice constant even slightly decreases due to the shrinking d orbitals.
plt.figure(figsize=(8, 5))
for p, (us, mags, gaps, elat) in plot_data.items():
plt.plot(us, mags, marker='o', label=p)
plt.xlabel("U value (in eV)")
plt.ylabel("Magnetic moment (in muB)")
plt.title("Magnetic moment vs U value for each PAW Potential")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
Here, we see that increasing the on-site Coulomb interaction can have a significant effect on the on-site magnetic moments at the equilibrium lattice constant. Yet, since the magnetic moment arises in the Ni d orbitals, the effect of the selected PAW potentials is very small.
plt.figure(figsize=(8, 5))
for p, (us, mags, gaps, elat) in plot_data.items():
plt.plot(us, gaps, marker='o', label=p)
plt.xlabel("U value (in eV)")
plt.ylabel("Bandgap (in eV)")
plt.title("Fundamental bandgap vs U value for each PAW Potential")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
We finally see that the choice of the PAW potential has vanishing effects on the electronic bandgap at the equilibrium lattice constant as a function of U. This emphasizes that one has to consider multiple aspects when selecting a PAW potential. In this example, the DFT+U method plays a central role, as well as the quantity of interest. The most accurate results are always obtained with the more expensive PAW potentials, which may, e.g., treat more electrons as valence, or use a higher cutoff energy. However, if the quantity of interest is not sensitive to the choice of the PAW potential it is better to opt for a cheaper pseudopotential to avoid unnecessary computational cost.
Also note that the quantity of interest may depend on the volume. Many publications use the experimental lattice constant, but you should keep in mind that the equilibrium lattice constant depends amongst other aspects on the selected exchange-correlation functional. See below how the bandgap changes with the lattice constant.
import plotly.express as px
for p, group in df.groupby('PAW potential'):
fig = px.line(
group,
x='lattice constants',
y='fundamental bandgap (in eV)',
color='U (in eV)',
markers=True,
title=p+': Bandgap vs lattice constant for different U values'
)
fig.update_layout(
xaxis_title='Lattice constant (Å)',
yaxis_title='Fundamental bandgap (in eV)',
legend_title='U (eV)',
template='plotly_white'
)
fig.show()
4.4 Questions¶
- What is an equation of state?
- Is there a particular quantity that the choice of the PAW potential should be based on?
- What options does the LDAUTYPE tag offer?