fcc Si

Based on the VASP wiki example in this link

Task: Run a self-consistent calculation for fcc Si. Thereafter, optimize the lattice constant. Use some basic tools and scripts.

  • An introductory example for working with VASP

  • In general, examples are chosen for fast calculation

  • Use different job scripts, submit to the Slurm job scheduler on Tetralith or LEONARDO

  • Use VASP, gnuplot, python, ASE

  • bash shell scripts

System-specific instructions

Select instructions for the system you are using:

Instructions for use on the NAISS cluster Tetralith (NSC)

First, copy the example folder which contains the input files INCAR, POSCAR, KPOINTS and some useful scripts

cp -r /software/sse2/tetralith_el9/manual/vasp/training/ws2024/fcc_Si .
cd fcc_Si

also copy the latest POTCAR (PBE GGA) file for Si

cp /software/sse2/tetralith_el9/manual/vasp/POTCARs/PBE/2024-03-19/Si/POTCAR .

Input files

POSCAR

fcc Si
 3.9
 0.5 0.5 0.0
 0.0 0.5 0.5
 0.5 0.0 0.5
Si
   1
cartesian
0 0 0

  • Fcc Si lattice constant a=3.9 Å

  • “Si” on line 6 isn’t needed for the run, but for viewing e.g. using VESTA and practical as a reminder if you have several elements. Note that POTCAR decides the actual order of the elements in POSCAR, check with grep PAW POTCAR

  • 1 atom in the cell

INCAR

System = fcc Si

ISTART = 0       
ICHARG = 2       
ENCUT  = 240     
ISMEAR = 0       
SIGMA  = 0.1

  • ISMEAR determines the partial occupancies of each orbital. For band gap systems -5 (tetrahedron method with Blöchl corrections) or 0 (Gaussian smearing) are suitable.

  • For very accurate total energy calculations or DOS, use ISMEAR = -5

KPOINTS

k-points
 0
Monkhorst Pack
 11 11 11
 0  0  0

  • Face-centered-cubic structure, so use an equally spaced k-mesh

  • Odd number of k-points in each direction -> Gamma centered mesh

  • Alternatively, use the Gamma centered k-mesh, replacing the line “Monkhorst Pack” (M) with “Gamma” (G)

POTCAR

  PAW_PBE Si 05Jan2001                   
   4.00000000000000     
 parameters from PSCTR are:
   VRHFIN =Si: s2p2
   LEXCH  = PE
   EATOM  =   103.0669 eV,    7.5752 Ry

   TITEL  = PAW_PBE Si 05Jan2001
   LULTRA =        F    use ultrasoft PP ?
   IUNSCR =        1    unscreen: 0-lin 1-nonlin 2-no
   RPACOR =    1.500    partial core radius
   POMASS =   28.085; ZVAL   =    4.000    mass and valenz
   RCORE  =    1.900    outmost cutoff radius
   RWIGS  =    2.480; RWIGS  =    1.312    wigner-seitz radius (au A)
   ENMAX  =  245.345; ENMIN  =  184.009 eV
   ICORE  =        2    local potential
   LCOR   =        T    correct aug charges
   LPAW   =        T    paw PP
   EAUG   =  322.069
   DEXC   =    0.000
   RMAX   =    1.950    core radius for proj-oper
   RAUG   =    1.300    factor for augmentation sphere
   RDEP   =    1.993    radius for radial grids
   RDEPT  =    1.837    core radius for aug-charge

   Atomic configuration
    6 entries
     n  l   j            E        occ.
     1  0  0.50     -1785.8828   2.0000
     2  0  0.50      -139.4969   2.0000
     2  1  1.50       -95.5546   6.0000
     3  0  0.50       -10.8127   2.0000
     3  1  0.50        -4.0811   2.0000
     3  2  1.50        -4.0817   0.0000
...

  • The top lines of the POTCAR file

  • Never (almost) edit

Check that correct POTCARs are used and the ENMAX values for setting ENCUT in INCAR

grep PAW POTCAR
grep ENMAX POTCAR

1. Calculation

