Finding lattice constants using EOS and the stress tensor

Finding lattice constants using EOS and the stress tensor#

Note

We are currently moving to a new way to display our examples. For this example we have an updated version, which you can find here. The example on this page is deprecated and will be removed once all examples have been moved to the new format.

See also

Equation of state (EOS).

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.build import bulk
from ase.calculators.emt import EMT
from ase.io import Trajectory, read

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

and create a trajectory for the results:


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

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


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)
        ni.calc = EMT()
        ni.get_potential_energy()
        traj.write(ni)

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]:


configs = read('Ni.traj@:')
energies = [config.get_potential_energy() for config in configs]
a = np.array([config.cell[0, 0] for config in configs])
c = np.array([config.cell[2, 2] 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:


functions = np.array([a**0, a, c, a**2, a * c, c**2])
p = np.linalg.lstsq(functions.T, energies, rcond=-1)[0]

and we can find the minimum like this:

p0 = p[0]
p1 = p[1:3]
p2 = np.array([(2 * p[3], p[4]), (p[4], 2 * p[5])])
a0, c0 = np.linalg.solve(p2.T, -p1)

with open('lattice_constant.csv', 'w') as fd:
    fd.write(f'{a0:.3f}, {c0:.3f}\n')

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)
Contents