Part 2: More silicon ¶
By the end of this tutorial, you will be able to:
- perform a volume relaxation manually
- obtain the density of states (DOS) and band structure
of a material with minimal guidance in the framework of density-functional theory (DFT).
4.1 Task¶
Perform a volume relaxation and compute the DOS and band structure along a high symmetry path of your choice for cubic-diamond silicon!
Recall what you did in Part 1 and apply the same procedure to another stable structure of silicon: the cubic diamond. It is possible to do this by choosing your input close to the expected structure. The result will serve as a benchmark for the discussion of the following example!
Here, we take advantage of the fact that the algorithm converges to a nearby local minimum in the energy landscape, but let us seize the opportunity to express a word of caution. Sometimes you might be tempted to claim your ab-initio calculation corresponds to the system's ground state, but in fact, it may correspond to a metastable state, even if the underlying physics is well captured within the framework of DFT. This might happen because the initial parameters are too far from the ground state and your calculation converged to a local minimum. The upside of this issue is that different crystal structures, magnetic configurations, etc. can be compared to each other. In order to arrive at the correct interpretation of your calculation, it is good practice to actually compare your result to nature, or in other words to the corresponding experiment, if it is available.
4.2 Input¶
The input files to run this example are prepared at $TUTORIALS/bulk/e04_cd-Si. Check them out and add any files you need to complete this example!
What are direct coordinates? If you are not sure, you may want to read the POSCAR article on the VASP Wiki!
The above input files are a good starting point to perform the volume relaxation. Recall what settings are possible for ISMEAR, ISTART, and ICHARG depending on the purpose of your calculation. What does the EDIFF tag control?
4.3 Calculation¶
Let's do this in steps! Change to this example's directory:
cd $TUTORIALS/bulk/e04_*
Step 1.)¶
Determine the optimal lattice constant up to 2 decimal places starting from a guess interval $[5.3,5.5]$ Ångström.
Click to get a hint!
You can use the provided script loop_lattice_constant.sh
for a in 5.3 5.4 5.5
do
cat > POSCAR << EOF
cubic diamond Si
$a
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.125
EOF
mpirun -np 2 vasp_std
en=$(awk '/F=/ {print $0}' OSZICAR)
echo $a $en >> loop_lattice_constant.dat
done
Then, plot the result with loop_lattice_constant.gp
set term png
set output "loop_lattice_constant.png"
set title "lattice constant of cd Si"
set xlabel "lattice constant [Angstrom]"
set ylabel "total free energy"
plot "loop_lattice_constant.dat" using 1:4 w lp
by entering the following into the terminal
gnuplot loop_lattice_constant.gp
and open loop_lattice_constant.png to check the result. Change the looped lattice constants until you reach the desired precision.
What is the optimal lattice constant? Back up the charge density of the self-consistent solution at the optimal lattice constant.
Click to see the answer!
Compute the energy of three lattice constants, then select the lowest two. The optimal minimum must be in that new interval. Make your next guess at the center of the two and repeat. The final lattice constant is $5.46(8)$ Å. Do not forget to make a final run with the optimal lattice constant in order to write the CHGCAR file and make a copy as a backup.
cubic diamond Si
5.468
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.125
mpirun -np 2 vasp_std
cp CHGCAR CHGCAR.bk
Click to see the answer!
cubic diamond Si
5.468
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.125
System = cd Si DOS
ENCUT = 400
ISMEAR = -5
LORBIT = 11
K-Points
0
Monkhorst Pack
15 15 15
0 0 0
The DOS is not as sensitive to the cutoff energy as the volume relaxation, so you can use a smaller value for ENCUT here. ISMEAR = -5 is the tetrahedron method with Blöchl correction and LORBIT = 11 allows DOSCAR to be written.
Plot the DOS using py4vasp! Could you have used ISMEAR = -5 for the volume relaxation?
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e04_cd-Si" )
mycalc.dos.plot()
Click to see the answer!
Yes, because cd silicon is not metallic and it is no problem that the energy is not variational with respect to partial occupancies, because there are no partially filled bands.
Click to see the answer!
System = cd Si band structure
ENCUT = 540
ISMEAR = 0
SIGMA = 0.1
ICHARG = 11
kpoints for band structure L-G-X-U K-G
10
line
reciprocal
0.50000 0.50000 0.50000 ! L
0.00000 0.00000 0.00000 ! G
0.00000 0.00000 0.00000 ! G
0.50000 0.00000 0.50000 ! X
0.50000 0.00000 0.50000 ! X
0.62500 0.25000 0.62500 ! U
0.37500 0.37500 0.75000 ! K
0.00000 0.00000 0.00000 ! G
Plot the band structure using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e04_cd-Si" )
mycalc.band.plot()
4.4 Questions¶
- What does LORBIT = 11 control?
- Is silicon more stable in its face-centered-cubic structure looked at in Part 1, or in its cubic-diamond structure computed here? Compare their total energy!
By the end of this tutorial, you will be able to:
- avoid Pulay stresses and forces
- relax the unit cell volume, cell shape, and ionic positions with a single VASP run
- set the ISIF tag
5.1 Task¶
Relax the volume, cell shape, and ionic positions of cd Si using the conjugate-gradient algorithm.
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.
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{k}\sigma}(\mathbf{G}) \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{k}\sigma}(\mathbf{G})$ 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}}.
$$
In the context of volume relaxation, the incompleteness of the plane-wave expansion in Equation $(2)$ causes a bias towards smaller unit cells, which is called Pulay stress. Intuitively, this is so because for a smaller volume the energy cutoff is effectively larger and the total energy decreases for increasing energy cutoff. The connection between volume and energy cutoff is given by Equation $(3)$, where the reciprocal vector $\mathbf{G}$ is defined with respect to a specific cell in real space. In other words, two independent VASP calculations with different volumes have a different basis, even when the energy cutoff is the same.
Pulay stress is unphysical stress and can be avoided by restarting your calculation at every volume as done in Example 4. Conveniently, it is also possible to optimize the volume in a single run. To avoid Pulay stress in that case, you need to set ENCUT appropriately high and restart your calculation when the volume changed significantly.
Check the linked article to learn about how to avoid Pulay stress, the PAW method and the energy cutoff on the VASP Wiki.
5.2 Input¶
The input files to run this example are prepared at $TUTORIALS/bulk/e05_cd-Si-relaxation. Check them out!
cubic diamond Si
5.5
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.130
System = cd Si
ISMEAR = 0
SIGMA = 0.1
ENCUT = 500
IBRION = 2
ISIF = 3
NSW = 5
EDIFFG = -0.001
EDIFF = 1E-06
K-Points
0
Monkhorst Pack
11 11 11
0 0 0
Pseudopotential of Si.
Click to see the answer!
By performing a DFT calculation with ISIF = 4 and IBRION = 2 after each calculation where the volume is changed!
5.3 Calculation¶
Click to see the answer!
The optimal lattice constant in Ångström is $5.5 \cdot 0.4972118640938737 / 0.5=5.469$. The calculation stops because NSW = 5. So you need to copy CONTCAR to replace POSCAR and restart the calculation.
Pulay stress appears as negative external pressure in the OUTCAR file and is negligibly small, because the cutoff energy is set high enough.
Mind that, while it is convenient to submit only one job, changing so many degrees of freedom simultaneously makes this procedure unstable. It is more robust to change the volume as done in Example 4 and perform additional calculations with ISIF = 4 and IBRION = 2 to relax the cell shape and ionic position. What's more, using the ISIF tag, the volume is relaxed at a constant basis, rather than a constant energy cutoff. If the volume or cell shape changes significantly, the basis needs to be changed as well, which is possible by restarting the calculation from the intermittent result of the geometry relaxation. That is the POSCAR file needs to be replaced by the CONTCAR file.
Visualize the relaxed structure using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e05_cd-Si-relaxation" )
mycalc.structure.plot()
5.4 Questions¶
- For which algorithm used to relax the ionic degrees of freedom during a volume optimization is Pulay stress a source of error?
- What reason could there be not to relax the lattice constant, cell shape and ionic positions in a single VASP run?
- What does the ISIF tag control?
By the end of this tutorial, you will be able to:
- stop a calculation by use of a STOPCAR file
- restart an ionic relaxation
6.1 Task¶
Relax ionic positions of a perturbed cd Si structure using the conjugate-gradient algorithm.
Depending on your system, you will most likely run VASP jobs at a high-performance-computing (hpc) center or at a hpc cluster at your institution. When you run heavy jobs, that take a lot of numerical effort, it can be convenient to prompt a command to control VASP without forcefully aborting the job so that the intermittent results are not lost. This can be achieved by adding files to the calculation's directory.
What's more, restarting a job from intermittent results can be important, not just when interrupting a job, but also if you want to improve a calculation's precision or compute additional quantities on top of a converged result.
Read more about aborting and restarting calculations in the following articles: STOPCAR, ICHARG and ISTART.
6.2 Input¶
The input files to run this example are prepared at $TUTORIALS/bulk/e06_perturbed-cd-Si. Check them out!
perturbed cd Si
5.5
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.130
System = perturbed cd Si
ISTART = 0
ICHARG = 2
ISMEAR = 0
SIGMA = 0.1
ENCUT = 360
IBRION = 2
ISIF = 2
NSW = 10
EDIFFG = -0.0001
EDIFF = 1E-06
K-Points
0
Monkhorst Pack
11 11 11
0 0 0
Pseudopotential of Si.
In the POSCAR file, one coordinate of the second Si atom is slightly perturbed.
Add comments after each tag in the INCAR file, in order to recall what they control! Which tag enables the ionic position to be changed? How many ionic steps will this run perform at most? What is the stopping criterion?
6.3 Calculation¶
Go ahead and run VASP:
cd $TUTORIALS/bulk/e06_*
mpirun -np 2 vasp_std
The stdout will show details about the calculation. What stopping criterion lead to writing the final output?
Click to see the answer!
The calculation stops because of NSW = 10 and did not reach full convergence.
Run a fresh calculation avoiding this criterion by changing the input so that VASP would continue to run for more iterations. Then, write a STOPCAR file to this example's directory so that you exit after a completed ionic step at approximately the same iteration as before! Hint: Open a second terminal window and navigate to this example's directory.
Click to see the answer!
Change the INCAR file as follows.
INCAR
System = perturbed cd Si
ISTART = 0
ICHARG = 2
ISMEAR = 0
SIGMA = 0.1
ENCUT = 360
IBRION = 2
ISIF = 2
NSW = 15
EDIFFG = -0.0001
EDIFF = 1E-06
Then, open two terminals and navigate to this example's directory! In one terminal start the calculation by entering the following.
cd $TUTORIALS/bulk/e06_*
vasp_rm
mpirun -np 4 vasp_std
In the other terminal, write a STOPCAR file containing the command to stop before entering the next ionic step by entering the following.
cd $TUTORIALS/bulk/e06_*
echo LSTOP = .TRUE. > STOPCAR
Click to see the answer!
Stopping criteria are always arbitrary to some extent. It is therefore important that you take control of it, and reflect upon why your calculation stopped and how you could continue it. If for instance, you managed to stop your calculation after 3 ionic iteration steps using the STOPCAR file, then you would need to perform 7 more steps to reach the same number of ionic steps. In order to set the stopping criteria to slightly better values as you have already reached, you could set NSW = 8 in the INCAR file.
Then, enter the following into the terminal:
cd $TUTORIALS/bulk/e06_*
rm STOPCAR
cp CONTCAR POSCAR
mpirun -np 2 vasp_std
Notice to what degree the crystal symmetry can be restored in the converged calculation.
6.4 Questions¶
- What is the meaning of a negative EDIFFG?
- Can you use the wavefunction and charge density of a calculation that was stopped using the STOPCAR file?
- The content of which file needs to replace the content of POSCAR to continue an ionic relaxation?
By the end of this tutorial, you will be able to:
- choose the appropriate $\mathbf{k}$-mesh density
- create a POTCAR file for single-element calculations
7.1 Task¶
Perform a convergence study for the relaxation of the volume with respect to the $\mathbf{k}$-mesh density and obtain the DOS for the $\beta$-tin structure of silicon.
The choice of the $\mathbf{k}$ mesh mainly depends on the structure defined in the POSCAR file, but it also depends on the specific system, e.g. the presence of long-range interaction, and the property of interest. To set the appropriate $\mathbf{k}$-mesh density for your calculation, you need to perform a convergence study.
7.2 Input¶
Check out the input files at $TUTORIALS/bulk/e07_beta-tin-Si!
beta-tin Si
4.8
1.0 0.0 0.0
0.0 1.0 0.0
0.5 0.5 0.26
2
Direct
-0.125 -0.375 0.25
0.125 0.375 -0.25
Recall that, the POTCAR file contains the pseudopotential, the solution of the nonrelativistic Schrödinger equation for the non-spinpolarized reference atom including the single-particle wavefunctions and charge densities, as well as 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 file at $TUTORIALS/bulk/e07_beta-tin-Si/pot. Follow the instructions on the VASP Wiki on how to create a POTCAR file for your calculation!
Click to see the answer!
cd $TUTORIALS/bulk/e07_*
cat pot/Si/POTCAR > POTCAR
The unit cell volume in the provided POSCAR file cannot be trusted, so you will need to relax the structure. Search on the web to find the defining features of a $\beta$-tin structure of silicon and relax the volume at different $\mathbf{k}$-mesh density!
Create your own INCAR and KPOINTS file depending on the task.
7.3 Calculation¶
First, find the optimal cell volume by repeatedly starting a relaxation with ISIF = 7 for different $\mathbf{k}$-mesh densities up to a precision of approximately 0.005 Å in the length of the first lattice vector and 2 meV in the total energy. In other words, plot both, the length of the first lattice vector and the total energy, as a function of the number of $\mathbf{k}$ points. Finally, change the tags appropriate for a DOS calculation. After you have prepared the input, go ahead and run VASP for each task:
cd $TUTORIALS/bulk/e07_*
mpirun -np 2 vasp_std
Click to see one possible solution!
Create input for the volume relaxation:
cp POSCAR POSCAR.bk
System = volume relaxation beta-tin Si
ISMEAR = 0
SIGMA = 0.1
ENCUT = 500
IBRION = 2
ISIF = 7
NSW = 6
EDIFFG = -0.001
EDIFF = 1E-06
loop.sh
for k in 03 04 05 06 07 08 09 10 11 12 13
do
cp POSCAR.bk POSCAR
vasp_rm
cat > KPOINTS << EOF
K-Points
0
Monkhorst Pack
$k $k $k
0 0 0
EOF
mpirun -np 2 vasp_std
cp CONTCAR POSCAR
mpirun -np 2 vasp_std
cp CONTCAR POSCAR.k$k
v=$(awk 'NR==3' CONTCAR)
en=$(tail -n 1 OSZICAR)
echo $k $en $v >> loop.dat
done
Then, plot the result with the script loop.gp.
set term png
set output "loop_energy.png"
set title "Convergence study beta-tin Si"
set xlabel "no of k points"
set ylabel "Total energy (eV)"
plot "loop.dat" using 1:4 w lp
set output "loop_lattice.png"
set ylabel "length of first lattice vector (Angstrom)"
plot "loop.dat" using 1:(4.8*$10) w lp
and open the files loop_energy.png and loop_lattice.png.
Notice that increasing the number of $\mathbf{k}$ points does not lead to a monotonous trend towards a certain value, but certainly, the variance decreases with increasing number of $\mathbf{k}$ points. What's more, a significant difference is introduced by either using a Γ-centered grid or not. For information how to choose a Γ-centered grid read about the KPOINTS file on the VASP Wiki!
Use py4vasp to visualize multiple unit cells of the converged structure!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e07_beta-tin-Si" )
mycalc.structure.plot(2)
Recall that the DOS can be obtained by adding LORBIT = 11 to the INCAR file.
System = beta-tin Si, Gaussian smearing
ISMEAR = 0
SIGMA = 0.1
ENCUT = 500
EDIFF = 1E-06
LORBIT = 11
Start from the converged structure and increase the $\mathbf{k}$-point mesh density to $18\times18\times18$. After the VASP calculation, back up the resulting vaspout.h5 file.
cd $TUTORIALS/bulk/e07_*
cp CONTCAR POSCAR
vasp_rm
mpirun -np 2 vasp_std
cp vaspout.h5 backup.h5
Compare the resulting DOS with one obtained with the tetrahedron method with Blöchl correction!
System = beta-tin Si, tetrahedron method
ISMEAR = -5
ENCUT = 500
EDIFF = 1E-06
LORBIT = 11
cd $TUTORIALS/bulk/e07_*
mpirun -np 2 vasp_std
Use py4vasp to plot the different DOS in a single figure.
import py4vasp
gaussian = py4vasp.Calculation.from_file("./e07_beta-tin-Si/backup.h5" )
tetrahedron = py4vasp.Calculation.from_path("./e07_beta-tin-Si/")
gaussian.dos.plot().label("Gaussian") + tetrahedron.dos.plot().label("Tetrahedron")
7.4 Questions¶
- Does the total energy show a monotonous behavior for increasing the number of $\mathbf{k}$ points?
- How can you determine the appropriate number of $\mathbf{k}$ points for your calculation?