Skip to main content
`Ctrl`+`K`

# Lattice constants using EOS and the stress tensor

Note

Go to the end to download the full example code.

# Lattice constants using EOS and the stress tensor

# Introduction

The bulk lattice constant for a material modelled with DFT is often different from the experimental lattice constant. For example, for DFT-GGA functionals, lattice constants can be on the order of 4 % larger than the experimental value. To model materials’ properties accurately, it is important to use the optimized lattice constant corresponding the method/functional used.

In this tutorial, we will first look at obtaining the lattice constant by fitting the equation of state and then obtaining it from the stress tensor.

## Finding lattice constants using the equation of state

### FCC

The lattice constant \(a\) for FCC bulk metal can be obtained using the equation of state as outline in the tutorial EOS: Introduction to the Equation of state by calculating \(a^3 = V\), where \(V\) is the volume of the unit cell.

### HCP

For HCP bulk metals, we need to account for two lattice constants, \(a\) and \(c\). Let’s try to find \(a\) and \(c\) for HCP nickel using the `EMT` potential.

First, we make an 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.filters import StrainFilter  # Import for stress tensor calculation
    from ase.io import Trajectory, read
    from ase.optimize import BFGS  # Import for stress tensor calculation

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

We create a trajectory to save the results:

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

Next, we do the 9 calculations (three values for \(a\) and three for \(c\)). Note that for a real-world case, we would want to try more values.

    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:

    ni = bulk('Ni', 'hcp', a=2.5, c=4.0)
    print(ni.cell)

    Cell([[2.5, 0.0, 0.0], [-1.25, 2.1650635094610964, 0.0], [0.0, 0.0, 4.0]])

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 \(E\) to an expression for the equation of state, e.g., a polynomial:

\[E = p_0 + p_1 a + p_2 c + p_3 a^2 + p_4 ac + p_5 c^2\]

\(p_i\) are the parameters. We can find the parameters for the best fit, e.g., through least squares fitting:

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

then, 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)

    # Save the lattice constants to a file
    with open('lattice_constant.csv', 'w') as fd:
        fd.write(f'{a0:.3f}, {c0:.3f}\n')

    # Show the results on screen:
    print('The optimized lattice constants are:')
    print(f'a = {a0:.3f} Å, c = {c0:.3f} Å')

    The optimized lattice constants are:
    a = 2.466 Å, c = 4.023 Å

    import matplotlib.pyplot as plt

    a_int = a0 * np.linspace(1 - eps, 1 + eps, 100)
    c_int = c0 * np.linspace(1 - eps, 1 + eps, 100)

    a_int_grid, c_int_grid = np.meshgrid(a_int, c_int)

    E = (
        p[0]
        + p[1] * a_int_grid
        + p[2] * c_int_grid
        + p[3] * a_int_grid**2
        + p[4] * a_int_grid * c_int_grid
        + p[5] * c_int_grid**2
    )

    plt.imshow(
        E.T,
        cmap='viridis',
        origin='lower',
        extent=[min(c_int), max(c_int), min(a_int), max(a_int)],
    )
    plt.colorbar()
    plt.title('Interpolated energy, eV')
    plt.xlabel('c lattice constant, Å')
    plt.ylabel('a lattice constant, Å')
    plt.show()

# Using the stress tensor

One can also use the stress tensor to optimize the unit cell. Using ‘StrainFilter’, the unit cell is relaxed until the stress is below a given threshold.

Note that if the initial guesses for \(a\) and \(c\) are far from the optimal values, the optimization can get stuck in a local minimum. Similarly, if the threshold is not chosen tight enough, the resulting lattice constants can again be far from the optimal values.

    ni = bulk('Ni', 'hcp', a=3.0, c=5.0)
    ni.calc = EMT()

rather then directly to the atoms. The filter presents an alternative view of the atomic degrees of freedom: instead of modifying atomic positions to minimise target forces, the BFGS algorithm is allowed to modify lattice parameters to minimise target stress.

    sf = StrainFilter(ni)
    opt = BFGS(sf)

    # If you want the optimization path in a trajectory, comment in these lines:
    # traj = Trajectory('path.traj', 'w', ni)
    # opt.attach(traj)

    # Set the threshold and run the optimization::
    opt.run(0.005)  # run until forces are below 0.005 eV/Å³

          Step     Time          Energy          fmax
    BFGS:    0 09:18:50        2.697838       11.456177
    BFGS:    1 09:18:50        0.885752        8.925609
    BFGS:    2 09:18:50        0.094401        4.475456
    BFGS:    3 09:18:50       -0.007973        1.687704
    BFGS:    4 09:18:50       -0.029392        0.206563
    BFGS:    5 09:18:50       -0.029737        0.023818
    BFGS:    6 09:18:50       -0.029743        0.006381
    BFGS:    7 09:18:50       -0.029743        0.000020

    np.True_

Analyze the results print the unit cell

    print('The optimized lattice constants from the stress tensor are:')
    print(f'a = {ni.cell[0, 0]:.3f} Å, c = {ni.cell[2, 2]:.3f} Å')

    The optimized lattice constants from the stress tensor are:
    a = 2.466 Å, c = 4.025 Å

Gallery generated by Sphinx-Gallery