A first static calculation for Si with the lattice constant 3.9 Å. Create a new folder “scf0”, move the input files + job script, copy the folder to “scf1” (we will use it later) and go to “scf0”

mkdir scf0
mv INCAR POSCAR KPOINTS POTCAR run.sh scf0
cp -r scf0 scf1
cd scf0

Submit the job script “run.sh” to the job queue with

sbatch run.sh

The job script looks like below

#!/bin/bash
#SBATCH -A naiss2024-22-241
#SBATCH -t 0:30:00
#SBATCH -n 4 
#SBATCH -J vaspjob

module load VASP/6.4.3.19032024-omp-hpc1-intel-2023a-eb
mpprun vasp_std

The job should finish quite quickly. Check the status of the job (if submitted to the queue) with

squeue -u USERNAME

it might already have finished, in that case it’ll not be shown. Among the output files in particular note

OUTCAR
OSZICAR
slurm-JOBID.out

in which JOBID is a string of numbers. Note that

  • Some warning messages will only show in the standard output of the job, slurm-JOBID.out if submitted to the queue. Therefore, also carefully check the output when running interactively.

You can quickly check the iteration steps and convergence with “cat”

cat OSZICAR

For checking into the main output file OUTCAR and other large files, it’s convenient to use the fast and safe (no writing) text viewer “less”

less OUTCAR

“less” can be quit with pressing q.

2. Total energy vs volume calculation

Go back to the main folder “fcc_Si” with

cd ..
pwd

“pwd” shows the present folder, confirm that you’re back in “fcc_Si”, otherwise use “cd” to go there.

In this part we will calculate the total energies of fcc Si between 3.5 and 4.3 Å. To do this a bit quicker, a special job script “run-vol.sh” is prepared (bash shell script):

#!/bin/bash
#SBATCH -A naiss2024-22-241
#SBATCH -t 0:30:00
#SBATCH -n 4
#SBATCH -J vaspjob

module load VASP/6.4.3.19032024-omp-hpc1-intel-2023a-eb

for i in  3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 ; do
mkdir -p $i
cd $i
cp /software/sse2/tetralith_el9/manual/vasp/POTCARs/PBE/2024-03-19/Si/POTCAR .
cp /software/sse2/tetralith_el9/manual/vasp/training/ws2024/fcc_Si/INCAR .
cp /software/sse2/tetralith_el9/manual/vasp/training/ws2024/fcc_Si/KPOINTS .
cat >POSCAR <<EOF
fcc:
   $i
 0.5 0.5 0.0
 0.0 0.5 0.5
 0.5 0.0 0.5
Si
   1
cartesian
0 0 0
EOF
mpprun vasp_std
E=`awk '/F=/ {print $0}' OSZICAR` ; echo $i $E  >>../SUMMARY.fcc
cd ..
done

In brief, the above script creates a new folder with the same name as the lattice constant, copying POTCAR, INCAR and KPOINTS, while creating a new POSCAR file. It also collects the total energies from all OSZICAR files into a new file “SUMMARY.fcc”.

Now, submit the job script “run-vol.sh” to the queue

sbatch run-vol.sh

After a successful run, there should now be folders 3.5 - 4.3 with the different calculations and a file SUMMARY.fcc which can be checked e.g. with “less” or “cat”

3.5 1 F= -.44256909E+01 E0= -.44234194E+01 d E =-.454310E-02
3.6 1 F= -.46614931E+01 E0= -.46600638E+01 d E =-.285864E-02
3.7 1 F= -.47980122E+01 E0= -.47959555E+01 d E =-.411323E-02
3.8 1 F= -.48645323E+01 E0= -.48630343E+01 d E =-.299612E-02
3.9 1 F= -.48774144E+01 E0= -.48758835E+01 d E =-.306173E-02
4.0 1 F= -.48487753E+01 E0= -.48481407E+01 d E =-.126921E-02
4.1 1 F= -.47852971E+01 E0= -.47845193E+01 d E =-.155573E-02
4.2 1 F= -.46937300E+01 E0= -.46922884E+01 d E =-.288335E-02
4.3 1 F= -.45831536E+01 E0= -.45812206E+01 d E =-.386597E-02

