Finding the volume of a unit cell at a fixed pressure

| categories: uncategorized | tags: | View Comments

A typical unit cell optimization in DFT is performed by minimizing the total energy with respect to variations in the unit cell parameters and atomic positions. In this approach, a pressure of 0 GPa is implied, as well as a temperature of 0K. For non-zero pressures, the volume that minimizes the total energy is not the same as the volume at P=0.

Let $$x$$ be the unit cell parameters that can be varied. For P ≠ 0, and T = 0, we have the following

$$G(x; p) = E(x) + p V(x)$$

and we need to minimize this function to find the groundstate volume. We will do this for fcc Cu at 5 GPa of pressure. We will assume there is only one degree of freedom in the unit cell, the lattice constant. First we get the $$E(x)$$ function, and then add the analytical correction.

from jasp import *
from ase import Atom, Atoms
from ase.utils.eos import EquationOfState

LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]
volumes, energies = [], []

P = 5.0 / 160.2176487  # pressure in eV/ang**3

for a in LC:
atoms = Atoms([Atom('Cu',(0, 0, 0))],
cell=0.5 * a*np.array([[1.0, 1.0, 0.0],
[0.0, 1.0, 1.0],
[1.0, 0.0, 1.0]]))

with jasp('../bulk/Cu-{0}'.format(a),
xc='PBE',
encut=350,
kpts=(8,8,8),
atoms=atoms) as calc:

try:
e = atoms.get_potential_energy()
energies.append(e)
volumes.append(atoms.get_volume())
except (VaspSubmitted, VaspQueued):

import sys; sys.exit()

import numpy as np
energies = np.array(energies)
volumes = np.array(volumes)

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print 'V0 at 0 GPa = {0:1.2f} ang^3'.format(v0)

eos5 = EquationOfState(volumes, energies + P * volumes)
v0_5, e0, B = eos5.fit()
print 'V0 at 5 GPa = {0:1.2f} ang^3'.format(v0_5)

V0 at 0 GPa = 12.02 ang^3
V0 at 5 GPa = 11.62 ang^3


You can see here that apply pressure decreases the equilibrium volume, and increases the total energy.