.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/findlatticeconstants/04-lattice.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_generated_findlatticeconstants_04-lattice.py: Finding lattice constants using EOS and the stress tensor ######################################################### .. GENERATED FROM PYTHON SOURCE LINES 8-18 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 20-23 Finding lattice constants using the equation of state ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 25-30 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` by calculating :math:`a^3 = V`, where :math:`V` is the volume of the unit cell. .. GENERATED FROM PYTHON SOURCE LINES 32-40 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 ` 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 40-52 .. 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 53-54 We create a trajectory to save the results: .. GENERATED FROM PYTHON SOURCE LINES 54-57 .. code-block:: Python traj = Trajectory('Ni.traj', 'w') .. GENERATED FROM PYTHON SOURCE LINES 58-60 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 60-68 .. 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 69-72 Analysis -------- Now, we need to extract the data from the trajectory. Try this: .. GENERATED FROM PYTHON SOURCE LINES 72-75 .. 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 76-78 So, we can get :math:`a` and :math:`c` from ``ni.cell[0, 0]`` and ``ni.cell[2, 2]``: .. GENERATED FROM PYTHON SOURCE LINES 78-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-121 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 121-125 .. code-block:: Python ni = bulk('Ni', 'hcp', a=3.0, c=5.0) ni.calc = EMT() .. GENERATED FROM PYTHON SOURCE LINES 126-132 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. .. GENERATED FROM PYTHON SOURCE LINES 132-143 .. 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 14:41:05 2.697838 11.456177 BFGS: 1 14:41:05 0.885752 8.925609 BFGS: 2 14:41:05 0.094401 4.475456 BFGS: 3 14:41:05 -0.007973 1.687704 BFGS: 4 14:41:05 -0.029392 0.206563 BFGS: 5 14:41:05 -0.029737 0.023818 BFGS: 6 14:41:05 -0.029743 0.006381 BFGS: 7 14:41:05 -0.029743 0.000020 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 144-146 Analyze the results print the unit cell .. GENERATED FROM PYTHON SOURCE LINES 146-148 .. 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_findlatticeconstants_04-lattice.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04-lattice.ipynb <04-lattice.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04-lattice.py <04-lattice.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04-lattice.zip <04-lattice.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_