Part 2: Molecules in VASP ¶
By the end of this tutorial, you will be able to:
- run a geometry relaxation using a conjugate-gradient algorithm to find the bond length of a dimer
- explain geometry relaxation considering Hellmann-Feyman forces and stress on the level of pseudocode
- explain the conjugate-gradient algorithm on the level of pseudocode
- set the step size for the conjugate-gradient algorithm and the number of ionic steps
4.1 Task¶
Relax the bond length of an oxygen dimer by means of geometry relaxation using a conjugate-gradient algorithm.
The full many-body Hamiltonian to describe two oxygen atoms bound to form a dimer comprises the electronic and the ionic degrees of freedom. Because ions are much heavier and, thus, slower, these degrees of freedom can be decoupled (Born-Oppenheimer approximation).
The term geometry-relaxation problem refers to the task of finding the ionic positions with zero forces and stress acting on them. On the level of pseudocode it reads
- given the ionic position, the electronic ground state can be computed within DFT,
- forces and stress acting on the ions can be computed as the expectation value of the gradient of the electronic Hamiltonian (Hellmann-Feyman theorem),
- ions are relaxed to their instantaneous ground state,
- repeat 1.-3. until convergence criterion is met.
For more information read the article about forces and the tags IBRION and ISIF.
4.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e04_O2-bond. Check them out!
O2 molecule in a box
1.0 ! universal scaling parameters
8.0 0.0 0.0 ! lattice vector a(1)
0.0 8.0 0.0 ! lattice vector a(2)
0.0 0.0 8.0 ! lattice vector a(3)
2 ! number of atoms
cart ! positions in cartesian coordinates
0 0 0 ! first atom
0 0 1.22 ! second atom
SYSTEM = O2 dimer in a box
ISMEAR = 0 ! Gaussian smearing
ISPIN = 2 ! spin-polarized calculation
NSW = 5 ! 5 ionic steps
IBRION = 2 ! use the conjugate-gradient algorithm
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Pseudopotential of O
The POSCAR file defines the position of two atoms in a large cubic box. They are 1.22 Å apart.
Concerning the INCAR file, in addition to the tags, that were set for the SDFT calculation discussed in 2 Spin-polarized oxygen atom, NSW and IBRION are set. Go ahead and check the corresponding articles on the VASP Wiki. You will learn about the conjugate-gradient algorithm and related tags such as POTIM.
Remind yourself why it is sufficient to consider one $\mathbf{k}$ point or revisit Section 1.2 to find the answer.
4.3 Calculation¶
Go ahead and run VASP for this example:
cd $TUTORIALS/molecules/e04_*
mpirun -np 2 vasp_gam
On the level of pseudocode, the conjugate-gradient algorithm reads
- make a steepest descent step starting from one ionic position:
- set the search direction equal to the direction of the largest gradient
- do line minimization until forces along the search direction for this ion become small
- update ionic position
- make a conjugate gradient step starting from the updated position of the same ion:
- find the direction of the largest gradient
- set the search direction equal to the direction of the conjugate gradient, that is obtained by orthogonalizing the largest gradient w.r.t. the search direction of the steepest descent step
- do line minimization until forces along the search direction for this ion become small
- update ionic position
- repeat 1.-2. until gradient becomes small or the maximum number of ionic steps is reached
Line minimization means that the minimum is determined only along the specific search direction, i.e., along a line. It is done using Brent's method with step size POTIM. Note that here force and gradient are synonyms. The maximum number of ionic steps is set by the NSW tag. The conjugate-gradient algorithm works well for few degrees of freedom ($\lesssim 4$) and if the initial guess is close to the ground state. Otherwise, you should employ a different method using the IBRION.
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_O2-bond")
my_calc.structure.print()
4.4 Questions¶
- How does the conjugate-gradient method work? When does it fail?
- What does POTIM and NSW control?
By the end of this tutorial, you will be able to:
- run VASP for multiple atomic species
- choosing the appropriate pseudopotential
- adjust the step size for the conjugate-gradient algorithm
- visualize structures using py4vasp
5.1 Task¶
Relax the bond length of a CO by means of geometry relaxation using a conjugate-gradient algorithm and visualize the result.
Pseudopotentials are introduced as an approximation to the real Coulomb potential created by the ionic sites which valence electrons observe. This approximation is necessary in practice in order to avoid the stark contrast between diverging Coulomb potentials at the nuclei and the slowly varying Coulomb potential in between ions. The complete potential landscape is constructed based on element specific pseudopotentials. Apart from this dependence on the chemical element, pseudopotentials also differ based on the number of electrons considered to be valence and based on how hard or soft they are, which are choices you have to make to some degree.
Read the article about available PAW potentials to learn about general recommendations. In which situation should you opt for hard pseudopotentials?
Click to see the answer!
If dimers with short bonds are present in the compound (O$_2$, CO, N$_2$, F$_2$, P$_2$, S$_2$, Cl$_2$), hard pseudopotentials are recommended.
5.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e05_CO-bond. Check the input files! You will find one file is missing and you will need to generate it as described below. Which file is missing?
Click to see the answer!
The POTCAR file.
CO dimer in a box
1.0 ! universal scaling parameters
8.0 0.0 0.0 ! lattice vector a(1)
0.0 8.0 0.0 ! lattice vector a(2)
0.0 0.0 8.0 ! lattice vector a(3)
1 1 ! number of atoms for each species
sel ! selective degrees of freedom are changed
cart ! positions in cartesian coordinates
0 0 0 F F T ! first atom
0 0 1.145 F F T ! second atom
SYSTEM = CO dimer in a box
ISMEAR = 0 ! Gaussian smearing
NSW = 5 ! 5 ionic steps
IBRION = 2 ! use the conjugate gradient algorithm
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Recall that, the POTCAR file contains the pseudopotential, the wavefunctions, the charge densities and other details about the atomic species. This is precomputed information that can be downloaded from the VASP portal, which is only accessible to license holders. However, in this example, you can find the relevant POTCAR files at $TUTORIALS/molecules/e05_CO-bond/pot.
Open the terminal and use cat as descibed in the VASP Wiki of the POTCAR file, in order to create the necessary POTCAR file for this example!
Click to get a hint!
In order to account for oxygen and carbon in this example, we need to concatenate the POTCAR of oxygen and carbon. Ask yourself: Which pseudopotentials are recommended according to the article available PAW potentials? Where does the final POTCAR file need to be located? Does the order of oxygen and carbon matter in this example?
Click to see the answer!
cd $TUTORIALS/molecules/e05_*
cat pot/C_h/POTCAR pot/O_h/POTCAR > POTCAR
Standard pseudopotentials are not a sufficient approximation for compounds featuring dimers with short bonds, so you should use hard pseudopotentials. Hard pseudopotentials have the disadvantage of needing a larger plane wave basis and, thus, are computationally demanding. For low accuracy calculations, such as in this tutorial, we may opt for the less demanding standard pseudopotentials instead.
5.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/molecules/e05_*
mpirun -np 2 vasp_gam
The stdout will show details about the conjugate-gradient algorithm.
Try running a fresh calculation with a smaller step size in the conjugate-gradient algorithm!
Click to get a hint!
In the directory of this example, we remove the WAVECAR file of the previous calculation to start a fresh calculation. Then, POTIM = 0.2 must be added to the INCAR file, which is smaller than the default value 0.5. Finally, run VASP.
Click to see how to complete the task!
cd $TUTORIALS/molecules/e05_*
rm WAVECAR
echo "POTIM = 0.2" >> INCAR
mpirun -np 2 vasp_gam
Vizualize the structure at each step using py4vasp!
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e05_CO-bond")
my_calc.structure[:].plot()
5.4 Questions¶
- In which order do you need to concatenate the POTCAR files?
- Which tag changes the step size in the conjugate-gradient algorithm?
- In which situation should you opt for hard pseudopotentials? Why?
By the end of this tutorial, you will be able to:
- compute vibration frequencies of molecules with VASP
- explain how the Hessian matrix and the phonon frequency are connected
6.1 Task¶
Compute the vibration frequency of a CO molecule using the method of finite differences.
In a real system, the bond length between two atoms is not fixed, but rather oscillates at frequencies ranging from approximately $10-100$ THz. These vibration modes of molecules are computed in the same way as phonons at the $\Gamma$ point. In principle it is necessary to compute the first and second derivative of the Born–Oppenheimer energy surface $E(\mathbf{R})$ that is defined by the ionic eigenvalue problem:
$$\tag{1}
\left( \sum_{\mu} \frac{\hbar^2}{2 M_{\mu}} \nabla_{\mathbf{R}_\mu}^2 + E(\mathbf{R}) \right) \zeta(\mathbf{R}) = \mathcal{E} \zeta(\mathbf{R}),
$$
where $\mathbf{R}_\mu$ is the coordinate of the $\mu$th ion with mass $M_{\mu}$ and $\zeta(\mathbf{R})$ and $\mathcal{E}$ are eigenfunctions and eigenvalues of the ionic system.
In practice, $E(\mathbf{R})$ is the groundstate energy of the electronic system, that parametrically depends on all ionic positions $\mathbf{R} = \left\{ \mathbf{R}_\mu \right\}$. Determining the vibration frequency $\omega$ then amounts to solving:
$$\tag{2}
\det \left| \frac{1}{\sqrt{ M_\mu M_\nu}} \frac{ \partial^2 E(\mathbf{R}) }{\partial R_{\mu i} \partial R_{\nu j} } - \omega^2 \right| = 0,
$$
where $\frac{ \partial^2 E(\mathbf{R}) }{\partial R_{\mu i} \partial R_{\nu j} }$ is the Hessian with respect to the ionic site $\mu, \nu$ and direction $i,j=x, y, z$.
Check out this article about phonons on the VASP Wiki for more information.
6.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e06_CO-vibration.
For the INCAR, check out the meaning of IBRION = 5, NFREE and POTIM in this context. You will learn how the Hessian matrix is computed and its connection to phonons.
CO dimer in a box
1.0 ! universal scaling parameters
8.0 0.0 0.0 ! lattice vector a(1)
0.0 8.0 0.0 ! lattice vector a(2)
0.0 0.0 8.0 ! lattice vector a(3)
1 1 ! number of atoms for each species
sel ! selective degrees of freedom are changed
cart ! positions in cartesian coordinates
0 0 0 F F T ! first atom
0 0 1.143 F F T ! second atom
SYSTEM = CO dimer in a box
ISMEAR = 0 ! Gaussian smearing
IBRION = 5 ! use the conjugate gradient algorithm
NFREE = 2 ! central differences
POTIM = 0.02 ! 0.02 A stepwidth
NSW = 1 ! ionic steps > 0
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Hard pseudopotentials of C and O.
6.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/molecules/e06_*
mpirun -np 2 vasp_gam
Try to find the normal mode of the CO molecule! To this end, open the OUTCAR file, at the end you will find the dynamical information below
Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------
Click to see the answer!
According to Mina-Camilde et al., J Chem Educ (1996) the experimentally observed stretching vibration of CO has a frequency of 2143 cm$^{-1}$, i.e., 64.25 THz, and, thus, the computed value 63.887522 THz is quite close to reality.
6.4 Questions¶
- How is the Hessian matrix defined?
- How is the Hessian matrix computed with finite difference?
- What does the POTIM tag set?
- Which tags must be set in order to compute the vibrational frequencies?
- In how many units is the vibrational frequency written to the output?
By the end of this tutorial, you will be able to:
- plot the density of states (DOS) using py4vasp
- explain the difference between the partial DOS, local DOS and the total DOS
7.1 Task¶
Calculate the partial DOS of a CO molecule.
The DOS describes the number of electronic states in some energy interval: $N = \int_{-\infty}^\infty f(\epsilon) D(\epsilon)\, \text{d}\epsilon$, where $N$ is the total number of states, $f(\epsilon)$ is the Fermi function, $D(\epsilon)$ is the DOS and $\epsilon$ is the energy variable. For interacting electrons at finite temperature this corresponds to the spectral function, that can be measured in experiments such as photo-electron spectroscopy (PES).
Generally, one can distinguish the partial DOS, local DOS and the total DOS. The partial DOS is the projection onto a specific orbital character. For instance, this could be the partial DOS due to oxygen $2p$ orbitals. The total DOS refers to the entire system. In this example, this is the entire CO molecule. Finally, in a crystal rather than an isolated molecule the electronic wavefunction depends on the wavevector $\mathbf{k}$. In that case, also the DOS is $\mathbf{k}$-dependent and the local DOS is obtained by integrating over $\text{d}\mathbf{k}$.
7.2 Input¶
The input files to run this example are prepared at $TUTORIALS/molecules/e07_CO-partial-dos.
Check them out!
CO dimer in a box
1.0 ! universal scaling parameters
8.0 0.0 0.0 ! lattice vector a(1)
0.0 8.0 0.0 ! lattice vector a(2)
0.0 0.0 8.0 ! lattice vector a(3)
1 1 ! number of atoms for each species
sel ! selective degrees of freedom are changed
cart ! positions in cartesian coordinates
0 0 0 F F T ! first atom
0 0 1.143 F F T ! second atom
SYSTEM = CO dimer in a box
ISMEAR = 0 ! Gaussian smearing
LORBIT = 11
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
Hard pseudopotentials of C and O.
In the INCAR file, LORBIT = 11 is set. Get familiar with the output generated due to this setting on the VASP Wiki!
Notice that setting selective degrees of freedom in the POSCAR file has no effect on the calculation if you only perform a DFT calculation.
7.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/molecules/e07_*
mpirun -np 2 vasp_gam
Use py4vasp to plot the DOS as follows. Select the partial DOS by parsing a string to the plot function.
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e07_CO-partial-dos")
my_calc.dos.plot()
# uncomment the line below to plot the partial DOS
# my_calc.dos.plot("p")
7.4 Questions¶
- What output is generated by seeting LORBIT = 11?
- How can you plot the partial DOS?