Part 2: Coupling constants and two-center corrections ¶
By the end of this tutorial, you will be able to:
- Perform a hyperfine coupling calculation
- Compare the hyperfine coupling parameter to experiment and literature
- See how the hypferfine coupling changes by using a GGA or a hybrid DFT functional
5.1 Task¶
Calculate the hyperfine coupling tensor for a gas-phase CH3$^\bullet$ radical in 10x10x10 Å3 cell using PBE and HSE06
The nuclei are not the only part of an atom that have spin, the electrons do too. The interaction between internally generated magnetic field (from the electrons) and the magnetic dipole moment of the nucleus result in the splitting of otherwise degenerate energy levels. This splitting is known as hyperfine splitting.
The Hamiltonian $\hat{H}_D$ describing this hyperfine splitting is defined for a nuclear magnetic dipole moment $\boldsymbol{\mu}_I$ in a magnetic field $\boldsymbol{B}$:
$$\tag{1} \begin{equation} \hat{H}_D = - \boldsymbol{\mu}_I \cdot \boldsymbol{B}. \end{equation} $$
In this case, the magnetic field is an internal one, described by contributions from the electronic orbital $\boldsymbol{B}_{el}^l$ and spin $\boldsymbol{B}_{el}^s$ angular momenta:
$$\tag{2} \begin{equation} \boldsymbol{B} = \boldsymbol{B}_{el}^l + \boldsymbol{B}_{el}^s. \end{equation} $$
The nuclear spin $S^I$ and the electronic spin $S^e$ are connected by the hyperfine tensor $A^I$ to give the energy of the hyperfine splitting [Phys. Rev. 88 (2013) 075202]:
$$\tag{3} \begin{equation} E = \sum_{ij} S^e_i A^I_{ij} S^I_j. \end{equation} $$
In this exercise, you will calculate the hyperfine tensor $A^I$ for the methyl (CH3$^\bullet$) radical using PBE and then compare this to using HSE06.
5.2 Input¶
CH3 radical
1.0
10.0000000000 0.0000000000 0.0000000000
0.0000000000 10.0000000000 0.0000000000
0.0000000000 0.0000000000 10.0000000000
C H
1 3
Cartesian
5.129648629 4.775283014 5.091733012
14.947834482 5.090380970 6.116588891
6.035344916 5.090363074 14.578585301
14.404303973 4.148275943 14.578762796
INCAR.pbe
ENCUT = 500 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-6 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LHYPERFINE = .TRUE. # Turns on calculating the hyperfine coupling tensor
NGYROMAG = 10.7084 42.577478461 # Specifies the nuclear gyromagnetic ratios for the ions
ISPIN = 2 # Turns on spin-polarization
LASPH = .TRUE. # Includes non-spherical contributions to the gradient inside the PAW spheres
INCAR.hse06
ENCUT = 500 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-6 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LHYPERFINE = .TRUE. # Turns on calculating the hyperfine coupling tensor
NGYROMAG = 10.7084 42.577478461 # Specifies the nuclear gyromagnetic ratios for the ions
ISPIN = 2 # Turns on spin-polarization
LASPH = .TRUE. # Includes non-spherical contributions to the gradient inside the PAW spheres
LHFCALC = .TRUE. # Sets the density functional to HSE06
GGA = PE
HFSCREEN = 0.2
ALGO = Damped # Sets a damped velocity friction algorithm - recommended for use with LHFCALC
Gamma-point k-mesh
0
Gamma
1 1 1
0 0 0
Pseudopotential of C and H.
The INCAR file is similar to that used in Example 1, with a few additional tags: LHYPERFINE = .TRUE. sets the hyperfine coupling tensor to be calculated, while NGYROMAG provides the nuclear gyromagnetic ratios of the ions in order to calculate the hyperfine coupling. In this case, 13C and 1H have been chosen. Since there is an unpaired electron, ISPIN = 2 is used which enables spin-polarized calculations.
The remaining tags in INCAR.hse06: LHFCALC = .TRUE., GGA = PE, and HFSCREEN = 0.2 set up the HSE06 functional. ALGO = Damped sets a damped velocity friction algorithm to be used, which is recommended for use with LHFCALC.
The KPOINTS file specifies a Γ-point $\mathbf{k}$-mesh.
5.3 Calculation¶
Open a terminal and navigate to this example's directory:
cd $TUTORIALS/NMR/e05_*
Try running doing a test-run to see what a hyperfine calculation looks like. Copy the INCAR.pbe file to INCAR and run it with the following:
cp INCAR.pbe INCAR
mpirun -np 4 vasp_std
Congratulations on calculating the hyperfine coupling parameter! The hyperfine coupling parameter is split into an isotropic part, the Fermi contact term, and an anistropic part, the dipole contributions. You can find it printed out in the OUTCAR file after Fermi contact (isotropic) hyperfine coupling parameter (MHz):
Fermi contact (isotropic) hyperfine coupling parameter (MHz)
-------------------------------------------------------------
ion A_pw A_1PS A_1AE A_1c A_tot
-------------------------------------------------------------
1 ... ... ... ... ...
Followed by the contributions for each ion listed in the POSCAR. These contributions in MHz are broken down in the table below:
| tag | meaning |
|---|---|
A_pw |
plane-wave contribution to the Fermi contact term |
A_1PS |
pseudo one-center contribution to the Fermi contact term |
A_1AE |
all-electron one-center contribution to the Fermi contact term |
A_1c |
one-center core contribution to the Fermi contact term [Phys. Rev. 71 (2005) 115110] |
A_tot |
total contribution to the Fermi contact term, omitting A_1c |
The hyperfine output continues after the line where the dipolar contributions to the hyperfine tensor $A_{ij}$ are given:
Dipolar hyperfine coupling parameters (MHz)
---------------------------------------------------------------------
ion A_xx A_yy A_zz A_xy A_xz A_yz
---------------------------------------------------------------------
1 ... ... ... ... ... ...
Finally, the total hyperfine tensors are written out after diagonalization:
Total hyperfine coupling parameters after diagonalization (MHz)
(convention: |A_zz| > |A_xx| > |A_yy|)
----------------------------------------------------------------------
ion A_xx A_yy A_zz asymmetry (A_yy - A_xx)/ A_zz
----------------------------------------------------------------------
1 ... ... ... ...
Now that you have a sense of what to look for in the OUTCAR file after the hyperfine calculation. Notice that the three H protons (ions 2-4) have near identical values. Take a look at the structure to see if you can work out why:
import py4vasp
struc = "./e05_hyperfine/"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
The CH3$^\bullet$ radical has three-fold symmetry, so the 1H nuclei are equivalent.
Enter the next into the terminal to find the hyperfine coupling parameters in the output:
cd $TUTORIALS/NMR/e05_*
en=$(awk '/free e/ {print $5}' OUTCAR)
hyperfine=$(grep -A 2 A_1AE OUTCAR | tail -1 | awk '{print $6 " " $5}')
echo $en $hyperfine >> hyperfine.dat
We will focus on the 13C nucleus. Take a look at the Fermi contact (isotropic) hyperfine coupling parameter in the OUTCAR. The total hyperfine coupling parameter $A_{tot}$ and the core contributions $A_{1c}$ have been saved to hyperfine.dat. $A_{tot}$ does not include $A_{1c}$ so they need to be added in explicitly. Using the following script, print out $A_{tot}$ (line 1), $A_{1c}$ (line 2), and their sum $A_{tot} + A_{1c}$ (line 3):
import numpy as np
hf = np.genfromtxt("./e05_hyperfine/hyperfine.dat",unpack=True,usecols=range(3))
print("Converged A_tot: " + str(np.round(hf[1], 1)) + " MHz")
print("Converged A_1c: " + str(np.round(hf[2], 1)) + " MHz")
print("Converged A_tot with A_1c: " + str(np.round((hf[1] + hf[2]), 1)) + " MHz")
Though not generally the case, the hyperfine coupling parameter converges rapidly with respect to the energy cutoff, converging by 500 eV, in stark contrast to the chemical shielding. In solid-state applications, the $\textbf{k}$-mesh would also need to be converged, e.g. for the negatively charged nitrogen-vacancy (NV) defect in diamond [Phys. Rev. B 88 (2013) 075202]. Given that the coupling parameters are already converged, we will take a look at the influence of the method on the parameter.
You can compare your total hyperfine coupling parameter $A_{tot}$ and the core contributions $A_{1c}$ to literature and experimental values. The experimental reference is measured using electron paramagnetic resonance (EPR), the unpaired electron analog to NMR. In the above literature reference, $A_{tot}$ is calculated to be 182.5 MHz. How do your results compare?
The difference between your value and the literature is explained by the inclusion of LASPH = .TRUE., which was not included in the above reference. LASPH turns on the non-spherical contributions to the gradient inside the PAW spheres. If you were to set LASPH = .FALSE. then your results would match the literature except for a small difference due to their use of a larger cell for the molecule, $20 \times 20 \times 20$ Å3, compared to your $10 \times 10 \times 10$ Å3 cell.
The experimental value is 107.4 MHz. The experimental is very different from the calculated value of $A_{tot}$, due to the exclusion of the core contributions. Once the core is added in ($A_{tot} + A_{1c}$), the literature value changes to 86.1 MHz, much closer to experiment. Your value should be ~70 MHz.
However, the remaining difference of 20-30 MHz is still significant. Let's see how using a hybrid functional HSE06 changes things [J. Chem. Phys. 125 (2006) 224106].
Take the INCAR.hse06 file and submit an HSE06 calculation using the following:
cp INCAR.hse06 INCAR
mpirun -np 4 vasp_std
Enter the next into the terminal to find the hyperfine coupling parameters in the output:
cd $TUTORIALS/NMR/e05_*
en=$(awk '/free e/ {print $5}' OUTCAR)
hyperfine=$(grep -A 2 A_1AE OUTCAR | tail -1 | awk '{print $6 " " $5}')
echo $en $hyperfine >> hyperfine.dat
How does HSE06 compare to experiment? Use the following to print out the hyperfine coupling parameter:
import numpy as np
hf = np.genfromtxt("./e05_hyperfine/hyperfine.dat",unpack=True,usecols=range(3))
print("Converged A_tot: " + str(np.round(hf[1][1], 1)) + " MHz")
print("Converged A_1c: " + str(np.round(hf[2][1], 1)) + " MHz")
print("Converged A_tot with A_1c: " + str(np.round((hf[1][1] + hf[2][1]), 1)) + " MHz")
The literature value of $A_{tot} + A_{1c}$ was 101.9 MHz. How closely does yours match?
Can you think of a reason why your value differs?
The difference between your value and the literature is explained by the inclusion of LASPH = .TRUE., which was not included in the above reference. LASPH turns on the non-spherical contributions to the gradient inside the PAW spheres. If you were to set LASPH = .FALSE. then your results would match the literature except for a small difference due to their use of a larger cell for the molecule, $20 \times 20 \times 20$ Å3, compared to your $10 \times 10 \times 10$ Å3 cell.
Compared to the experimental value of 107.4 MHz, it is clear that HSE06 offers a significant improvement over PBE.
Can you think of why HSE06 offers so much improvement?
Hybrid functionals include exact exchange, which helps to localize lone electrons or defects [J. Chem. Phys. 125 (2006) 224106].
5.4 Questions¶
- What term should be added to $A_{tot}$ to improve comparison with experiment?
- How can you define the nuclear gyromagnetic ratios?
- Why do hybrid functionals offer an improvement over GGA when calculating the coupling constant?
By the end of this tutorial, you will be able to:
- Calculate the electric field gradient (EFG)
- Converge the quadrupolar coupling constant by tightening the energetic break condition for the self-consistency loop
6.1 Task¶
Calculate the nuclear quadrupole moments for tetragonal-phased methylammonium lead iodide using a $\Gamma$-only $\textbf{k}$-point mesh and test impact of the energetic break condition for the self-consistency loop on the quadrupolar moments
The electric field gradient (EFG) is the rate of change of the electric field at the nucleus. Quadrupolar nuclei (with spin $\gt \frac{1}{2}$) have a quadrupolar electric moment, originating from non-spherical nuclei, which couples together with the EFG. The EFG is the second derivative of the electrostatic potential $V$ at the position of the nucleus $x$:
$$\tag{4} \begin{equation} V_{ij} = \frac{\partial^2 V}{\partial x_i \partial x_j}. \end{equation} $$
The EFG itself is not measurable, but the nuclear quadrupolar coupling constant $C_q$ is, which may be measured using various spectroscopic techniques such as NMR, EPR, or nuclear quadrupole resonance (NQR) spectroscopy. $C_q$ is defined as:
$$\tag{5} \begin{equation} C_q = \frac{eQV_{zz}}{h}. \end{equation} $$
where $e$ is the charge of an electron, $Q$ is the element and isotope specific quadrupole moment, and $h$ is the Planck constant.
The nuclear quadrupolar moments are important for methylammonium lead iodide MAPbI3, a perovskite used in solar cells [J. Phys. Chem. Lett. 8 (2017) 61]. In this exercise, you will calculate the nuclear quadrupole moments for MAPbI3 for I and N and see how they are impacted by tightening the energetic break condition for the self-consistency loop (EDIFF).
6.2 Input¶
Click to show the POSCAR!
Pnma Weller
1.0000000000000000
8.8657400000000006 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.629300000000001 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.5768900000000006
Pb I N C H
4 12 4 4 24
Direct
0.00000000000000 0.00000000000000 0.50000000000000
0.50000000000000 0.00000000000000 0.00000000000000
0.00000000000000 0.50000000000000 0.50000000000000
0.50000000000000 0.50000000000000 0.00000000000000
0.48420000000000 0.25000000000000 -0.05620000000000
0.01580000000000 0.75000000000000 0.44380000000000
-0.48420000000000 0.75000000000000 0.05620000000000
0.98420000000000 0.25000000000000 0.55620000000000
0.18860000000000 0.01470000000000 0.18440000000000
0.31140000000000 -0.01470000000000 0.68440000000000
-0.18860000000000 0.51470000000000 -0.18440000000000
0.68860000000000 0.48530000000000 0.31560000000000
-0.18860000000000 -0.01470000000000 -0.18440000000000
0.68860000000000 0.01470000000000 0.31560000000000
0.18860000000000 0.48530000000000 0.18440000000000
0.31140000000000 0.51470000000000 0.68440000000000
-0.94210000000000 0.25000000000000 -0.02970000000000
1.44210000000000 0.75000000000000 0.47030000000000
0.94210000000000 0.75000000000000 0.02970000000000
-0.44210000000000 0.25000000000000 0.52970000000000
0.43720000000000 0.25000000000000 0.44250000000000
0.06280000000000 0.75000000000000 0.94250000000000
-0.43720000000000 0.75000000000000 -0.44250000000000
0.93720000000000 0.25000000000000 0.05750000000000
0.93720000000000 0.25000000000000 0.18740000000000
-0.43720000000000 0.75000000000000 0.68740000000000
-0.93720000000000 0.75000000000000 -0.18740000000000
1.43720000000000 0.25000000000000 0.31260000000000
0.86610000000000 0.17010000000000 0.02900000000000
-0.36610000000000 -0.17010000000000 0.52900000000000
-0.86610000000000 0.67010000000000 -0.02900000000000
1.36610000000000 0.32990000000000 0.47100000000000
-0.86610000000000 -0.17010000000000 -0.02900000000000
1.36610000000000 0.17010000000000 0.47100000000000
0.86610000000000 0.32990000000000 0.02900000000000
-0.36610000000000 0.67010000000000 0.52900000000000
0.12750000000000 0.18910000000000 -0.00850000000000
0.37250000000000 -0.18910000000000 0.49150000000000
-0.12750000000000 0.68910000000000 0.00850000000000
0.62750000000000 0.31090000000000 0.50850000000000
-0.12750000000000 -0.18910000000000 0.00850000000000
0.62750000000000 0.18910000000000 0.50850000000000
0.12750000000000 0.31090000000000 -0.00850000000000
0.37250000000000 0.68910000000000 0.49150000000000
-0.95430000000000 0.25000000000000 -0.14590000000000
1.45430000000000 0.75000000000000 0.35410000000000
0.95430000000000 0.75000000000000 0.14590000000000
-0.45430000000000 0.25000000000000 0.64590000000000
INCAR.efg
ENCUT = 400 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-4 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LEFG = .TRUE. # Electric field gradient calculations
QUAD_EFG = 0. -696. 20.44 0. 2.860 # Nuclear quadrupolar moments for Pb I N O D
Gamma-point k-mesh
0
Monkhorst-Pack scheme
1 1 1
0 0 0
Pseudopotential of Pb, I, N, C, and H.
Check the tags set in the INCAR file! The INCAR tags are identical to those used in Example 1, with the addition of two more tags:
LEFG enables the electric field gradient calculations and QUAD_EFG defines the nuclear quadrupole moments (in millbarn) for the atomic types on the POTCAR file.
The KPOINTS file specifies a $\Gamma$-only $\mathbf{k}$-point mesh.
6.3 Calculation¶
Open a terminal and navigate to this example's directory by entering the following:
cd $TUTORIALS/NMR/e06_*
You will change the energetic convergence parameter for the electronic step, EDIFF, and see how the calculated electric field gradient changes depending on it. In order to run VASP using different convergence settings, the INCAR file needs to be adjusted, EDIFF = 1E-4 eV is replaced with: 1E-5, 1E-6, 1E-7, and 1E-8 eV, which is done by the ediff_efg.sh. In this example's directory enter the following into the terminal to run the script:
bash ediff_efg.sh
ediff_efg.sh runs the calculations and finds the quadrupole coupling constant for Cq in the OUTCAR file using the following grep command:
grep -A 30 'Cq(MHz)' OUTCAR
Click to see the provided script to run the jobs and extract the quadrupolar coupling constants!
rm -f ediff_efg.dat
for ediff in 4 5 6 7 8
do
cp INCAR.efg INCAR
sed -i "s/EDIFF = 1E-4/EDIFF = 1E-$ediff/g" INCAR
mpirun -np 4 vasp_std
Cq=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $2}')
eta=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $3}')
CqI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $2}')
etaI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $3}')
CqI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $2}')
etaI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $3}')
echo $ediff $Cq $eta $CqI1 $etaI1 $CqI2 $etaI2 >> ediff_efg.dat
done
Congratulations on running your first EFG calculation!
The quadrupole coupling constants from multiple VASP calculations at different energy cutoffs have been stored in ediff_efg.dat using the ediff_efg.sh script! The electric field gradient is found in the OUTCAR file:
Electric field gradients after diagonalization (V/A^2)
(convention: |V_zz| > |V_xx| > |V_yy|)
----------------------------------------------------------------------
ion V_xx V_yy V_zz asymmetry (V_yy - V_xx)/ V_zz
----------------------------------------------------------------------
1 --- --- --- ---
2 --- --- --- ---
----------------------------------------------------------------------
Where V_ii is the electric field gradient component along the i axis.
The quadrupolar coupling constants are found further down in the OUTCAR file:
----------------------------------------------------------------------
NMR quadrupolar parameters
Cq : quadrupolar parameter Cq=e*Q*V_zz/h
eta: asymmetry parameters (V_yy - V_xx)/ V_zz
Q : nuclear electric quadrupole moment in mb (millibarn)
----------------------------------------------------------------------
ion Cq(MHz) eta Q (mb)
----------------------------------------------------------------------
1 --- --- ---
2 --- --- ---
----------------------------------------------------------------------
Using py4vasp, take a look at the structure of MAPbI3:
import py4vasp
struc = "./e06_efg_ediff/"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
I: purple
Pb: black
N: blue
C: grey
H: white
There are two types of I atoms in the cell, those with two I in a layer, axial to the Pb, (in the same layer as the CH3NH3 molecules, methyl amine MA) and those between the methyl amine layers, with four I in a layer, equatorial to the Pb. Let's plot the quadrupole coupling constants for each of these and the N in the methyl amine to see how they vary with different EDIFF values. First, print out the Cq values:
import numpy as np
import pandas as pd
cq = np.genfromtxt("./e06_efg_ediff/ediff_efg.dat", unpack=True, usecols=range(7))
df = pd.DataFrame(cq).transpose()
df.rename(columns={df.columns[0]: "log(EDIFF in eV)", df.columns[1]: "C_q(N) in MHz", df.columns[3]: "C_q(axial I) in MHz", df.columns[5]: "C_q(equatorial I) in MHz"}, inplace = True)
df = df.set_index(list(df)[0])
df = df.iloc[:, 0:6:2]
df
Then plot them using py4vasp:
import py4vasp
import numpy as np
import plotly.graph_objs as go
cq = np.genfromtxt("./e06_efg_ediff/ediff_efg.dat", unpack=True, usecols=range(7))
py4vasp.plot(cq[0], cq[3], y2=True, xlabel="log(EDIFF in eV)", y2label="Quadrupole coupling constant in MHz", label="Axial I (right y-axis)")\
+ py4vasp.plot(cq[0], cq[5], ylabel="Quadrupole coupling constant in MHz", label="Equatorial I (left y-axis)")
You can see that the quadrupolar coupling constants of Axial I and Equatorial I change a lot going from a self-consistency loop break condition of 10-4 eV to 10-5 eV. After 10-5 eV, it changes by less than 0.1 MHz, and convergence is considered achieved. N is already converged using 10-4 eV, cf. the table above. The degree to which the coupling constant (or more generally, your property of interest) changes in response to a technical parameter will depend a lot on your indiviudal system and must be tested for each; the EFG is often strongly affected by the setting of EDIFF. These sorts of convergence tests are necessary to perform for any new set of calculations that you do, here the EFG and the coupling constant to describe the quadrupolar electric field of spin $\gt \frac{1}{2}$ nuclei.
Note, though not explored here, the choice of POTCAR can have a significant impact on the value of the quadrupolar coupling constant and should be considered in real calculations. For example, GW pseudopotentials and additional semi-core electrons (i.e., *_pv and *_sv POTCAR) can be important. Before converging the quadrupolar coupling constant, you should check whether or not the POTCAR has a large impact.
Having tested the importance of the SCF convergence parameter, in the next exercise, you will converge the quadrupolar coupling constants with respect to the energy cutoff and the $\textbf{k}$-point mesh.
By the end of this tutorial, you will be able to:
- Converge the electric field gradient with respect to k-points and plane-wave energy cutoff
- Determine suitable convergence limits for the quadrupolar coupling constant for different elements
7.1 Task¶
Calculate the nuclear quadrupole moments for tetragonal-phased methylammonium lead iodide, converging with respect to $\textbf{k}$-point mesh and plane-wave energy cutoff.
In methylammonium lead iodide MAPbI3, there are two elements with high proportions of quadrupolar nuclei: 14N (99.9(6) % abundant) and 127I (100 % abundant). Experiments such as nuclear quadrupole resonance spectroscopy (NQR) allow for the observation of these nuclei in the solid state.
There are two distinct types of I nuclei arranged in octahedra around Pb in MAPb3, the axial, those in plane with the Pb atoms (4 per cell) and the equatorial, those in plane with the methylamine molecules (2 per cell). Together, the I form an octahedron around the Pb atoms. In NQR, these two different I environments can be detected [J. Phys. Chem. Lett. 8 (2017) 61]. The N, however, are in very similar environments, so are difficult to distinguish by experiment.
In this exercise, you will continue on from Example 6 to converge the quadrupolar coupling constant for 14N and 127I with respect to $\textbf{k}$-point mesh and plane-wave energy cutoff, comparing to experiment where possible.
7.2 Input¶
Click to show the POSCAR!
Pnma Weller
1.0000000000000000
8.8657400000000006 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.629300000000001 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.5768900000000006
Pb I N C H
4 12 4 4 24
Direct
0.00000000000000 0.00000000000000 0.50000000000000
0.50000000000000 0.00000000000000 0.00000000000000
0.00000000000000 0.50000000000000 0.50000000000000
0.50000000000000 0.50000000000000 0.00000000000000
0.48420000000000 0.25000000000000 -0.05620000000000
0.01580000000000 0.75000000000000 0.44380000000000
-0.48420000000000 0.75000000000000 0.05620000000000
0.98420000000000 0.25000000000000 0.55620000000000
0.18860000000000 0.01470000000000 0.18440000000000
0.31140000000000 -0.01470000000000 0.68440000000000
-0.18860000000000 0.51470000000000 -0.18440000000000
0.68860000000000 0.48530000000000 0.31560000000000
-0.18860000000000 -0.01470000000000 -0.18440000000000
0.68860000000000 0.01470000000000 0.31560000000000
0.18860000000000 0.48530000000000 0.18440000000000
0.31140000000000 0.51470000000000 0.68440000000000
-0.94210000000000 0.25000000000000 -0.02970000000000
1.44210000000000 0.75000000000000 0.47030000000000
0.94210000000000 0.75000000000000 0.02970000000000
-0.44210000000000 0.25000000000000 0.52970000000000
0.43720000000000 0.25000000000000 0.44250000000000
0.06280000000000 0.75000000000000 0.94250000000000
-0.43720000000000 0.75000000000000 -0.44250000000000
0.93720000000000 0.25000000000000 0.05750000000000
0.93720000000000 0.25000000000000 0.18740000000000
-0.43720000000000 0.75000000000000 0.68740000000000
-0.93720000000000 0.75000000000000 -0.18740000000000
1.43720000000000 0.25000000000000 0.31260000000000
0.86610000000000 0.17010000000000 0.02900000000000
-0.36610000000000 -0.17010000000000 0.52900000000000
-0.86610000000000 0.67010000000000 -0.02900000000000
1.36610000000000 0.32990000000000 0.47100000000000
-0.86610000000000 -0.17010000000000 -0.02900000000000
1.36610000000000 0.17010000000000 0.47100000000000
0.86610000000000 0.32990000000000 0.02900000000000
-0.36610000000000 0.67010000000000 0.52900000000000
0.12750000000000 0.18910000000000 -0.00850000000000
0.37250000000000 -0.18910000000000 0.49150000000000
-0.12750000000000 0.68910000000000 0.00850000000000
0.62750000000000 0.31090000000000 0.50850000000000
-0.12750000000000 -0.18910000000000 0.00850000000000
0.62750000000000 0.18910000000000 0.50850000000000
0.12750000000000 0.31090000000000 -0.00850000000000
0.37250000000000 0.68910000000000 0.49150000000000
-0.95430000000000 0.25000000000000 -0.14590000000000
1.45430000000000 0.75000000000000 0.35410000000000
0.95430000000000 0.75000000000000 0.14590000000000
-0.45430000000000 0.25000000000000 0.64590000000000
INCAR.efg
ENCUT = 200 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-4 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LEFG = .TRUE. # Electric field gradient calculations
QUAD_EFG = 0. -696. 20.44 0. 2.860 # Nuclear quadrupolar moments for Pb I N O D
KPOINTS.efg
Gamma-point k-mesh
0
Monkhorst-Pack scheme
2 2 2
0 0 0
Pseudopotential of Pb, I, N, C, and H.
7.3 Calculation¶
Open a terminal and navigate to this example's directory:
cd $TUTORIALS/NMR/e07_*
We will first converge the chemical shielding with respect to the plane-wave energy cutoff. In order to run VASP using different energy cutoffs, the INCAR file needs to be adjusted. In the following script, ENCUT = 200 eV in INCAR is replaced with a new energy cutoff: 300 to 600 eV using a step-size of 100 eV (to save time, results from 500 eV and 600 eV are provided). Note you should never use an ENCUT value lower than the ENMAX defined in the POTCAR file; they are used in this example only to illustrate convergence.
In this example's directory, enter the following into the terminal to run the encut_efg.sh script:
bash encut_efg.sh
encut_efg.sh runs the calculations and finds the coupling constants for the N, the axial I, and the equatorial I in the OUTCAR file using the following grep commands:
grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $2}' # N
grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $2}' # Axial I
grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $2}' # Equatorial I
Click to see the provided script to run the jobs and extract the quadrupolar coupling constants!
rm -f encut_efg.dat
for encut in 200 300 400
do
cp INCAR.efg INCAR
sed -i "s/ENCUT = 200/ENCUT = $encut/g" INCAR
cp KPOINTS.efg KPOINTS
mpirun -np 4 vasp_std
Cq=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $2}')
eta=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $3}')
CqI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $2}')
etaI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $3}')
CqI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $2}')
etaI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $3}')
echo $encut $Cq $eta $CqI1 $etaI1 $CqI2 $etaI2 >> encut_efg.dat
done
echo "500 0.423 0.053 -548.255 0.331 -530.985 0.227" >> encut_efg.dat
echo "600 0.422 0.050 -548.273 0.331 -531.025 0.227" >> encut_efg.dat
To keep the structure in mind, visualize it in py4vasp using the following:
import py4vasp
struc = "./e07_efg_converge"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot() # take a look at the structure
I: purple
Pb: black
N: blue
C: grey
H: white
The quadrupolar coupling constants from multiple VASP calculations at different energy cutoffs have been stored in encut_efg.dat using the encut_efg.sh script!
Using the next two python scripts, visualize the quadrupolar coupling constants that you have calculated and then plot it at each cutoff using py4vasp!
import numpy as np
import pandas as pd
cq = np.genfromtxt("./e07_efg_converge/encut_efg.dat", unpack=True, usecols=range(7))
df = pd.DataFrame(cq).transpose()
df.rename(columns={df.columns[0]: "ENCUT in eV", df.columns[1]: "C_q(N) in MHz", df.columns[3]: "C_q(Axial I) in MHz", df.columns[5]: "C_q(Equatorial I) in MHz"}, inplace = True)
df = df.set_index(list(df)[0])
df = df.iloc[:, 0:6:2]
df
import py4vasp
import numpy as np
cq = np.genfromtxt("./e07_efg_converge/encut_efg.dat", unpack=True, usecols=range(7))
py4vasp.plot(cq[0], cq[3], y2=True, xlabel="Energy cutoff in eV", y2label="Quadrupolar coupling constant (I) in MHz", label="Axial I (right y-axis)") \
+ py4vasp.plot(cq[0], cq[5], y2=True, y2label="Quadrupolar coupling constant (I) in MHz", label="Equatorial I (right y-axis)")\
+ py4vasp.plot(cq[0], cq[1], ylabel="Quadrupolar coupling constant (N) in MHz", label="N (left y-axis)")
How quickly do the different quadrupolar coupling constants converge with respect to the energy cutoff?
Cq converges rapidly for the equatorial I, by 300 eV. The axial I converge at a similar rate. Note that this is to within ~0.1 MHz. Due to the magnitude, the N requires much tighter convergence settings, to within 0.01 MHz. This is not achieved until 400 eV and 500 eV for complete convergence. From these tests, we would recommend a cutoff of 400 eV as a minimum, preferably 500 eV.
The plane-wave energy cutoff is only one of the important parameters, alongside the $\textbf{k}$-point mesh. Here, we use a cutoff of 400 eV while testing the $\textbf{k}$-point mesh. In real applications, parameters should be converged sequentially, so a cutoff of 500 eV would be used to converge the $\textbf{k}$-point mesh.
In order to run VASP using different $\textbf{k}$-point meshes, the KPOINTS file needs to be updated. For each run, replace 2 2 2 in KPOINTS (where N N N corresponds to a NxNxN $\textbf{k}$-point mesh) with 1 1 1 and 2 2 2, in addition the results from 4 4 4 and 6 6 6 will be added in, though they are not calculated here due to time constraints. The script kpoints_efg.sh will set up and run these calculations for you, extract the quadrupolar coupling constants, and store them in kpoints_efg.dat. In this example's directory enter the following into the terminal:
bash kpoints_efg.sh
Click to see the provided script to run the jobs and extract the quadrupolar coupling constants!
rm -f kpoints_efg.dat
cp INCAR.efg INCAR
sed -i "s/ENCUT = 200/ENCUT = 400/g" INCAR
for kpt in 1 2
do
cp KPOINTS.efg KPOINTS
sed -i "s/2 2 2/$kpt $kpt $kpt/g" KPOINTS
mpirun -np 4 vasp_std
Cq=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $2}')
eta=$(grep -A 30 'Cq(MHz)' OUTCAR | grep '20.440' | tail -1 | awk '{ print $3}')
CqI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $2}')
etaI1=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | head -1 | awk '{ print $3}')
CqI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $2}')
etaI2=$(grep -A 30 'Cq(MHz)' OUTCAR | grep ' -696' | tail -1 | awk '{ print $3}')
echo $kpt $Cq $eta $CqI1 $etaI1 $CqI2 $etaI2 >> kpoints_efg.dat
done
echo "4 0.435 0.052 -528.959 0.342 -549.491 0.218" >> kpoints_efg.dat
echo "6 0.435 0.051 -528.579 0.342 -549.714 0.218" >> kpoints_efg.dat
Using the next two python scripts, visualize the quadrupolar coupling constants that you have calculated and then plot it at each $\textbf{k}$-mesh using py4vasp!
import numpy as np
import pandas as pd
cq = np.genfromtxt("./e07_efg_converge/kpoints_efg.dat", unpack=True, usecols=range(7))
df = pd.DataFrame(cq).transpose()
df.rename(columns={df.columns[0]: "No. k-points along lattice parameter", df.columns[1]: "C_q(N) in MHz", df.columns[3]: "C_q(Axial I) in MHz", df.columns[5]: "C_q(Equatorial I) in MHz"}, inplace = True)
df = df.set_index(list(df)[0])
df = df.iloc[:, 0:6:2]
df
import py4vasp
import numpy as np
cq = np.genfromtxt("./e07_efg_converge/kpoints_efg.dat", unpack=True, usecols=range(7))
py4vasp.plot(cq[0], cq[3], y2=True, xlabel="No. k-points along lattice parameter", y2label="Quadrupolar coupling constant (I) in MHz", label="Axial I (right y-axis)") \
+ py4vasp.plot(cq[0], cq[5], y2=True, y2label="Quadrupolar coupling constant (I) in MHz", label="Equatorial I (right y-axis)")\
+ py4vasp.plot(cq[0], cq[1], ylabel="Quadrupolar coupling constant (N) in MHz", label="N (left y-axis)")
It is immediately clear from the plot that the $\Gamma$-point is far from convergence. The N atom is converged to within 0.001 MHz by a 2x2x2 $\textbf{k}$-point mesh, while the I atoms require a much tighter 6x6x6 grid to have convergence within 1 MHz, which is more explicit in the tabulated data, i.e. it requires a denser mesh to describe the Brillouin zone. An important point to take from here is that different elements will have different tightnesses for convergence, 0.001 MHz for N, compared to below 1 MHz for I, i.e. the magnitude of the quadrupolar coupling constant is the key factor.
In Franssen et al., the experimental Cq measured using Nuclear Quadrupolar Resonance spectroscopy (NQR) are reported for 14N in Figure 2(B) [J. Phys. Chem. Lett. 8 (2017) 61]. From -100 to 75 °C, Cq decreases from 0.11 to 0.00 MHz, which is very different from out calculated values. The structure of methylamine of MAPbI3 is taken from the experimental structure and is calculated in the paper to be 0.45 MHz (cf. Figure S2), close to ours. So why is there a mismatch between experiment and calculation?
Sometimes having converged calculations is only part of the issue when comparing to experiment and further considerations must be taken into account. The methylamine can reorientate itself on the order of milliseconds, which is fast on the NMR timescale. If the average of four configurations is taken, then the calculated Cq follows the trend of the experiment, decreasing from 0.09 to 0.00 MHz between -100 and 50 °C. When these unanticipated dynamics are included, experiment and theory once again match.
7.4 Questions¶
- How dependent on the Brillouin zone is the quadrupolar coupling constant?
- Is the quadrupolar coupling constant more influenced by the plane-wave basis or the description of the Brillouin zone for MAPbI3?
By the end of this tutorial, you will be able to:
- Calculate the chemical shielding for LiH with and without two-center corrections
- Converge the shieldings with respect to plane-wave energy cutoff
- Compare the impact of two-center contributions to the shielding for Li and H
8.1 Task¶
Calculate chemical shieldings with and without two-center corrections for LiH in a 15x15x15 Å3cell converging with respect to plane-wave energy cutoff
Chemical shieldings are typically assumed to be one-center properties. Generally, this is valid and the contributions from the augmentation currents in neighboring PAW spheres is neglected. However, these two-center contributions can sometimes be substantial [J. Chem. Phys. 139 (2013) 014109]. Specifically, they are important when using harder PAW pseudopotentials (*_h) on atoms from the top rows of the periodic table (e.g. H, B, C, N, O, F), with short bonds.
In this exercise, you will calculate the chemical shieldings with and without two-center corrections for LiH.
8.2 Input¶
LiH molecule
1.0
15.00000000 0.00000000 0.00000000
0.00000000 15.00000000 0.00000000
0.00000000 0.00000000 15.00000000
Li H
1 1
Cartesian
7.50000000 7.50000000 6.70206900
7.50000000 7.50000000 8.29793200
INCAR.llraug
ENCUT = 400 # Plane-wave energy cutoff in eV
ISMEAR = 0; SIGMA = 0.01 # Defines the type of smearing; smearing width in eV
EDIFF = 1E-8 # Energy cutoff criterion for the SCF loop, in eV
PREC = Accurate # Sets the "precision" mode
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
LCHIMAG = .TRUE. # Turns on linear response for chemical shifts
LNMR_SYM_RED = .TRUE. # Consistent symmetry with star and k-space derivatives
NLSPLINE = .TRUE. # Differentiable projectors in reciprocal space
#LLRAUG = .TRUE. # Turns on two-center contributions to the chemical shifts
LASPH = .TRUE. # Non-spherical contributions to the gradient of the density in the PAW spheres
Gamma-point k-mesh
0
Monkhorst-Pack scheme
1 1 1
0 0 0
Pseudopotential of Li and H.
8.3 Calculation¶
Open a terminal and navigate to this example's directory:
cd $TUTORIALS/NMR/e08_*
In this directory, you can see the llraug.sh script. This script runs calculations at an energy cutoff of 400 eV, with and without the two-center contributions to the chemical shielding tensor (LLRAUG = .FALSE. and .TRUE., respectively. The chemical shieldings for the calculation with the following grep commands for the Li and H, respectively:
grep -A 5 "iso_shield" OUTCAR | tail -1 | awk '{print $6}'
grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}'
Results from 500 and 600 eV are then added to the end to show the convergence. The shieldings are written to llraug.dat after each calculation. Run llraug.sh with the following:
bash llraug.sh
Click to see the provided script to run the jobs and extract the chemical shieldings!
rm -f llraug.dat
for encut in 400
do
cp INCAR.llraug INCAR
sed -i "s/ENCUT = 400/ENCUT = $encut/g" INCAR
mpirun -np 4 vasp_std
shieldLi=$(grep -A 5 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
shieldH=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
echo "LLRAUG = .TRUE." >> INCAR
mpirun -np 4 vasp_std
shieldLi_AUG=$(grep -A 5 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
shieldH_AUG=$(grep -A 10 "iso_shield" OUTCAR | tail -1 | awk '{print $6}')
echo $encut $shieldLi $shieldLi_AUG $shieldH $shieldH_AUG >> llraug.dat
done
echo "500 88.6826 88.6863 24.6228 26.1131" >> llraug.dat
echo "600 88.8883 88.8918 24.6947 26.2643" >> llraug.dat
While the calculations are running, take a look at the structure of the molecule using py4vasp:
import py4vasp
struc = "./e08_two_center"
calc = py4vasp.Calculation.from_path(struc)
calc.structure.plot(selection="POSCAR")
#calc.structure.plot() # take a look at the structure
Congratulations on calculating the chemical shieldings! Let's plot them in py4vasp to see how much of a difference including the two-center contributions (2-c) makes:
import py4vasp
import numpy as np
import pandas as pd
shieldings = np.genfromtxt("./e08_two_center/llraug.dat", unpack=True, usecols=range(5))
py4vasp.plot(shieldings[0], shieldings[1], xlabel="Energy cutoff in eV", ylabel="Chemical shielding in ppm (Li)", label="Li") \
+ py4vasp.plot(shieldings[0], shieldings[2], ylabel="Chemical shielding in ppm (Li)", label="Li with 2-c")\
+ py4vasp.plot(shieldings[0], shieldings[3], y2=True, y2label="Chemical shielding in ppm (H)", label="H") \
+ py4vasp.plot(shieldings[0], shieldings[4], y2=True, y2label="Chemical shielding in ppm (H)", label="H with 2-c")
The difference in the chemical shieldings for Li with and without two-center contributions is small. However, for H it is much larger. Run the following script to tabulate the data and find out the exact values:
import numpy as np
import pandas as pd
df = pd.read_csv('./e08_two_center/llraug.dat', sep=' ', names=["ENCUT in eV", "σ(Li) in ppm", "σ(Li) with 2-c in ppm", "σ(H) in ppm", "σ(H) with 2-c in ppm"])
df = df.set_index(list(df)[0])
df['σ(Li) diff in ppm'] = df['σ(Li) with 2-c in ppm'] - df['σ(Li) in ppm']
df['σ(H) diff in ppm'] = df['σ(H) with 2-c in ppm'] - df['σ(H) in ppm']
df
The chemical shieldings for Li are largely converged by a cutoff of 600 eV, regardless of whether or not two-center contributions are included. This is also true for H. The main difference between Li and H is the magnitude of the contribution from the two-center part.
The two-center terms change little for Li, remaining at approximately ~0.004 ppm, almost negligible. For H, however, it is three orders of magnitude larger, ranging from ~1.1 to ~1.6 ppm going from a cutoff of 400 eV to 600 eV. This huge difference highlights how important it can be to consider the two-center contributions.
These two-center corrections are only necessary due to the use of PAW spheres. In a complete atom-centered basis, the shieldings are: 88.12 ppm (for Li) and 26.34 ppm (for H) [Phys. Chem. Chem. Phys. 18 (2016) 21145]. You can see how closely these match to your calculated values when the two-center approach is used.
8.4 Questions¶
- Which rows of the periodic table are two-center contributions to the current important?
- Are the two-center contributions more significant for Li or H?