
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples_generated/tutorials/lattice.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_generated_tutorials_lattice.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_generated_tutorials_lattice.py:

.. _lattice_constant_example:

Lattice constants using EOS and the stress tensor
=================================================

.. GENERATED FROM PYTHON SOURCE LINES 9-19

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.

.. GENERATED FROM PYTHON SOURCE LINES 21-24

Finding lattice constants using the equation of state
-----------------------------------------------------


.. GENERATED FROM PYTHON SOURCE LINES 26-31

FCC
^^^
The lattice constant :math:`a` for FCC bulk metal can be obtained using the
equation of state as outline in the tutorial :ref:`eos_example` by calculating
:math:`a^3 = V`, where :math:`V` is the volume of the unit cell.

.. GENERATED FROM PYTHON SOURCE LINES 33-41

HCP
^^^
For HCP bulk metals, we need to account for two lattice constants, :math:`a`
and :math:`c`. Let's try to find :math:`a` and :math:`c` for HCP nickel
using the :mod:`EMT <ase.calculators.emt>` potential.

First, we make an initial guess for :math:`a` and :math:`c` using the FCC
nearest neighbor distance and the ideal :math:`c/a` ratio:

.. GENERATED FROM PYTHON SOURCE LINES 41-53

.. code-block:: Python


    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








.. GENERATED FROM PYTHON SOURCE LINES 54-55

We create a trajectory to save the results:

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: Python


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








.. GENERATED FROM PYTHON SOURCE LINES 59-61

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

.. GENERATED FROM PYTHON SOURCE LINES 61-69

.. code-block:: Python

    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)








.. GENERATED FROM PYTHON SOURCE LINES 70-73

Analysis
--------
Now, we need to extract the data from the trajectory.  Try this:

.. GENERATED FROM PYTHON SOURCE LINES 73-76

.. code-block:: Python

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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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




.. GENERATED FROM PYTHON SOURCE LINES 77-79

So, we can get :math:`a` and :math:`c` from ``ni.cell[0, 0]``
and ``ni.cell[2, 2]``:

.. GENERATED FROM PYTHON SOURCE LINES 79-83

.. code-block:: Python

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







.. GENERATED FROM PYTHON SOURCE LINES 84-91

We fit the energy :math:`E` to an expression for the equation of state,
e.g., a polynomial:

.. math:: E = p_0 + p_1 a + p_2 c + p_3 a^2 + p_4 ac + p_5 c^2

:math:`p_i` are the parameters. We can find the parameters for the
best fit, e.g., through least squares fitting:

.. GENERATED FROM PYTHON SOURCE LINES 91-94

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 95-96

then, we can find the minimum like this:

.. GENERATED FROM PYTHON SOURCE LINES 96-109

.. code-block:: Python

    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} Å')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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




.. GENERATED FROM PYTHON SOURCE LINES 110-138

.. code-block:: Python


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



.. image-sg:: /examples_generated/tutorials/images/sphx_glr_lattice_001.png
   :alt: Interpolated energy, eV
   :srcset: /examples_generated/tutorials/images/sphx_glr_lattice_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 139-150

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 :math:`a` and :math:`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.

.. GENERATED FROM PYTHON SOURCE LINES 150-154

.. code-block:: Python


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








.. GENERATED FROM PYTHON SOURCE LINES 155-161

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
<https://en.wikipedia.org/wiki/Broyden–Fletcher–Goldfarb–Shanno_algorithm>`_
algorithm is allowed to modify lattice parameters to minimise target
stress.

.. GENERATED FROM PYTHON SOURCE LINES 161-172

.. code-block:: Python


    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/Å³





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

          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_



.. GENERATED FROM PYTHON SOURCE LINES 173-175

Analyze the results
print the unit cell

.. GENERATED FROM PYTHON SOURCE LINES 175-177

.. code-block:: Python

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




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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





.. _sphx_glr_download_examples_generated_tutorials_lattice.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: lattice.ipynb <lattice.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: lattice.py <lattice.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: lattice.zip <lattice.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