You can make a quick plot of the total energy as a function of the lattice constant, e.g. by using gnuplot, start it with

gnuplot

thereafter, at the “gnuplot>” prompt type

plot "SUMMARY.fcc" using ($1):($4) w lp

From the plot we can see that the equilibrium lattice constant is around 3.9 Å.

To find a more exact value, we can use an equation of state method. For example, a popular method is Birch-Murnaghan. To quickly check, you can try a python script which uses functions from the Atomic Simulation Environment (ASE) collection of tools based on Python3.

First, we need to load ASE and a suitable Python3 module,

module load ASE/3.22.1-hpc1-python
module load Python/3.10.4-env-hpc1-gcc-2022a-eb

As input a file with volume (Å^3) and total energy (eV) is needed. To save some time, there is a script which can be used, so run:

./get-vol-etot.sh

The bash script “get-vol-etot.sh” looks like

#!/bin/bash
for i in  3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 ; do
#V = `grep volume $i/vasprun.xml | head -1 | awk '{print $3}'`
cd $i
V=`awk '/volume/ {print $3; exit}' vasprun.xml`
E=`awk '/F=/ {print $3}' OSZICAR`  
echo $V $E  >>../vol_etot.dat
cd ..
done

and produces the file “vol_etot.dat”. We will use a small python script “eqos.py” which uses ASE and looks like

import numpy
from ase.units import GPa
from ase.eos import EquationOfState

# Read data
volumes,energies = numpy.loadtxt('vol_etot.dat', unpack=True)

# Fit EOS
eos = EquationOfState(volumes,energies, eos='birchmurnaghan')

# There are several other choices for eos, see doc at:
# https://wiki.fysik.dtu.dk/ase/ase/eos.html
# set eos to e.g. birchmurnaghan, sjeos (ASE default), p3

v0,e0,B = eos.fit()

print ()
print ("Equation of state:")
print ("-------------------")
print ("E0: %2.6f eV" % (e0))
print ("V0: %2.6f A^3" % (v0))
print ("B: %2.6f GPa" % (B/GPa))
print ()

run it to calculate the equilibrium volume using an equation of state (first loading the modules or activating the environment, as described above)

python eqos.py

the output might look like

Equation of state:
-------------------
E0: -4.880229 eV
V0: 14.546730 A^3
B: 84.365350 GPa

The EOS methods are expected to work better closer to the equilibrium. The chosen method can be changed by editing eos='birchmurnaghan' in the line

eos = EquationOfState(volumes,energies, eos='birchmurnaghan')

Note that the volume V0 is for a single atom, so to obtain the lattice constant a, one needs to calculate (since there are 4 atoms per fcc cell)

a = (4 x V0)^1/3

To quickly compute a, you might e.g. use python directly in the terminal

python

and input at the “>>>” prompt

(4*14.546730) ** (1./3.)

and exit with ctrl d (what happens adding 0.7 + 0.6 using python?).

3. Calculation

A second static calculation, now using the obtained equilibrium volume. Go to the folder “scf1” which was created previously

cd scf1

and insert the lattice parameter a obtained in 2. by editing POSCAR, e.g. with

vi POSCAR

or use your favourite text editor (emacs, nano, gedit, …). Finally, submit the job

sbatch run.sh

and compare the total energy with previous results. For example, try

grep "free  energy" OUTCAR
grep "free energy" OUTCAR
grep F= OSZICAR
cat OSZICAR
grep F= slurm-JOBID.out
cat slurm-JOBID.out

  • Comparing with the EOS result, does it look reasonable?

Extra tasks

  • Try different methods to compute the equilibrium by editing the script “eqos.py”

  • Compute the equilibrium volumes for a less dense k-mesh 5 5 5 in KPOINTS and with ENCUT set to ENMIN from POTCAR, and (2) for a denser k-mesh 29 29 29 and with ENCUT = ENMAX x 1.3. Suggestion: create new folders, copy and edit the “run_vol.sh” job script