.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/tutorials/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_tutorials_04-lattice.py: .. _lattice_constant_example: Finding 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` 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 ` 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-84 .. 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 85-92 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 92-95 .. 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 96-97 then, we can find the minimum like this: .. GENERATED FROM PYTHON SOURCE LINES 97-110 .. 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 111-122 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 122-126 .. code-block:: Python ni = bulk('Ni', 'hcp', a=3.0, c=5.0) ni.calc = EMT() .. GENERATED FROM PYTHON SOURCE LINES 127-133 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 133-144 .. 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 08:55:05 2.697838 11.456177 BFGS: 1 08:55:05 0.885752 8.925609 BFGS: 2 08:55:05 0.094401 4.475456 BFGS: 3 08:55:05 -0.007973 1.687704 BFGS: 4 08:55:05 -0.029392 0.206563 BFGS: 5 08:55:05 -0.029737 0.023818 BFGS: 6 08:55:05 -0.029743 0.006381 BFGS: 7 08:55:05 -0.029743 0.000020 np.True_ .. GENERATED FROM PYTHON SOURCE LINES 145-147 Analyze the results print the unit cell .. GENERATED FROM PYTHON SOURCE LINES 147-149 .. 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_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 `_