Equation of state tutorial

Requirements

Software components

  • nanotools

  • RESCU+

Pseudopotentials

I will need the following pseudopotential file

  • C_AtomicData.mat

Let’s copy it from the pseudo database to the current working directory and export RESCUPLUS_PSEUDO=$PWD.

Briefing

In this tutorial, I show how to calculate the lattice constant of gallium arsenide via equation of state.

Code

from nanotools import Atoms, Cell, System, TotalEnergy
from nanotools.eos import EquationOfState as EOS
import numpy as np
a = 3.562 / 2.
cell = Cell(avec=[[0.,a,a],[a,0.,a],[a,a,0.]], resolution=0.12)
fxyz = [[0.00,0.00,0.00],[0.25,0.25,0.25]]
atoms = Atoms(fractional_positions=fxyz, formula="C(2)")
sys = System(cell=cell, atoms=atoms)
ecalc = TotalEnergy(sys)
eos = EOS.from_totalenergy(ecalc)
eos.relative_volumes = np.linspace(0.95, 1.05, 5)
eos.solve()
eos.plot_eos(filename="eos-C.png")

Explanations

I will compute the equilibrium volume of diamond by fitting energy versus volume data to a Birch-Murnaghan equation of state.

\[P = \frac{3}{2} K_0 \left[\left(\frac{V}{V_0}\right)^{-\frac{7}{3}} - \left(\frac{V}{V_0}\right)^{-\frac{5}{3}}\right]\left[1 + \frac{3}{4}\left(K_0' - 4\right)\left[\left(\frac{V}{V_0}\right)^{-\frac{2}{3}}-1\right]\right]\]

Here is a high level view of the calculation workflow:

  1. Create a TotalEnergy calculator and specify calculation parameters like resolution, k-sampling, tolerance, etc.

  2. Create an EquationOfState calculator.

  3. Solve the Kohn-Sham equation for a few volumes to obtain the EOS.

Total energy calculator

I first define a diamond crystal system with a lattice constant of 3.562 Ang. I then create a total energy calculator (ecalc) from this System object. For more detail on how to setup a total energy calculation, refer to total energy tutorial.

Equation of state calculator

I can then initialize a EOS calculator from the total energy calculator. I use the numpy function linspace to create an array of relative volumes which I set in relative_volumes. The line

eos.relative_volumes = np.linspace(0.95, 1.05, 5)

is equivalent to

eos.relative_volumes = [0.95 , 0.975, 1.   , 1.025, 1.05 ]

Upon calling solve, the EOS calculator then scales the volume of the original TotalEnergy calculator keeping the fractional coordinates fixed, and it computes the total energy of each crystal. The calculator will save results to nano_eos_out.json as the solving progresses. The results can then be loaded with the class method read. You can finally see the data and EOS fit typing

fig = calc.plot_eos()

which will produce

EOS fit for diamond