Finding lattice constants using EOS and the stress tensor¶

HCP¶

Let’s try to find the $$a$$ and $$c$$ lattice constants for HCP nickel using the EMT potential.

First, we make a good initial guess for $$a$$ and $$c$$ using the FCC nearest neighbor distance and the ideal $$c/a$$ ratio:

import numpy as np

from ase.io import Trajectory

and create a trajectory for the results:

from ase.calculators.emt import EMT

Finally, we do the 9 calculations (three values for $$a$$ and three for $$c$$):

a0 = 3.52 / np.sqrt(2)
c0 = np.sqrt(8 / 3.0) * a0

traj = Trajectory('Ni.traj', 'w')

eps = 0.01
for a in a0 * np.linspace(1 - eps, 1 + eps, 3):
for c in c0 * np.linspace(1 - eps, 1 + eps, 3):
ni = bulk('Ni', 'hcp', a=a, c=c)

Analysis¶

Now, we need to extract the data from the trajectory. Try this:

>>> from ase.build import bulk
>>> ni = bulk('Ni', 'hcp', a=2.5, c=4.0)
>>> ni.cell
array([[ 2.5  ,  0.   ,  0.   ],
[-1.25 ,  2.165,  0.   ],
[ 0.   ,  0.   ,  4.   ]])

So, we can get $$a$$ and $$c$$ from ni.cell[0, 0] and ni.cell[2, 2]:

ni.get_potential_energy()
traj.write(ni)

energies = [config.get_potential_energy() for config in configs]
a = np.array([config.cell[0, 0] for config in configs])

We fit the energy to this expression:

$p_0 + p_1 a + p_2 c + p_3 a^2 + p_4 ac + p_5 c^2$

The best fit is found like this:

c = np.array([config.cell[2, 2] for config in configs])

and we can find the minimum like this:

p = np.linalg.lstsq(functions.T, energies, rcond=-1)

p0 = p
p1 = p[1:3]
p2 = np.array([(2 * p, p),

Results:

a

c

2.466

4.023

Using the stress tensor¶

One can also use the stress tensor to optimize the unit cell. For this we cannot use the EMT calculator.:

from ase.optimize import BFGS
from ase.constraints import StrainFilter
from gpaw import GPAW, PW
ni = bulk('Ni', 'hcp', a=a0,c=c0)
calc = GPAW(mode=PW(200),xc='LDA',txt='Ni.out')
ni.calc = calc
sf = StrainFilter(ni)
opt = BFGS(sf)
opt.run(0.005)

If you want the optimization path in a trajectory, add these lines before calling the run() method:

traj = Trajectory('path.traj', 'w', ni)
opt.attach(traj)