{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Calculating Delta-values\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this tutorial we compare the equation-of-state (EOS) calculated for 7 FCC\nmetals using values from :class:`~ase.calculators.emt.EMT`, WIEN2k and\nexperiment. Each EOS is described by three parameters:\n\n* volume per atom\n* bulk modulus\n* pressure derivative of the bulk modulus\n\nDifferences between two EOS's can be measured by a single $\\Delta$ value\ndefined as:\n\n\\begin{align}\\Delta = \\sqrt{\\frac{\\int_{V_a}^{V_b}\n      \\left( E_1(V) - E_2(V) \\right)^2 \\, dV}{V_b - V_a}}\\end{align}\n\nwhere $E_n(V)$ is the energy per atom as a function of volume.\nThe $\\Delta$ value can be calculated using the\n:func:`ase.utils.deltacodesdft.delta` function:\n\n\n.. autofunction:: ase.utils.deltacodesdft.delta\n   :noindex:\n\n.. seealso::\n\n   * Collection of ground-state elemental crystals: `dcdft`\n   * Equation-of-state module: :mod:`ase.eos`\n\nWe get the WIEN2k and experimental numbers from\nthe `dcdft` ASE-collection\nand we calculate the EMT EOS using this script:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pathlib import Path\n\nimport numpy as np\n\nfrom ase.calculators.emt import EMT\nfrom ase.collections import dcdft\nfrom ase.io import Trajectory\n\ntrajfiles = []\n\nfor symbol in ['Al', 'Ni', 'Cu', 'Pd', 'Ag', 'Pt', 'Au']:\n    trajfile = Path(f'{symbol}.traj')\n    with Trajectory(trajfile, 'w') as traj:\n        for s in range(94, 108, 2):\n            atoms = dcdft[symbol]\n            atoms.set_cell(atoms.cell * (s / 100) ** (1 / 3), scale_atoms=True)\n            atoms.calc = EMT()\n            atoms.get_potential_energy()\n            traj.write(atoms)\n    trajfiles.append(trajfile)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And fit to a Birch-Murnaghan EOS:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import json\n\nfrom ase.eos import EquationOfState as EOS\nfrom ase.io import read\n\n\ndef fit(symbol: str) -> tuple[float, float, float, float]:\n    V = []\n    E = []\n    for atoms in read(f'{symbol}.traj@:'):\n        V.append(atoms.get_volume() / len(atoms))\n        E.append(atoms.get_potential_energy() / len(atoms))\n\n    eos = EOS(V, E, 'birchmurnaghan')\n    eos.fit(warn=False)\n    e0, B, Bp, v0 = eos.eos_parameters\n    return e0, v0, B, Bp\n\n\ndata = {}  # Dict[str, Dict[str, float]]\nfor path in trajfiles:\n    symbol = path.stem\n    e0, v0, B, Bp = fit(symbol)\n    data[symbol] = {\n        'emt_energy': e0,\n        'emt_volume': v0,\n        'emt_B': B,\n        'emt_Bp': Bp,\n    }\n\nPath('fit.json').write_text(json.dumps(data))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Result for Pt using EMT:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nV, E = [], []\nfor atoms in read('Pt.traj@:'):\n    V.append(atoms.get_volume() / len(atoms))\n    E.append(atoms.get_potential_energy() / len(atoms))\n\neos = EOS(V, E, 'birchmurnaghan')\neos.fit(warn=False)\n\nfig = plt.figure()\nax = fig.gca()\neos.plot(ax=ax)  # draw onto the current axes\nax.set_xlim(14.0, 16.0)\nax.set_xlabel('volume [\u00c5^3/atom]')\nax.set_ylabel('energy [eV/atom]')\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Result for Pt using EMT compared to experiment and WIEN2k\nEquilibrium volumes (\u00c5^3/atom):\n\n.. csv-table::\n   :header: symbol, emt, exp, wien2k\n\n   Pt, 15.08, 15.02, 15.64\n   Al, 15.93, 16.27, 16.48\n   Ni, 10.60, 10.81, 10.89\n   Au, 16.68, 16.82, 17.97\n   Pd, 14.59, 14.56, 15.31\n   Ag, 16.77, 16.85, 17.85\n   Cu, 11.57, 11.65, 11.95\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Bulk moduli in GPa:\n\n.. csv-table::\n   :header: symbol, emt, exp, wien2k\n\n  Pt, 278.67, 285.51, 248.71\n  Al, 39.70, 77.14, 78.08\n  Ni, 176.23, 192.46, 200.37\n  Au, 174.12, 182.01, 139.11\n  Pd, 180.43, 187.19, 168.63\n  Ag, 100.06, 105.71, 90.15\n  Cu, 134.41, 144.28, 141.33\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Pressure derivative of bulk-moduli:\n\n.. csv-table::\n   :header: symbol, emt, exp, wien2k\n\n  Pt, 5.31, 5.18, 5.46\n  Al, 2.72, 4.45, 4.57\n  Ni, 3.76, 4.00, 5.00\n  Au, 5.46, 6.40, 5.76\n  Pd, 5.17, 5.00, 5.56\n  Ag, 4.75, 4.72, 5.42\n  Cu, 4.21, 4.88, 4.86\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we can calculate $\\Delta$ between EMT and WIEN2k for Pt:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.units import kJ\nfrom ase.utils.deltacodesdft import delta\n\ndelta(15.08, 278.67 * 1e-24 * kJ, 5.31, 15.64, 248.71 * 1e-24 * kJ, 5.46)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here are all the `\\Delta` values (in meV/atom) calculated\nwith the script below:\n\n.. csv-table::\n   :header: symbol, emt-exp, emt-wien2k, exp-wien2k\n\n  Pt, 3.5, 32.2, 35.9\n  Al, 5.9, 8.6, 3.6\n  Ni, 8.6, 12.5, 3.7\n  Au, 5.9, 43.7, 39.4\n  Pd, 1.0, 27.6, 29.0\n  Ag, 1.9, 22.4, 21.3\n  Cu, 2.7, 11.9, 9.5\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from ase.eos import birchmurnaghan\n\n# Read EMT data:\ndata = json.loads(Path('fit.json').read_text())\n# Insert values from experiment and WIEN2k:\nfor symbol in data:\n    dcdft_dct = dcdft.data[symbol]\n    dcdft_dct['exp_B'] *= 1e-24 * kJ\n    dcdft_dct['wien2k_B'] *= 1e-24 * kJ\n    data[symbol].update(dcdft_dct)\n\nfor name in ['volume', 'B', 'Bp']:\n    with open(name + '.csv', 'w') as fd:\n        print('# symbol, emt, exp, wien2k', file=fd)\n        for symbol, dct in data.items():\n            values = [\n                dct[code + '_' + name] for code in ['emt', 'exp', 'wien2k']\n            ]\n            if name == 'B':\n                values = [val * 1e24 / kJ for val in values]\n            print(\n                f'{symbol},',\n                ', '.join(f'{value:.2f}' for value in values),\n                file=fd,\n            )\n\nwith open('delta.csv', 'w') as fd:\n    print('# symbol, emt-exp, emt-wien2k, exp-wien2k', file=fd)\n    for symbol, dct in data.items():\n        # Get v0, B, Bp:\n        emt, exp, wien2k = (\n            (dct[code + '_volume'], dct[code + '_B'], dct[code + '_Bp'])\n            for code in ['emt', 'exp', 'wien2k']\n        )\n        print(\n            f'{symbol}, {delta(*emt, *exp) * 1000:.1f}, '\n            f'{delta(*emt, *wien2k) * 1000:.1f}, '\n            f'{delta(*exp, *wien2k) * 1000:.1f}',\n            file=fd,\n        )\n\n        if symbol == 'Pt':\n            va = min(emt[0], exp[0], wien2k[0])\n            vb = max(emt[0], exp[0], wien2k[0])\n            v = np.linspace(0.94 * va, 1.06 * vb)\n            for (v0, B, Bp), code in [\n                (emt, 'EMT'),\n                (exp, 'experiment'),\n                (wien2k, 'WIEN2k'),\n            ]:\n                plt.plot(v, birchmurnaghan(v, 0.0, B, Bp, v0), label=code)\n            e0 = dct['emt_energy']\n            V = []\n            E = []\n            for atoms in read('Pt.traj@:'):\n                V.append(atoms.get_volume() / len(atoms))\n                E.append(atoms.get_potential_energy() / len(atoms) - e0)\n            plt.plot(V, E, 'o')\n            plt.legend()\n            plt.xlabel('volume [Ang^3]')\n            plt.ylabel('energy [eV/atom]')\n            plt.savefig('Pt.png')"
      ]
    }
  ],
  "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
}