{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Finding lattice constants using EOS and the stress tensor\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\nThe bulk lattice constant for a material modelled with DFT is often different\nfrom the experimental lattice constant. For example, for DFT-GGA functionals,\nlattice constants can be on the order of 4 % larger than the experimental\nvalue. To model materials' properties accurately, it is important to use\nthe optimized lattice constant corresponding the method/functional used.\n\nIn this tutorial, we will first look at obtaining the lattice constant by\nfitting the equation of state and then obtaining it from the stress tensor.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding lattice constants using the equation of state\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### FCC\nThe lattice constant $a$ for FCC bulk metal can be obtained using the\nequation of state as outline in the tutorial `eos` by calculating\n$a^3 = V$, where $V$ is the volume of the unit cell.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### HCP\nFor HCP bulk metals, we need to account for two lattice constants, $a$\nand $c$. Let's try to find $a$ and $c$ for HCP nickel\nusing the :mod:`EMT ` potential.\n\nFirst, we make an initial guess for $a$ and $c$ using the FCC\nnearest neighbor distance and the ideal $c/a$ ratio:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n\nfrom ase.build import bulk\nfrom ase.calculators.emt import EMT\nfrom ase.filters import StrainFilter # Import for stress tensor calculation\nfrom ase.io import Trajectory, read\nfrom ase.optimize import BFGS # Import for stress tensor calculation\n\na0 = 3.52 / np.sqrt(2)\nc0 = np.sqrt(8 / 3.0) * a0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a trajectory to save the results:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "traj = Trajectory('Ni.traj', 'w')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we do the 9 calculations (three values for $a$ and three for\n$c$). Note that for a real-world case, we would want to try more values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "eps = 0.01\nfor a in a0 * np.linspace(1 - eps, 1 + eps, 3):\n for c in c0 * np.linspace(1 - eps, 1 + eps, 3):\n ni = bulk('Ni', 'hcp', a=a, c=c)\n ni.calc = EMT()\n ni.get_potential_energy()\n traj.write(ni)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Analysis\nNow, we need to extract the data from the trajectory. Try this:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ni = bulk('Ni', 'hcp', a=2.5, c=4.0)\nprint(ni.cell)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, we can get $a$ and $c$ from ``ni.cell[0, 0]``\nand ``ni.cell[2, 2]``:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "configs = read('Ni.traj@:')\nenergies = [config.get_potential_energy() for config in configs]\na = np.array([config.cell[0, 0] for config in configs])\nc = np.array([config.cell[2, 2] for config in configs])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We fit the energy $E$ to an expression for the equation of state,\ne.g., a polynomial:\n\n\\begin{align}E = p_0 + p_1 a + p_2 c + p_3 a^2 + p_4 ac + p_5 c^2\\end{align}\n\n$p_i$ are the parameters. We can find the parameters for the\nbest fit, e.g., through least squares fitting:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "functions = np.array([a**0, a, c, a**2, a * c, c**2])\np = np.linalg.lstsq(functions.T, energies, rcond=-1)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "then, we can find the minimum like this:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p0 = p[0]\np1 = p[1:3]\np2 = np.array([(2 * p[3], p[4]), (p[4], 2 * p[5])])\na0, c0 = np.linalg.solve(p2.T, -p1)\n\n# Save the lattice constants to a file\nwith open('lattice_constant.csv', 'w') as fd:\n fd.write(f'{a0:.3f}, {c0:.3f}\\n')\n\n# Show the results on screen:\nprint('The optimized lattice constants are:')\nprint(f'a = {a0:.3f} \u00c5, c = {c0:.3f} \u00c5')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the stress tensor\n\nOne can also use the stress tensor to optimize the unit cell.\nUsing 'StrainFilter', the unit cell is relaxed until the stress is\nbelow a given threshold.\n\nNote that if the initial guesses for $a$ and $c$ are far\nfrom the optimal values, the optimization can get stuck in a local minimum.\nSimilarly, if the threshold is not chosen tight enough, the resulting\nlattice constants can again be far from the optimal values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ni = bulk('Ni', 'hcp', a=3.0, c=5.0)\nni.calc = EMT()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "rather then directly to the atoms. The filter presents an\nalternative view of the atomic degrees of freedom: instead of\nmodifying atomic positions to minimise target forces, the [BFGS](https://en.wikipedia.org/wiki/Broyden\u2013Fletcher\u2013Goldfarb\u2013Shanno_algorithm)\nalgorithm is allowed to modify lattice parameters to minimise target\nstress.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sf = StrainFilter(ni)\nopt = BFGS(sf)\n\n# If you want the optimization path in a trajectory, comment in these lines:\n# traj = Trajectory('path.traj', 'w', ni)\n# opt.attach(traj)\n\n# Set the threshold and run the optimization::\nopt.run(0.005) # run until forces are below 0.005 eV/\u00c5\u00b3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Analyze the results\nprint the unit cell\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print('The optimized lattice constants from the stress tensor are:')\nprint(f'a = {ni.cell[0, 0]:.3f} \u00c5, c = {ni.cell[2, 2]:.3f} \u00c5')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.7" } }, "nbformat": 4, "nbformat_minor": 0 }