Part 3: Water ¶
By the end of this tutorial, you will be able to:
- explain the residual minimization method with direct inversion in the iterative subspace (RMM-DIIS) on the level of pseudocode
- judge whether to use the RMM-DIIS or conjugate-gradient (CG) algorithm
- write a POSCAR file from scratch
- use the scaling parameter in the POSCAR file
- perform a geometry relaxation with two degrees of freedom
- set the convergence criterion for the ionic relaxiation
8.1 Task¶
Relax the bond length of a water molecule by means of geometry relaxation using a RMM-DIIS.
RMM-DIIS, also known as generalized minimal residual method (GMRES), is a quasi-Newton algorithm, like the CG algorithm, which relaxes the ionic position to the instantaneous ground state. However there are fundamental differences in the minimization method.
Essentially in DIIS, the minimization problem is iteratively solved for a few iterations. Then, the final result is obtained in that subspace using a direct method. In particular, RMM is the direct method here. Mind that the quantity minimized is the force acting on the ions and not the total energy. For a differentiable potential $V(\mathbf{r}; \mathbf{R})$, that acts on the electrons and depends parametrically on all ionic positions $\mathbf{R}=\left\{ \mathbf{R}_\mu \right\}$, the force and Hessian are obtained using the Hellmann–Feynman theorem:
$$\tag{1a}
F_{\mu i}= \frac{\partial E}{\partial R_{\mu i}} = \int \frac{\partial V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i}} n(\mathbf{r}; \mathbf{R}) \text{d} \mathbf{r},
$$
$$\tag{1b}
H_{\mu i \nu j}=\frac{\partial^2 E}{\partial R_{\mu i} \partial R_{\nu j} } = \int \frac{\partial^2 V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i} \partial R_{\nu j}} n(\mathbf{r}; \mathbf{R}) \text{d} \mathbf{r} + \int \frac{\partial V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i}} \frac{n(\mathbf{r}; \mathbf{R})}{\partial R_{\nu j}} \text{d} \mathbf{r}.
$$
Here, $\mu, \nu$ labels the ionic site, $i,j=x,y,z$ and $n(\mathbf{r}; \mathbf{R})$ is the electronic charge density. RMM-DIIS then minimizes the norm of the residual vector, i.e., $|\mathbf{F}|^2$, where
$$\tag{2}
F_{\mu i}= \sum_{\nu j} H_{\mu i \nu j} \partial R_{\nu j}.
$$
Note that RMM-DIIS breaks the minimization up into subspaces according to the spectrum of the Hessian. It makes this algorithm very fast, but also sensitive to initialization problems and sporadically other difficulties, that need special attention by the user. It is a good choice for systems with ionic positions close to their minimum and $\gtrsim 3$ degrees of freedom. However, for $\gtrsim 20$ degrees of freedom and very broad vibration spectra, other methods such as molecular dynamics are recommended.
8.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e08_H2O-relaxation.
Check them out!
Set up your own POSCAR file for H$_2$O! Place the oxygen atom at the origin, use an approximate H - O bond length of $0.96$ Å and $105^{\circ}$ as an angle between the H - O bonds. Consider which box size might be sufficiently large to obtain the result of an isolated H$_2$O molecule. Finally, you need to allow two degrees of freedom for each H atom using the selective dynamics mode.
Click to get a hint!
You can use the template provided in the directory of this example.
cd $TUTORIALS/molecules/e08_*
cp POSCAR_template POSCAR
Click to see the answer!
H2O
0.52918 ! scaling parameter
12 0 0
0 12 0
0 0 12
1 2
select
cart
0.00 0.00 0.00 F F F ! O
1.10 -1.43 0.00 T T F ! H
1.10 1.43 0.00 T T F ! H
Here, we used the scaling parameter in the second line of the POSCAR file. It is a multiplicative factor for the lattice matrix and atomic coordinates.
The KPOINTS and INCAR files are provided in this example's directory. Concerning the INCAR file, check out the meaning of IBRION = 1, NFREE, and EDIFFG! Do you recall what the remaining tags in the INCAR file used here control? ENCUT is increased by $30\%$ from its default value. This is recommended when performing ionic relaxation, and generally when computing the stress tensor and pressure.
Click to get a hint!
Mind the order in which you concatenate the POTCAR files.
Click to see the answer!
In this example's directory, enter the following into the terminal:
cd $TUTORIALS/molecules/e08_*
cat pot/O/POTCAR pot/H/POTCAR > POTCAR
8.3 Calculation¶
What H - O bond length and angle between bonds do you obtain?
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e08_H2O-relaxation")
my_calc.structure.print()
Click to see the answer!
H2O
0.529180000000000
12.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
O H
1 2
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.0937992998266598 0.8791568206830546 0.0000000000000000 T T F
0.0937992998266598 0.1208431793169455 0.0000000000000000 T T F
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
$$ \cos(\alpha) = \frac{\mathbf{b}_1 \cdot \mathbf{b}_2}{|\mathbf{b}_1| |\mathbf{b}_2|} $$
Click to see H - O bond length and angle!
H - O bond length: (0.0937992998266598^2+ 0.1208431793169455^2)^0.5120.52918 = 0.9714163904 Å
angle:
arcos((0.09379929982665981^2-0.12084317931694552^2)/(0.09379929982665981^2+0.12084317931694552^2)) = 104.3622835581412°
Alternatively, you can use py4vasp and extract the information directly from the structure plot. First, plot a 2x2x2 supercell. Then, right click the first atom and right click the second atom twice to get the distance. For the angle, right click atom 1, then atom 2 and right click atom 3 twice. Check that the result agrees with the discussion above!
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e08_H2O-relaxation")
my_calc.structure.plot(2)
8.4 Questions¶
- What does EDIFFG set?
- Describe the RMM-DIIS algorithm on the level of pseudocode?
- How can the ionic degrees of freedom be set in the POSCAR file?
- What is the effect of the scaling parameter given in the second line of the POSCAR file?
By the end of this tutorial, you will be able to:
- set the energy cutoff for the plane wave basis of the pseudo-orbitals
- set the precision of a calculation using PREC
9.1 Task¶
Determine the precision for the DFT calculation of an isolated water molecule.
In the framework of the projector-augmented wave (PAW) method, the ansatz for the Kohn-Sham (KS) orbitals,
$$\tag{1}
| \psi_{i\mathbf{k}\sigma} \rangle = |\tilde{\psi}_{i\mathbf{k}\sigma} \rangle + \sum_\eta \left( |\phi_{\eta} \rangle - |\tilde{\phi}_{\eta} \rangle \right) \langle \tilde{p}_\eta |\tilde{\psi}_{i\mathbf{k}\sigma} \rangle,
$$
comprise a delocalized and an atom-centered part.
Let us focus on the first term, $|\tilde{\psi}_{i\mathbf{k}\sigma} \rangle$, the delocalized pseudo-orbitals, that are expressed in terms of plane waves. In real space this reads
$$\tag{2}
\langle \mathbf{r} | \tilde{\psi}_{i\mathbf{k}\sigma} \rangle = \frac{1}{\sqrt{\Omega_r}} \sum_{\mathbf{G}} C_{i\mathbf{G}\mathbf{k}}^\sigma(\mathbf{r}) \text{e}^{\text{i} \left(\mathbf{G}+\mathbf{k}\right) \cdot \mathbf{r}},
$$
where $\mathbf{r}$ is the real space coordinate, $i$ is the band index enumerating the KS equations, $\sigma$ is the spin index, $\Omega_r$ is the real space volume, $\mathbf{G}$ is the wave vector of the plane wave, $\mathbf{k}$ is the wave vector defined in the KPOINTS file and $C_{i\mathbf{G}\mathbf{k}}^\sigma(\mathbf{r})$ are the variational expansion coefficients.
In principle, Equation $(1)$ is a complete basis if infinitely many plane waves are taken into account in Equation $(2)$. But in practice, strongly oscillating components in real space are not expected to significantly contribute to the computed physical quantities and thus we can save a lot of computational effort by neglecting these components. This is done by introducing an energy cutoff:
$$\tag{3}
\frac{1}{2} | \mathbf{G} + \mathbf{k} |^2 < E_{\text{cutoff}}.
$$
The choice of the energy cutoff is mainly influenced by the choice of the pseudopotential, but it also depends on the specific system and the property of interest. To set the appropriate energy cutoff for your calculation, you need to perform a convergence study.
For more information read about the PAW method and the energy cutoff on the VASP Wiki.
9.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e09_H2O-energy-cutoff.
Check out the subdirectories! Each contains the same setup for a DFT calculation, except for the energy cutoff.
Concerning the INCAR file, check out the meaning of the PREC and the ENCUT tag!
H2O
0.529180000000000
12.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
O H
1 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0937992998266598 -0.1208431793169455 0.0000000000000000
0.0937992998266598 0.1208431793169455 0.0000000000000000
SYSTEM = H2O
PREC = Normal ! standard precision
ISMEAR = 0
SIGMA = 0.1
ENCUT = 400 ! up to 820
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Pseudopotentials of O and H.
9.3 Calculation¶
Perform multiple VASP calculations with different energy cutoff in each subdirectory named ENCUTXXX!
Click to see the answer!
In each subdirectory named ENCUTXXX, run vasp:
mpirun -np 2 vasp_gam
You can run all with these lines
cd $TUTORIALS/molecules/e09_*
for d in $(ls | grep ENCUT); do cd $d ; mpirun -np 2 vasp_std ; cd .. ; done
Plot the total energy as a function of the energy cutoff. What can you say about the precision of your calculation? Which cutoff energy should you choose?
Click to get a hint!
You can use the bash and gnuplot scripts provided in this example's directory by entering the following into the terminal:
bash convergence_study.sh
gnuplot convergence_study.gp
This generates convergence_study.dat and convergence_study.png. Open the PNG from the file browser and ask yourself again: What can you say about the precision of your calculation? Which cutoff energy should you choose?
Click to see the answer!
The total energy is converging to the value it would have, if the plane wave basis used in the PAW method were complete, as opposed to being cut off at the cutoff energy specified by ENCUT. The precision of the quantities you are computing is fundamentally limited by the completeness of your basis, but there is no direct way to infer what cutoff you need to compute a specific quantity. So, best practice would be, to compute the specific quantity at different values of ENCUT and make a convergence study for that quantity in analogy to the convergence study you performed in this example for the total energy.
9.4 Questions¶
- What are pseudo-orbitals?
- Is the basis used in the PAW method complete? Why?
- What is controlled by the PREC tag?
By the end of this tutorial, you will be able to:
- explain the difference between the method of finite differences and density-functional-perturbation theory (DFPT)
- set the number of KS orbitals/electronic bands
10.1 Task¶
Compute the vibration frequency of a H$_2$O molecule by finite differences and DFPT.
The method of finite differences and DFPT are both capable to compute the second derivatives, Hessian matrix and phonon frequencies. Recall that for a differentiable potential $V(\mathbf{r}; \mathbf{R})$, that acts on the electrons and depends parametrically on all ionic positions $\mathbf{R}=\left\{ \mathbf{R}_\mu \right\}$, the force and Hessian read as follows according to the Hellmann–Feynman theorem:
$$\tag{1a}
\frac{\partial E}{\partial R_{\mu i}} = \int \frac{\partial V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i}} n(\mathbf{r}; \mathbf{R}) \text{d} \mathbf{r},
$$
$$\tag{1b}
\frac{\partial^2 E}{\partial R_{\mu i} \partial R_{\nu j} } = \int \frac{\partial^2 V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i} \partial R_{\nu j}} n(\mathbf{r}; \mathbf{R}) \text{d} \mathbf{r} + \int \frac{\partial V(\mathbf{r}; \mathbf{R})}{\partial R_{\mu i}} \frac{n(\mathbf{r}; \mathbf{R})}{\partial R_{\nu j}} \text{d} \mathbf{r}.
$$
Here, $\mu, \nu$ is the ionic site, $i,j=x,y,z$ and $n(\mathbf{r}; \mathbf{R})$ is the electronic charge density. The linear change in the charge density yields
$$\tag{2}
\Delta^\mathbf{R} n(\mathbf{r}; \mathbf{R})= 2 \text{Re} \sum_{i,\sigma}^{occ} \psi_{i\sigma}^*(\mathbf{r}; \mathbf{R})\Delta^\mathbf{R} \psi_{i\sigma}(\mathbf{r}; \mathbf{R}),
$$
where $\Delta^\mathbf{R}$ is the finite-difference operator that describes the displacement of ions and is defined as
$$\tag{3}
\Delta^\mathbf{R} f = \sum_\mu \sum_{i=x,y,z} \frac{\partial f}{\partial R_{\mu i}} \Delta R_{\mu i}.
$$
Let us review DFPT in a nutshell. The variation of the KS orbitals $\Delta^\mathbf{R} \psi_{i\sigma}(\mathbf{r}; \mathbf{R})$ that appears in Equation (2) is obtained by standard first-order perturbation theory and yields the so-called Sternheimer equation:
$$\tag{4}
\left( H - \epsilon_{i\sigma} \right) | \Delta^{\mathbf{R}} \psi_{i\sigma} \rangle = - \left( \Delta^\mathbf{R} V - \Delta^{\mathbf{R}} \epsilon_{i\sigma} \right) | \psi_{i\sigma} \rangle,
$$
where
$$\tag{5}
H = - \frac{\hbar^2}{2m} \nabla^2_{\mathbf{r}} + V(\mathbf{r}; \mathbf{R})
$$
is the unperturbed KS Hamiltonian,
$$\tag{6}
\Delta^\mathbf{R} V(\mathbf{r}; \mathbf{R}) = \Delta^\mathbf{R} V_{\text{ext}}(\mathbf{r}; \mathbf{R}) + e^2 \int \frac{ \Delta^\mathbf{R} n(\mathbf{r}'; \mathbf{R})}{|\mathbf{r} - \mathbf{r}'|} \text{d}\mathbf{r}' + \left. \frac{\text{d} v_{xc}(n)}{\text{d} n } \right|_{n=n(\mathbf{r}; \mathbf{R})} \Delta^\mathbf{R} n(\mathbf{r}; \mathbf{R}),
$$
and $\Delta^{\mathbf{R}} \epsilon_{i\sigma} = \langle \psi_{i\sigma} | \Delta^\mathbf{R} V | \psi_{i\sigma} \rangle$. Note that given $\Delta^\mathbf{R} n$ in Equation (2), $\Delta^\mathbf{R} V$ in Equation (6) is defined. This allows to solve the Sternheimer equation in Equation (4) to obtain $\Delta^{\mathbf{R}} \psi_{i\sigma}$. Then, $\Delta^\mathbf{R} n$ in Equation (2) can be updated. After self-consistency is reached the forces and Hessian can be computed as given in Equation (1a) and (1b), respectively.
The advantage of DFPT over the method of finite differences is that only occupied bands are needed. Mind that in principal DFPT has another great advantage, namely that the response to perturbations of different wavelengths are decoupled, which allows the calculation of phonon frequencies at arbitrary wave vectors $\mathbf{q}$ and avoiding the use of supercells. However, this is not implemented in VASP and irrelevant to the vibration frequencies of isolated molecules.
10.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e10_H2O-vibration.
Check out the subdirectories! IBRION6 contains the finite difference calculation and IBRION8 the calculation using perturbation theory.
Concerning the INCAR file, check out the meaning of PREC = Accurate and set ENCUT based on the convergence study in the previous example! Check the meaning of the remaining tags, such as NBANDS.
H2O
0.529180000000000
12.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
O H
1 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0937992998266598 -0.1208431793169455 0.0000000000000000
0.0937992998266598 0.1208431793169455 0.0000000000000000
IBRION6/INCAR
SYSTEM = H2O vibration
PREC = Accurate ! recommended for computing forces
ENCUT = 520
ISMEAR = 0
SIGMA = 0.1
IBRION = 6 ! finite differences with symmetry
NFREE = 2 ! central differences (default)
POTIM = 0.015 ! step size (default)
EDIFF = 1E-8 ! convergence criterion
NSW = 1 ! ionic steps > 0
NBANDS = 6
IBRION8/INCAR
SYSTEM = H2O vibration
PREC = Accurate ! recommended for computing forces
ENCUT = 520
ISMEAR = 0 ! Gaussian smearing
IBRION = 8 ! perturbation theory
NFREE = 2 ! central differences (default)
POTIM = 0.015 ! step size (default)
EDIFF = 1E-8 ! convergence criterion
NSW = 1 ! ionic steps > 0
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Pseudopotentials of O and H.
10.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/molecules/e10_*/IBRION6
mpirun -np 2 vasp_gam
cd $TUTORIALS/molecules/e10_*/IBRION8
mpirun -np 2 vasp_gam
How many zero frequency modes should be observed and why? Use the linear response calculation (IBRION=8 and EDIFF=1E-8) as a reference result. For finite differences, are the results sensitive to the step width POTIM?
Click to see vibration modes computed with the method of finite differences!
IBRION6, POTIM = 0.015 and ENCUT = 800 yields:
Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------
1 f = 115.516136 THz 725.809286 2PiTHz 3853.203528 cm-1 477.736135 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 -0.269007 0.000000
0.595641 5.582786 0.000000 -0.417541 0.538031 0.000000
0.595641 0.767374 0.000000 0.417541 0.538031 -0.000000
2 f = 112.251954 THz 705.299826 2PiTHz 3744.322138 cm-1 464.236570 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.192188 0.000000 -0.000000
0.595641 5.582786 0.000000 -0.384378 0.577742 -0.000000
0.595641 0.767374 0.000000 -0.384378 -0.577742 -0.000000
3 f = 47.648480 THz 299.384231 2PiTHz 1589.382220 cm-1 197.058192 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.272326 0.000000 0.000000
0.595641 5.582786 0.000000 -0.544706 -0.407694 0.000000
0.595641 0.767374 0.000000 -0.544706 0.407694 0.000000
4 f = 0.013547 THz 0.085117 2PiTHz 0.451874 cm-1 0.056025 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.000000 -0.000000 -0.940920
0.595641 5.582786 0.000000 -0.000000 -0.000000 -0.239446
0.595641 0.767374 0.000000 -0.000000 -0.000000 -0.239446
5 f = 0.008845 THz 0.055572 2PiTHz 0.295023 cm-1 0.036578 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.000000 -0.942322 0.000000
0.595641 5.582786 0.000000 -0.001412 -0.236669 0.000000
0.595641 0.767374 0.000000 0.001412 -0.236669 0.000000
6 f/i= 0.001194 THz 0.007504 2PiTHz 0.039839 cm-1 0.004939 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.942816 -0.000000 -0.000000
0.595641 5.582786 0.000000 0.235688 -0.000010 -0.000000
0.595641 0.767374 0.000000 0.235688 0.000010 -0.000000
7 f/i= 2.409664 THz 15.140367 2PiTHz 80.377747 cm-1 9.965566 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.000000 -0.000000 0.338628
0.595641 5.582786 0.000000 -0.000000 -0.000000 -0.665331
0.595641 0.767374 0.000000 -0.000000 -0.000000 -0.665331
8 f/i= 2.838950 THz 17.837649 2PiTHz 94.697181 cm-1 11.740949 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.000000 0.000000 0.000000
0.595641 5.582786 0.000000 -0.000000 -0.000000 -0.707107
0.595641 0.767374 0.000000 0.000000 -0.000000 0.707107
9 f/i= 3.602905 THz 22.637719 2PiTHz 120.179968 cm-1 14.900410 meV
X Y Z dx dy dz
0.000000 0.000000 0.000000 0.000000 0.199158 0.000000
0.595641 5.582786 0.000000 -0.570664 -0.393078 0.000000
0.595641 0.767374 0.000000 0.570664 -0.393078 -0.000000
Click to see the answer!
In this specific case, the drift in the forces is too large to obtain the zero frequency modes "exactly", and it is simplest to increase the cutoff ENCUT to 800 eV. The important and physically meaningful frequencies are, however, insensitive to the choice of the cutoff.
10.4 Questions¶
- What is the advantage of density-functional-perturbation theory compared to the method of finite differences? Are there disadvatages?
- What is the default of NBANDS?
By the end of this tutorial, you will be able to:
- switch off the use of symmetry
- perform a standard ab-initio molecular dynamics (MD) calculation
- explain the basic concept of Nose-Hoover thermostat
11.1 Task¶
Perform a standard ab-initio MD calculation for water and plot the pair-correlation function.
In MD calculations the motion of atoms (or molecules) at a specific temperature is simulated by means of Newton's second law. In other words, each iteration simulated a time step, where atoms are treated as classical particles exhibiting forces. When these forces are computed quantum mechanically using ab-initio methods, one speaks of ab-initio MD.
Ab-initio MD can be used to solve the ionic relaxation problem, when slowly reducing the temperature. This so-called simulated annealing approach can for instance be realized by introducing friction into the system as is done in damped MD. It is particularly well-suited to treat systems with $\gtrsim 20$ degrees of freedom and very broad vibrational spectra. Thus, the ionic relaxation problem of this water molecule is better treated using RMM-DIIS, as you did in 8 Bond length in H$_2$O.
The advantage of ab-initio MD over the RMM-DIIS and CG algorithms is that the MD trajectories, i.e., the ionic position after each time step, are physically meaningful. It enables the computation of the pair-correlation function, which is the probability of finding the a particle at a given distance from the center of another particle.
In order to learn more about MD algorithms in VASP and how the effect of temperature is included by means of the Nosé-Hoover thermostat in this calculation, read the linked VASP Wiki articles.
11.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e11_H2O-MD.
H2O
0.529180000000000
12.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
O H
1 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0937992998266598 -0.1208431793169455 0.0000000000000000
0.0937992998266598 0.1208431793169455 0.0000000000000000
PREC = Normal ! standard precision
ENCUT = 520
ISMEAR = 0
SIGMA = 0.1
IBRION = 0 ! standard molecular dynamics (MD)
NSW = 500 ! 500 steps
POTIM = 0.5 ! timestep 0.5 fs
ISYM = 0 ! no imposed symmetry for MD
SMASS = -3
TEBEG = 2000 ! temperature at beginning
TEEND = 2000 ! temperature at end
NBANDS = 8
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Pseudopotentials of O and H.
11.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/molecules/e11_*
mpirun -np 2 vasp_gam
Plot the total energy of the structure at each time step using py4vasp:
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e11_H2O-MD")
my_calc.energy[:].plot()
The amplitude of the oscillation in the total energy is due to the temperature of the system. That is, the position of the ions are slightly moved in accordance with the Nosé-Hoover thermostat depending on the temperature and the time step. After the electronic relaxation, this results in a new total energy. As the temperature at the beginning and at the end of this calculation is the same, the oscillating motion and thus the oscillation in the total energy follows the same distribution throughout this calculation.
Plot the motion of the ionic positions using py4vasp:
# If the last code box you have executed is the code box just above,
# there is no need to run the following two commands again!
# import py4vasp
# my_calc = py4vasp.Calculation.from_path( "./e11_H2O-MD")
my_calc.structure[:].plot()
A less animated, but quantitative way to capture the probability of finding the a particle at a given distance from the center of another particle is the pair-correlation function. The data of the pair-correlation function is written to the PCDAT file. Use the provided scripts to plot the result by entering the following at this example's directory into the terminal:
cd $TUTORIALS/molecules/e11_*
bash pair_correlation.sh
gnuplot pair_correlation.gp
Then, open the generated PNG file named pair_correlation.png from the file browser. What is the interpretation of the pair-correlation function at large distances? Why does the long range behavior start around 6 Å?
Click to see the answer!
The pair-correlation function at large distances corresponds to the density of the system. In this calculation, the box size is $0.52918 \cdot 12 = 6.35016$ Å, so that could be a rough estimate where the long range behavior starts. The reason there is weight even below 6 Å, is that the distance between an H atom and a O atom of a neighboring H$_2$O molecule is smaller than 6 Å.
What is the interpretation of the pair-correlation function at short distances?
Click to see the answer!
In short range, the pair-correlation function is shaped by the interaction between pairs within one H$_2$O molecule.
Why is there a peak around 1 Å and a double peak at approximately 1.53 Å? Consider that the animation shows a symmetric stretch of the molecule as a dominant vibration mode. Use the vibration modes computed in the Example 10 to estimate the splitting of the peak around 1.53 Å!
Click to see the answer!
In the ground state, the H - O bond is close to $0.97$ Å and the distance between O atoms is $1.53$ Å. The vibration of the molecule broadens the peaks and the vibration modes introduces a peak splitting. From Example 10, we know that there is a $3853.203528$ cm$^{-1}$. Converting to Ångström this yields: $1/3.853203528 = 0.2595243134013865$ Å, which agrees with the distance between the two peaks!
11.4 Questions¶
- What is a pair-correlation function?
- Can molecular dynamics be used to solve the ionic relaxation problem? If yes, when do you opt for MD? If no, why not?
- What does the tag SMASS control?
- What is the unit of TEBEG?