Finding lattice constants using EOS and the stress tensor¶
See also
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:
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)