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
a0 = 3.52 / np.sqrt(2)
c0 = np.sqrt(8 / 3.0) * a0
and create a trajectory for the results:
from ase.io import Trajectory
traj = Trajectory('Ni.traj', 'w')
Finally, we do the 9 calculations (three values for \(a\) and three for \(c\)):
from ase.build import bulk
from ase.calculators.emt import EMT
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]
:
from ase.io import read
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)
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)