index | modules | gitlab | page source

# Finding lattice constants using EOS and the stress tensor

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
    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.set_calculator(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:

\[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)

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.set_calculator(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)

Documentation v: 3.19.0

Versions
  3.20.0b1 (development)
  3.19.0 (latest stable)