How to do volume optimization using ASE interface

 

The ASE (Atomic Simulation Environment)  interface s a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations. To install it click on the following link:

How to download and install the Atomic Simulation Environment ASE

The ASE interface create the input files and do the calculation using only a python script.

We will take the example of Magnesium

First we need to add the path for vasp executable and path for pseudopotentials as follows:

#~-----vasp-ase path ---------
export VASP_COMMAND="/home/algerien1970/abinitio/vasp.5.4.1/bin/vasp_std"
export VASP_PP_PATH=/home/algerien1970/abinitio/vasp.5.4.1/POTCARs 

 

The directoriy POTCARs will contain 3 subdirectories named: 

potpaw  for LDA pseudopotentials

potpaw_PBE   for PBE pseudopotentials

potpaw_GGA  for PW91 pseudopotentials

Calculation

algerien1970@linux-kq57:~/abinitio/vasp-ase-tutorials/Si> conda activate ase 

(ase) algerien1970@linux-kq57:~/abinitio/vasp-ase-tutorials> mkdir mg-hcp 

(ase) algerien1970@linux-kq57:~/abinitio/vasp-ase-tutorials> cd mg-hcp 

 

mg-hcp.py script

import numpy
from ase.build import bulk
from ase.calculators.vasp import Vasp
from ase.units import GPa
from ase.eos import EquationOfState

jobs = []

# Generate cell configurations
for volume in numpy.linspace(0.90,1.10,10):
  mg = bulk('Mg','hcp',a=3.1,covera=1.62)
  mg.set_cell(mg.get_cell()*volume,scale_atoms=True)

  calc = Vasp(prec = 'Accurate', 
            xc = 'PBE', 
            lreal = False,
            encut = 150.0,
            istart =0,
            icharg= 2,
            nelm = 20,
            kpts = [6,6,6],
            ismear = 1,
            sigma = 0.2,
            algo = "fast")

  mg.set_calculator(calc)
  
  jobs.append(mg)

# Run jobs
volumes = [job.get_volume() for job in jobs]
energies = [job.get_potential_energy() for job in jobs]

print ("Calculation output")
print ("------------------")
print (("Volumes:"), volumes)
print (("Energies:"), energies)
print

# Fit EOS
eos = EquationOfState(volumes,energies)

v0,e0,B = eos.fit()

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

(ase) algerien1970@linux-kq57:~/abinitio/vasp-ase-tutorials/mg-hcp> chmod +x mg-hcp.py

(ase) algerien1970@linux-kq57:~/abinitio/vasp-ase-tutorials/mg-hcp> python mg-hcp.py  

Calculation output
------------------
Volumes: [30.469003876435877, 32.782153276649034, 35.2095096268752, 37.753824901813466, 40.417851076162975, 43.20434012462278, 46.11604402189205, 49.155714742669836, 52.326104261655324, 55.62996455354755]
Energies: [-1.86460711, -2.27470326, -2.57574479, -2.78654683, -2.92182514, -2.99416926, -3.01374962, -2.9902221, -2.93151641, -2.84418069]
Equation of state parameters
----------------------------
E0: -3.013844 eV
V0: 45.9 A^3
B: 36.5 GPa