diff --git a/docs/source/apidoc/janus_core.rst b/docs/source/apidoc/janus_core.rst index 04a117360..506f31fb3 100644 --- a/docs/source/apidoc/janus_core.rst +++ b/docs/source/apidoc/janus_core.rst @@ -44,6 +44,16 @@ janus\_core.calculations.eos module :undoc-members: :show-inheritance: +janus\_core.calculations.elasticity module +------------------------------------------ + +.. automodule:: janus_core.calculations.elasticity + :members: + :special-members: + :private-members: + :undoc-members: + :show-inheritance: + janus\_core.calculations.geom\_opt module ----------------------------------------- @@ -117,6 +127,17 @@ janus\_core.cli.eos module :undoc-members: :show-inheritance: +janus\_core.cli.elasticity module +--------------------------------- + +.. automodule:: janus_core.cli.elasticity + :members: + :special-members: + :private-members: + :undoc-members: + :show-inheritance: + + janus\_core.cli.geomopt module ------------------------------ diff --git a/docs/source/tutorials/cli/elasticity.ipynb b/docs/source/tutorials/cli/elasticity.ipynb new file mode 100644 index 000000000..4a8a60e0c --- /dev/null +++ b/docs/source/tutorials/cli/elasticity.ipynb @@ -0,0 +1,373 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "view-in-github" + }, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/stfc/janus-core/blob/main/docs/source/tutorials/cli/elasticity.ipynb)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "4vcbRxmIhHLL" + }, + "source": [ + "# Elasticity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set up environment (optional)\n", + "\n", + "These steps are required to run this tutorial with Google Colab. To do so, uncomment and run the cell below.\n", + "\n", + "This will replace pre-installed versions of `numpy` and `torch` in Colab with versions that are known to be compatible with `janus-core`.\n", + "\n", + "It may be possible to skip the steps that uninstall and reinstall `torch`, which will save a considerable amount of time.\n", + "\n", + "These instructions but may work for other systems too, but it is typically preferable to prepare a virtual environment separately before running this notebook if possible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "TF-EiWxyuMc7" + }, + "outputs": [], + "source": [ + "# import locale\n", + "# locale.getpreferredencoding = lambda: \"UTF-8\"\n", + "\n", + "# ! pip uninstall numpy -y # Uninstall pre-installed numpy\n", + "\n", + "# ! pip uninstall torch torchaudio torchvision transformers -y # Uninstall pre-installed torch\n", + "# ! uv pip install torch==2.5.1 # Install pinned version of torch\n", + "\n", + "# ! uv pip install janus-core[mace,visualise] data-tutorials --system # Install janus-core with MACE and WeasWidget, and data-tutorials\n", + "\n", + "# get_ipython().kernel.do_shutdown(restart=True) # Restart kernel to update libraries. This may warn that your session has crashed." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To ensure you have the latest version of `janus-core` installed, compare the output of the following cell to the latest version available at https://pypi.org/project/janus-core/" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from janus_core import __version__\n", + "\n", + "print(__version__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Prepare modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from ase.build import nanotube, bulk\n", + "from ase.lattice.cubic import Diamond\n", + "\n", + "from ase.io import write\n", + "from weas_widget import WeasWidget\n", + "\n", + "Path(\"../data\").mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generation of samples\n", + "\n", + "In ```janus_core``` we can calculate the elasticity tensor via [pymatgen](https://github.com/materialsproject/pymatgen). As an example we will do this for three samples: Aluminium, Diamond, and a Carbon-nanotube.\n", + "\n", + "Using the ASE we can build each sample using the utility functions in ```ase.build```.\n", + "\n", + "Firstly Aluminium can be generated using ```ase.build.bulk```." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "al = bulk(\"Al\", crystalstructure=\"fcc\")*(3,3,3)\n", + "\n", + "write(\"../data/Aluminium.xyz\", al)\n", + "\n", + "v=WeasWidget()\n", + "v.from_ase(al)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we can build Diamond," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diamond = Diamond('C')*(2,2,2)\n", + "\n", + "write(\"../data/Diamond.xyz\", diamond)\n", + "\n", + "v=WeasWidget()\n", + "v.from_ase(diamond)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally the Carbon-nanotube," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nt = nanotube(6, 6, length=4)\n", + "# Place into a box\n", + "nt.cell = [nt.cell[2][2], nt.cell[2][2], nt.cell[2][2]]\n", + "nt.pbc = [True, True, True]\n", + "write(\"../data/Carbon-nanotube.xyz\", nt)\n", + "v=WeasWidget()\n", + "v.from_ase(nt)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Elasticity calculations\n", + "\n", + "```janus-core``` includes the elasticity CLI command for calculating the elasticity tensor. This is done by first creating a set of deformed (strained) structures. The stress tensor tensor ($\\sigma$) is then calculated for each of these deformed structures. Finally the elasticity tensor $C^{ijkl}$ is estimated from the relationship between stress and strain. Or mathematically:\n", + "\n", + "$$ \\sigma^{ij} = C^{ijkl} E^{kl} $$\n", + "\n", + "To find $C^{ijkl}$ we can use the ```janus elasticity``` CLI command As with other ```janus_core``` commands we can via its help message for more information:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! janus elasticity --help" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To calculate the elasticity tensor we can use the following CLI command for our Aluminium sample." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! janus elasticity --arch mace_mp --struct ../data/Aluminium.xyz --no-tracker" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The results will appear ```janus_results/Aluminium-elastic_tensor.dat```, which includes the components of the $3\\times3\\times3\\times3$ in (row-major) Voigt reduced form ($6\\times6$) which is preceded by various derived values such as the shear and bulk moduli." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cat janus_results/Aluminium-elastic_tensor.dat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The experimental bulk, shear, and Young's moduli are approximately 76 GPa, 26 GPa, and 68 GPa for Aluminium. Our results are in the right ball park.\n", + "\n", + "For elasticity the main options for user control are the magnitudes of the applied shear and normal strains as well as their number.\n", + "\n", + "By default 4 shear and 4 normal strains are applied. Split equally positively and negatively, and along each possible direction. This means $2\\times 4 \\times 3 = 24$ total stress calculations.\n", + "\n", + "For example we can increase the shear and normal magnitudes, and apply 16 strains of each type as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! janus elasticity --arch mace_mp --struct ../data/Aluminium.xyz --no-tracker --n-strains 16 --shear-magnitude 0.3 --normal-magnitude 0.2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With the results now being different." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cat janus_results/Aluminium-elastic_tensor.dat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the same way we can simply calculate the data for our Diamond and Carbon nano-tube samples," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! janus elasticity --arch mace_mp --struct ../data/Diamond.xyz --no-tracker" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cat janus_results/Diamond-elastic_tensor.dat" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! janus elasticity --arch mace_mp --struct ../data/Carbon-nanotube.xyz --no-tracker" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cat janus_results/Carbon-nanotube-elastic_tensor.dat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can examine the results to asses the comparative stiffness of the materials, do the results line up with experiments?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, layout='constrained')\n", + "for material, marker in zip(('Aluminium', 'Diamond', 'Carbon-nanotube'), \"o^s\"):\n", + " tensor = np.loadtxt(f\"janus_results/{material}-elastic_tensor.dat\")\n", + " bulk, shear, youngs = tensor[[0, 3, 6]]\n", + "\n", + " theta = [i*2.0*np.pi/3 for i in range(3)]\n", + " ax.scatter(theta, [bulk, shear,youngs], marker=marker, label=material)\n", + " ax.set_xticks(theta, (\"Bulk\", \"Shear\", \"Youngs\"))\n", + "fig.legend(loc='outside upper right')\n", + "fig.suptitle(\"Elastic modulii for our samples\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "accelerator": "GPU", + "colab": { + "authorship_tag": "ABX9TyNvtIsPHVgkt0NvUv51T6ZG", + "gpuType": "T4", + "include_colab_link": true, + "private_outputs": true, + "provenance": [] + }, + "kernelspec": { + "display_name": "janus-core", + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/docs/source/tutorials/python/elasticity.ipynb b/docs/source/tutorials/python/elasticity.ipynb new file mode 100644 index 000000000..0826dc1e5 --- /dev/null +++ b/docs/source/tutorials/python/elasticity.ipynb @@ -0,0 +1,334 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "view-in-github" + }, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/stfc/janus-core/blob/main/docs/source/tutorials/python/elasticity.ipynb)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "4vcbRxmIhHLL" + }, + "source": [ + "# Elasticity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set up environment (optional)\n", + "\n", + "These steps are required to run this tutorial with Google Colab. To do so, uncomment and run the cell below.\n", + "\n", + "This will replace pre-installed versions of `numpy` and `torch` in Colab with versions that are known to be compatible with `janus-core`.\n", + "\n", + "It may be possible to skip the steps that uninstall and reinstall `torch`, which will save a considerable amount of time.\n", + "\n", + "These instructions but may work for other systems too, but it is typically preferable to prepare a virtual environment separately before running this notebook if possible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "TF-EiWxyuMc7" + }, + "outputs": [], + "source": [ + "# import locale\n", + "# locale.getpreferredencoding = lambda: \"UTF-8\"\n", + "\n", + "# ! pip uninstall numpy -y # Uninstall pre-installed numpy\n", + "\n", + "# ! pip uninstall torch torchaudio torchvision transformers -y # Uninstall pre-installed torch\n", + "# ! uv pip install torch==2.5.1 # Install pinned version of torch\n", + "\n", + "# ! uv pip install janus-core[mace,visualise] data-tutorials --system # Install janus-core with MACE and WeasWidget, and data-tutorials\n", + "\n", + "# get_ipython().kernel.do_shutdown(restart=True) # Restart kernel to update libraries. This may warn that your session has crashed." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To ensure you have the latest version of `janus-core` installed, compare the output of the following cell to the latest version available at https://pypi.org/project/janus-core/" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from janus_core import __version__\n", + "()\n", + "print(__version__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Prepare data and modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from weas_widget import WeasWidget\n", + "\n", + "from ase.build import bulk, nanotube\n", + "from ase.lattice.cubic import Diamond\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from janus_core.calculations.elasticity import Elasticity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use `data_tutorials` to get the data required for this tutorial:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generation of samples\n", + "\n", + "In ```janus_core``` we can calculate the elasticity tensor via [pymatgen](https://github.com/materialsproject/pymatgen).\n", + "\n", + "As an example we will do this for three samples: Aluminium, Diamond, and a Carbon-nanotube. Then compare their relative stiffnesses.\n", + "\n", + "Using the ASE we can build each sample using the utility functions in ```ase.build```.\n", + "\n", + "Firstly Aluminium can be generated using ```ase.build.bulk```." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "al = bulk(\"Al\", crystalstructure=\"fcc\")*(3,3,3)\n", + "\n", + "v=WeasWidget()\n", + "v.from_ase(al)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we can build Diamond," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diamond = Diamond('C')*(2,2,2)\n", + "\n", + "v=WeasWidget()\n", + "v.from_ase(diamond)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally the Carbon-nanotube," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nt = nanotube(6, 6, length=4)\n", + "# Place into a box\n", + "nt.cell = [nt.cell[2][2], nt.cell[2][2], nt.cell[2][2]]\n", + "nt.pbc = [True, True, True]\n", + "v=WeasWidget()\n", + "v.from_ase(nt)\n", + "v" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Elasticity calculations\n", + "\n", + "```janus-core``` includes the elasticity calculation type for calculating the elasticity tensor. This is done by first creating a set of deformed (strained) structures. The stress tensor tensor ($\\sigma$) is then calculated for each of these deformed structures. Finally the elasticity tensor $C^{ijkl}$ is estimated from the relationship between stress and strain. Or mathematically:\n", + "\n", + "$$ \\sigma^{ij} = C^{ijkl} E^{kl} $$\n", + "\n", + "To find $C^{ijkl}$ we can use ```janus_core.calculations.elasticity.Elasticity```. First we construct the calculation object and run the calculation using the ```run()``` method. This will first geometry optimize the input structure, and then generate the set of deformed structures from it. Finally stresses are calculated for each deformed structure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "elasticity_aluminium = Elasticity(\n", + " struct=al.copy(),\n", + " arch=\"mace_mp\",\n", + " device=\"cpu\",\n", + " model=\"small\",\n", + " calc_kwargs={\"default_dtype\": \"float64\"}\n", + ").run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The results will be written to ```Al27-elastic_tensor.dat``` (automatically generated from our structures composition). The file includes the components of the $3\\times3\\times3\\times3$ elasticity tensor in (row-major) Voigt reduced form ($6\\times6$), which is preceded by various derived values such as the shear and bulk moduli." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cat janus_results/Al27-elastic_tensor.dat" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The experimental bulk, shear, and Young's moduli are approximately 76 GPa, 26 GPa, and 68 GPa for Aluminium. Our results are in the right ball park.\n", + "\n", + "For elasticity the main options for user control are the magnitudes of the applied shear and normal strains as well as their number.\n", + "\n", + "By default 4 shear and 4 normal strains are applied. Split equally positively and negatively, and along each possible direction. This means $2\\times 4 \\times 3 = 24$ total stress calculations.\n", + "\n", + "For example we can increase the shear and normal magnitudes, and apply 16 strains of each type as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "elasticity_aluminium = Elasticity(\n", + " struct=al.copy(),\n", + " arch=\"mace_mp\",\n", + " device=\"cpu\",\n", + " model=\"small\",\n", + " calc_kwargs={\"default_dtype\": \"float64\"},\n", + " shear_magnitude=0.3,\n", + " normal_magnitude=0.2,\n", + " n_strains=16\n", + ").run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "elasticity_diamond = Elasticity(\n", + " struct=diamond.copy(),\n", + " arch=\"mace_mp\",\n", + " device=\"cpu\",\n", + " model=\"small\",\n", + " calc_kwargs={\"default_dtype\": \"float64\"}\n", + ").run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "elasticity_nt = Elasticity(\n", + " struct=nt.copy(),\n", + " arch=\"mace_mp\",\n", + " device=\"cpu\",\n", + " model=\"small\",\n", + " calc_kwargs={\"default_dtype\": \"float64\"}\n", + ").run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, layout='constrained')\n", + "\n", + "for tensor, material, marker in zip((elasticity_aluminium, elasticity_diamond, elasticity_nt), (\"Aluminium\", \"Diamond\", \"Carbon-nanotube\"), \"o^s\"):\n", + " ecs = dict()\n", + " ecs[\"Bulk modulus\"] = tensor.k_voigt\n", + " ecs[\"Shear modulus\"] = tensor.g_voigt\n", + " # Note https://github.com/materialsproject/pymatgen/issues/4435\n", + " ecs[\"Young's modulus\"] = tensor.y_mod/1e9\n", + " theta = [i*2.0*np.pi/(len(ecs)) for i in range(len(ecs))]\n", + " ax.scatter(theta, ecs.values(), marker=marker, label=material)\n", + " ax.set_xticks(theta, ecs.keys())\n", + "\n", + "fig.legend(loc='outside upper right')\n", + "fig.suptitle(\"Elastic modulii for our samples\")" + ] + } + ], + "metadata": { + "accelerator": "GPU", + "colab": { + "authorship_tag": "ABX9TyNvtIsPHVgkt0NvUv51T6ZG", + "gpuType": "T4", + "include_colab_link": true, + "private_outputs": true, + "provenance": [] + }, + "kernelspec": { + "display_name": "janus-core", + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/docs/source/user_guide/command_line.rst b/docs/source/user_guide/command_line.rst index 2644d1660..4399daa22 100644 --- a/docs/source/user_guide/command_line.rst +++ b/docs/source/user_guide/command_line.rst @@ -71,6 +71,7 @@ This should return something similar to: Commands: descriptors Calculate MLIP descriptors. eos Calculate equation of state. + elasticity Calculate elasticity tensors. geomopt Perform geometry optimization and save optimized structure... md Run molecular dynamics simulation, and save trajectory and... neb Run Nudged Elastic Band method. @@ -480,6 +481,32 @@ Optimization at constant volume for all generated structures can also be perform For all options, run ``janus eos --help``. +Elasticity +---------- + +Calculate the elasticity tensor for a given structure (using the `MACE-MP `_ "small" force-field): + +.. code-block:: bash + + janus elasticity --struct tests/data/NaCl.cif --arch mace_mp + +This will use `pymatgen `_ to generate a set of deformed structures of the input structure. These +are used to estimate the stress-strain relationship and thereby calculate the elasticity tensor. + +The output, by default, is written to ```janus_results/NaCl-elastic_tensor.dat``` containing derived values of the bulk and shear modulus in +Reuss/Voigt/VRH schemes, Young's modulus, the Universal anisotropy, Homogeneous Poisson ratio, and finally the elasticty tensor in row-major +(Voigt form by default). The units are GPa throughout. + +By default 4 shear and 4 normal strains are applied. This means two positively and two negatively sheared structures in each possible direction. +Which results in 4*3*2 stress calculations. + +These can be adjusted in number and independently in magnitude. For example to use 10 shear and normal strains with maximum magnitudes 0.1 +and 0.02 the following command can be executed: + +.. code-block:: bash + + janus elasticity --struct tests/data/NaCl.cif --arch mace_mp --shear-magnitude 0.1 --normal-magnitude 0.02 --n-strains 10 + Phonons ------- diff --git a/janus_core/calculations/elasticity.py b/janus_core/calculations/elasticity.py new file mode 100644 index 000000000..5be345ece --- /dev/null +++ b/janus_core/calculations/elasticity.py @@ -0,0 +1,404 @@ +"""Elasticity Tensor.""" + +from __future__ import annotations + +from copy import copy +from typing import Any + +from ase import Atoms +from ase.units import GPa +import numpy as np +from pymatgen.analysis.elasticity import ElasticTensor +from pymatgen.analysis.elasticity.strain import ( + DeformedStructureSet, +) +from pymatgen.analysis.elasticity.stress import Stress +from pymatgen.io.ase import AseAtomsAdaptor + +from janus_core.calculations.base import BaseCalculation +from janus_core.calculations.geom_opt import GeomOpt +from janus_core.helpers.janus_types import ( + Architectures, + ASEReadArgs, + Devices, + OutputKwargs, + PathLike, +) +from janus_core.helpers.struct_io import output_structs +from janus_core.helpers.utils import build_file_dir, none_to_dict, set_minimize_logging + + +class Elasticity(BaseCalculation): + """ + Calculate the elasticity tensor. + + Parameters + ---------- + struct + ASE Atoms structure(s), or filepath to structure(s) to simulate. + arch + MLIP architecture to use for calculations. Default is `None`. + device + Device to run model on. Default is "cpu". + model + MLIP model label, path to model, or loaded model. Default is `None`. + model_path + Deprecated. Please use `model`. + read_kwargs + Keyword arguments to pass to ase.io.read. Default is {}. + calc_kwargs + Keyword arguments to pass to the selected calculator. Default is {}. + attach_logger + Whether to attach a logger. Default is True if "filename" is passed in + log_kwargs, else False. + log_kwargs + Keyword arguments to pass to `config_logger`. Default is {}. + track_carbon + Whether to track carbon emissions of calculation. Requires attach_logger. + Default is True if attach_logger is True, else False. + tracker_kwargs + Keyword arguments to pass to `config_tracker`. Default is {}. + file_prefix + Prefix for output filenames. Default is `None`. + minimize + Whether to optimize geometry for the initial structure. + Default is True. + minimize_all + Whether to optimize geometry for all generated structures. + Default is False. + minimize_kwargs + Keyword arguments to pass to optimize. Default is `None`. + write_results + Whether to write out the elasticity tensor. Default is True. + write_structures + Whether to write out all generated structures. Deault is False. + write_kwargs + Keyword arguments to pass to ase.io.write to save generated structures. + Default is {}. + write_voigt + Whether to write out in Voigt notation, Default is True. + shear_magnitude + The magnitude of shear strain to apply for deformed structures. + Default is 0.06. + normal_magnitude + The magnitude of normal strain to apply for deformed structures. + Default is 0.01. + n_strains + The number of normal and shear strains to apply for deformed structures. + Default is 4. + + Attributes + ---------- + elastic_tensor: ElasticTensor + The ElasticTensor (GPa). + """ + + def __init__( + self, + struct: Atoms | PathLike, + arch: Architectures | None = None, + device: Devices = "cpu", + model: PathLike | None = None, + model_path: PathLike | None = None, + read_kwargs: ASEReadArgs | None = None, + calc_kwargs: dict[str, Any] | None = None, + attach_logger: bool | None = None, + log_kwargs: dict[str, Any] | None = None, + track_carbon: bool | None = None, + tracker_kwargs: dict[str, Any] | None = None, + file_prefix: PathLike | None = None, + minimize: bool = True, + minimize_all: bool = False, + minimize_kwargs: dict[str, Any] | None = None, + write_results: bool = True, + write_structures: bool = False, + write_kwargs: OutputKwargs | None = None, + write_voigt: bool = True, + shear_magnitude: float = 0.06, + normal_magnitude: float = 0.01, + n_strains: int = 4, + ) -> None: + """ + Initialise class. + + Parameters + ---------- + struct + ASE Atoms structure(s), or filepath to structure(s) to simulate. + arch + MLIP architecture to use for calculations. Default is `None`. + device + Device to run model on. Default is "cpu". + model + MLIP model label, path to model, or loaded model. Default is `None`. + model_path + Deprecated. Please use `model`. + read_kwargs + Keyword arguments to pass to ase.io.read. Default is {}. + calc_kwargs + Keyword arguments to pass to the selected calculator. Default is {}. + attach_logger + Whether to attach a logger. Default is True if "filename" is passed in + log_kwargs, else False. + log_kwargs + Keyword arguments to pass to `config_logger`. Default is {}. + track_carbon + Whether to track carbon emissions of calculation. Requires attach_logger. + Default is True if attach_logger is True, else False. + tracker_kwargs + Keyword arguments to pass to `config_tracker`. Default is {}. + file_prefix + Prefix for output filenames. Default is `None`. + minimize + Whether to optimize geometry for the initial structure. + Default is True. + minimize_all + Whether to optimize geometry for all generated structures. + Default is False. + minimize_kwargs + Keyword arguments to pass to optimize. Default is `None`. + write_results + Whether to write out the elasticity tensor. Default is True. + write_structures + Whether to write out all generated structures. Deault is False. + write_kwargs + Keyword arguments to pass to ase.io.write to save generated structures. + Default is {}. + write_voigt + Whether to write out in Voigt notation, Default is True. + shear_magnitude + The magnitude of shear strain to apply for deformed structures. + Default is 0.06. + normal_magnitude + The magnitude of normal strain to apply for deformed structures. + Default is 0.01. + n_strains + The number of normal and shear strains to apply for deformed structures. + Default is 4. + """ + read_kwargs, minimize_kwargs, write_kwargs = none_to_dict( + read_kwargs, minimize_kwargs, write_kwargs + ) + + self.minimize = minimize + self.minimize_all = minimize_all + self.minimize_kwargs = minimize_kwargs + self.write_results = write_results + self.write_structures = write_structures + self.write_kwargs = write_kwargs + self.write_voigt = write_voigt + + strains = np.array( + [-1 + 2 * i / n_strains for i in range(n_strains // 2)] + + [2 * (i + 1) / n_strains for i in range(n_strains // 2)] + ) + self.normal_strains = strains * normal_magnitude + self.shear_strains = strains * shear_magnitude + + if ( + (self.minimize or self.minimize_all) + and "write_results" in self.minimize_kwargs + and self.minimize_kwargs["write_results"] + ): + raise ValueError( + "Please set the `write_structures` parameter to `True` to save " + "optimized structures, instead of passing `write_results` through " + "`minimize_kwargs`" + ) + + # Read last image by default + read_kwargs.setdefault("index", -1) + + # Initialise structures and logging + super().__init__( + struct=struct, + calc_name=__name__, + arch=arch, + device=device, + model=model, + model_path=model_path, + read_kwargs=read_kwargs, + sequence_allowed=False, + calc_kwargs=calc_kwargs, + attach_logger=attach_logger, + log_kwargs=log_kwargs, + track_carbon=track_carbon, + tracker_kwargs=tracker_kwargs, + file_prefix=file_prefix, + ) + + if not self.struct.calc: + raise ValueError("Please attach a calculator to `struct`.") + + set_minimize_logging( + self.logger, self.minimize_kwargs, self.log_kwargs, track_carbon + ) + + # Set output files + self.write_kwargs["filename"] = self._build_filename( + "generated.extxyz", filename=self.write_kwargs.get("filename") + ) + + self.elasticity_file = self._build_filename("elastic_tensor.dat") + + self.elastic_tensor = None + + @property + def output_files(self) -> None: + """ + Dictionary of output file labels and paths. + + Returns + ------- + dict[str, PathLike] + Output file labels and paths. + """ + return { + "log": self.log_kwargs["filename"] if self.logger else None, + "generated_structures": self.write_kwargs["filename"] + if self.write_structures + else None, + "elasticity": self.elasticity_file if self.write_results else None, + } + + def run(self) -> ElasticTensor: + """ + Calculate the elasticity tensor. + + Returns + ------- + ElasticTensor + The ElasticTensor (GPa). + """ + self._set_info_units() + + if self.minimize: + if self.logger: + self.logger.info("Minimising initial structure") + optimizer = GeomOpt(self.struct, **self.minimize_kwargs) + optimizer.run() + + # Optionally write structure to file + output_structs( + images=self.struct, + struct_path=self.struct_path, + write_results=self.write_structures, + write_kwargs=self.write_kwargs, + config_type="elasticity", + ) + + self._calculate_elasticity() + + if self.write_results: + build_file_dir(self.elasticity_file) + with open(self.elasticity_file, "w", encoding="utf8") as out: + print( + "# Bulk modulus (Reuss) [GPa] |" + " Bulk modulus (Voigt) [GPa] |" + " Bulk modulus (VRH) [GPa] |" + " Shear modulus (Reuss) [GPa] |" + " Shear modulus (Voigt) [GPa] |" + " Shear modulus (VRH) [GPa] |" + " Young's modulus [GPa] |" + " Universal anisotropy |" + " Homogeneous Poisson ratio |" + " Elastic constants (row-major) [GPa]", + file=out, + ) + + values = [ + self.elastic_tensor.property_dict[prop] + for prop in ( + "k_reuss", + "k_voigt", + "k_vrh", + "g_reuss", + "g_voigt", + "g_vrh", + ) + ] + + # y_mod multiplies by 1e9, which happens to convert + # to Pa if the tensor is actually supplied in GPa. + # See https://github.com/materialsproject/pymatgen/issues/4435 + values.append(self.elastic_tensor.y_mod / 1e9) + + for dimensionless in ("universal_anisotropy", "homogeneous_poisson"): + values.append(self.elastic_tensor.property_dict[dimensionless]) + + vals = ( + self.elastic_tensor.voigt + if self.write_voigt + else self.elastic_tensor + ) + for cijkl in vals.flatten(): + values.append(cijkl) + + print(" ".join(map(str, values)), file=out) + + return self.elastic_tensor + + def _calculate_elasticity(self) -> None: + """Generate deformed structures and calculate the elasticity.""" + if self.logger: + self.logger.info("Starting structure calculations") + if self.tracker: + self.tracker.start_task("Calculate configurations") + + self.deformed_structure_set = DeformedStructureSet( + AseAtomsAdaptor.get_structure(self.struct), + norm_strains=self.normal_strains, + shear_strains=self.shear_strains, + ) + + self.stresses = [] + self.strains = [ + i.green_lagrange_strain for i in self.deformed_structure_set.deformations + ] + self.deformed_structures = [ + AseAtomsAdaptor.get_atoms(struct) for struct in self.deformed_structure_set + ] + for i, deformed_structure in enumerate(self.deformed_structures): + deformed_structure.calc = copy(self.struct.calc) + progress = f"{i} / {len(self.deformed_structures)}" + if self.minimize_all: + if self.logger: + self.logger.info( + "Calculating stress for deformed structure " + progress + ) + optimizer = GeomOpt(deformed_structure, **self.minimize_kwargs) + optimizer.run() + + if self.logger: + self.logger.info( + "Calculating stress for deformed structure " + progress, + ) + + # Always append first original structure + self.write_kwargs["append"] = True + # Write structures, but no need to set info c_struct is not used elsewhere + output_structs( + images=deformed_structure, + struct_path=self.struct_path, + write_results=self.write_structures, + set_info=False, + write_kwargs=self.write_kwargs, + config_type="elasticity", + ) + + self.stresses.append(deformed_structure.get_stress(voigt=False) / GPa) + + if self.logger: + self.logger.info("Calculating stress for initial structure") + + self.eq_stress = self.struct.get_stress(voigt=False) / GPa + + self.elastic_tensor = ElasticTensor.from_independent_strains( + self.strains, [Stress(s) for s in self.stresses], self.eq_stress + ) + + if self.logger: + self.logger.info("Structure calculations complete.") + if self.tracker: + emissions = self.tracker.stop_task().emissions + self.struct.info["emissions"] = emissions diff --git a/janus_core/cli/elasticity.py b/janus_core/cli/elasticity.py new file mode 100644 index 000000000..bd22519ae --- /dev/null +++ b/janus_core/cli/elasticity.py @@ -0,0 +1,268 @@ +"""Set up elasticity commandline interface.""" + +from __future__ import annotations + +from typing import Annotated + +from typer import Context, Option, Typer +from typer_config import use_config + +from janus_core.cli.types import ( + Architecture, + CalcKwargs, + Device, + FilePrefix, + LogPath, + MinimizeKwargs, + Model, + ModelPath, + ReadKwargsLast, + StructPath, + Summary, + Tracker, + WriteKwargs, +) +from janus_core.cli.utils import yaml_converter_callback + +app = Typer() + + +@app.command() +@use_config(yaml_converter_callback, param_help="Path to configuration file.") +def elasticity( + # numpydoc ignore=PR02 + ctx: Context, + # Required + arch: Architecture, + struct: StructPath, + # Calculation + shear_magnitude: Annotated[ + float, + Option( + help="Magnitude of shear strain for deformed structures.", + rich_help_panel="Calculation", + ), + ] = 0.06, + normal_magnitude: Annotated[ + float, + Option( + help="Magnitude of normal strain for deformed structures.", + rich_help_panel="Calculation", + ), + ] = 0.01, + n_strains: Annotated[ + int, + Option( + help="Number of normal and shear strains to use for deformed structures.", + rich_help_panel="Calculation", + ), + ] = 4, + write_voigt: Annotated[ + bool, + Option( + help="Write the ElasticityTensor in Voigt notation.", + rich_help_panel="Calculation", + ), + ] = True, + minimize: Annotated[ + bool, + Option( + help="Whether to minimize initial structure before calculations.", + rich_help_panel="Calculation", + ), + ] = True, + minimize_all: Annotated[ + bool, + Option( + help="Whether to minimize all generated structures for calculations.", + rich_help_panel="Calculation", + ), + ] = False, + fmax: Annotated[ + float, + Option( + help="Maximum force for optimization convergence.", + rich_help_panel="Calculation", + ), + ] = 0.1, + minimize_kwargs: MinimizeKwargs = None, + write_structures: Annotated[ + bool, + Option( + help="Whether to write out all genereated structures.", + rich_help_panel="Calculation", + ), + ] = False, + # MLIP Calculator + device: Device = "cpu", + model: Model = None, + model_path: ModelPath = None, + calc_kwargs: CalcKwargs = None, + # Structure I/O + file_prefix: FilePrefix = None, + read_kwargs: ReadKwargsLast = None, + write_kwargs: WriteKwargs = None, + # Logging/summary + log: LogPath = None, + tracker: Tracker = True, + summary: Summary = None, +) -> None: + """ + Calculate equation of state and write out results. + + Parameters + ---------- + ctx + Typer (Click) Context. Automatically set. + arch + MLIP architecture to use for calculations. + struct + Path of structure to simulate. + shear_magnitude + The magnitude of shear strain to apply for deformed structures. + Default is 0.06. + normal_magnitude + The magnitude of normal strain to apply for deformed structures. + Default is 0.01. + n_strains + The number of normal and shear strains to apply for deformed structures. + Default is 4. + write_voigt + Whether to write out in Voigt notation, Default is True. + minimize + Whether to minimize initial structure before calculations. Default is True. + minimize_all + Whether to optimize geometry for all generated structures. Default is False. + fmax + Set force convergence criteria for optimizer in units eV/Å. + Default is 0.1. + minimize_kwargs + Other keyword arguments to pass to geometry optimizer. Default is {}. + write_structures + True to write out all genereated structures. Default is False. + device + Device to run model on. Default is "cpu". + model + Path to MLIP model or name of model. Default is `None`. + model_path + Deprecated. Please use `model`. + calc_kwargs + Keyword arguments to pass to the selected calculator. Default is {}. + file_prefix + Prefix for output files, including directories. Default directory is + ./janus_results, and default filename prefix is inferred from the input + stucture filename. + read_kwargs + Keyword arguments to pass to ase.io.read. By default, + read_kwargs["index"] is -1. + write_kwargs + Keyword arguments to pass to ase.io.write to save generated structures. + Default is {}. + log + Path to write logs to. Default is inferred from `file_prefix`. + tracker + Whether to save carbon emissions of calculation in log file and summary. + Default is True. + summary + Path to save summary of inputs, start/end time, and carbon emissions. Default + is inferred from `file_prefix`. + config + Path to yaml configuration file to define the above options. Default is None. + """ + from janus_core.calculations.elasticity import Elasticity + from janus_core.cli.utils import ( + carbon_summary, + check_config, + end_summary, + get_config, + get_struct_info, + parse_typer_dicts, + set_read_kwargs_index, + start_summary, + ) + + # Check options from configuration file are all valid + check_config(ctx) + + [read_kwargs, calc_kwargs, minimize_kwargs, write_kwargs] = parse_typer_dicts( + [read_kwargs, calc_kwargs, minimize_kwargs, write_kwargs] + ) + + # Set initial config + all_kwargs = { + "read_kwargs": read_kwargs.copy(), + "calc_kwargs": calc_kwargs.copy(), + "minimize_kwargs": minimize_kwargs.copy(), + "write_kwargs": write_kwargs.copy(), + } + config = get_config(params=ctx.params, all_kwargs=all_kwargs) + + # Read only first structure by default and ensure only one image is read + set_read_kwargs_index(read_kwargs) + + # Check fmax option not duplicated + if "fmax" in minimize_kwargs: + raise ValueError("'fmax' must be passed through the --fmax option") + minimize_kwargs["fmax"] = fmax + + log_kwargs = {"filemode": "w"} + if log: + log_kwargs["filename"] = log + + # Dictionary of inputs for Elasticity class + elasticity_kwargs = { + "struct": struct, + "arch": arch, + "device": device, + "model": model, + "model_path": model_path, + "read_kwargs": read_kwargs, + "calc_kwargs": calc_kwargs, + "attach_logger": True, + "log_kwargs": log_kwargs, + "track_carbon": tracker, + "minimize": minimize, + "minimize_all": minimize_all, + "minimize_kwargs": minimize_kwargs, + "write_structures": write_structures, + "write_kwargs": write_kwargs, + "file_prefix": file_prefix, + "write_voigt": write_voigt, + "shear_magnitude": shear_magnitude, + "normal_magnitude": normal_magnitude, + "n_strains": n_strains, + } + + # Initialise Elasticity + elasticity = Elasticity(**elasticity_kwargs) + + # Set summary and log files + summary = elasticity._build_filename("elasticity-summary.yml", filename=summary) + log = elasticity.log_kwargs["filename"] + + # Add structure, MLIP information, and log to info + info = get_struct_info( + struct=elasticity.struct, + struct_path=struct, + ) + + output_files = elasticity.output_files + + # Save summary information before calculations begin + start_summary( + command="elasticity", + summary=summary, + info=info, + config=config, + output_files=output_files, + ) + + # Run elasticity calculations + elasticity.run() + + # Save carbon summary + if tracker: + carbon_summary(summary=summary, log=log) + + # Time after calculations have finished + end_summary(summary) diff --git a/janus_core/cli/janus.py b/janus_core/cli/janus.py index 76b551942..a5278628f 100644 --- a/janus_core/cli/janus.py +++ b/janus_core/cli/janus.py @@ -8,6 +8,7 @@ from janus_core import __version__ from janus_core.cli.descriptors import descriptors +from janus_core.cli.elasticity import elasticity from janus_core.cli.eos import eos from janus_core.cli.geomopt import geomopt from janus_core.cli.md import md @@ -42,6 +43,10 @@ help="Calculate equation of state.", rich_help_panel="Calculations", )(eos) +app.command( + help="Calculate elasticity tensors.", + rich_help_panel="Calculations", +)(elasticity) app.command( help="Run Nudged Elastic Band method.", rich_help_panel="Calculations", diff --git a/janus_core/cli/types.py b/janus_core/cli/types.py index 44cb9a0f3..17055fd1b 100644 --- a/janus_core/cli/types.py +++ b/janus_core/cli/types.py @@ -353,6 +353,22 @@ def __str__(self) -> str: ), ] +ElasticityKwargs = Annotated[ + TyperDict | None, + Option( + parser=parse_dict_class, + help=( + """ + Keyword arguments to pass to elasiticty for user strains. Must be + passed as a dictionary wrapped in quotes, e.g. "[{'key' : values}]". + """ + ), + metavar="DICT", + rich_help_panel="Calculation", + ), +] + + LogPath = Annotated[ Path | None, Option( diff --git a/tests/test_elasticity.py b/tests/test_elasticity.py new file mode 100644 index 000000000..2f64199e1 --- /dev/null +++ b/tests/test_elasticity.py @@ -0,0 +1,123 @@ +"""Test elasticity calculations.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.build import bulk +import numpy as np +from pytest import approx + +from janus_core.calculations.elasticity import Elasticity +from janus_core.helpers.mlip_calculators import choose_calculator +from tests.utils import assert_log_contains + +DATA_PATH = Path(__file__).parent / "data" +MODEL_PATH = Path(__file__).parent / "models" / "mace_mp_small.model" + + +def test_calc_elasticity(tmp_path): + """Test calculating elasticity for Aluminium.""" + elasticity_path = tmp_path / "elasticity-elastic_tensor.dat" + log_file = tmp_path / "elasticity.log" + + struct = bulk("Al", crystalstructure="fcc") + struct.calc = choose_calculator(arch="mace_mp", model=MODEL_PATH) + + elasticity = Elasticity( + struct, + file_prefix=tmp_path / "elasticity", + log_kwargs={"filename": log_file}, + shear_magnitude=0.4, + normal_magnitude=0.2, + n_strains=8, + ) + + elastic_tensor = elasticity.run() + + assert len(elasticity.deformed_structure_set) == 48 + assert elastic_tensor.k_reuss == approx(97.07446854941564, rel=1e-3) + + # Check geometry optimization run by default + assert_log_contains( + log_file, + includes=["Using filter", "Using optimizer", "Starting geometry optimization"], + ) + + assert elasticity_path.exists() + written_elasticity = np.loadtxt(elasticity_path) + # Defaulting to 6x6 Voigt notation + assert len(written_elasticity) == 45 + + for i, prop in enumerate( + ( + "k_reuss", + "k_voigt", + "k_vrh", + "g_reuss", + "g_voigt", + "g_vrh", + "y_mod", + "universal_anisotropy", + "homogeneous_poisson", + ) + ): + units = 1.0 / 1e9 if prop == "y_mod" else 1.0 + assert elastic_tensor.property_dict[prop] * units == approx( + written_elasticity[i], rel=1e-3 + ) + + assert written_elasticity[9:] == approx(elastic_tensor.voigt.flatten(), rel=1e-3) + + +def test_no_optimize_no_write_voigt(tmp_path): + """Test calculating elasticity for Aluminium without optimization.""" + elasticity_path = tmp_path / "elasticity-elastic_tensor.dat" + log_file = tmp_path / "elasticity.log" + + struct = bulk("Al", crystalstructure="fcc") + struct.calc = choose_calculator(arch="mace_mp", model=MODEL_PATH) + + elasticity = Elasticity( + struct, + file_prefix=tmp_path / "elasticity", + log_kwargs={"filename": log_file}, + minimize=False, + write_voigt=False, + ) + elastic_tensor = elasticity.run() + + assert elastic_tensor.k_reuss == approx(80.34637735135381, rel=1e-3) + + # Check geometry optimization run by default + assert_log_contains( + log_file, + excludes=["Using filter", "Using optimizer", "Starting geometry optimization"], + ) + + assert elasticity_path.exists() + written_elasticity = np.loadtxt(elasticity_path) + # Selected to full tensor 3x3x3x3 + assert len(written_elasticity) == 90 + + assert written_elasticity[0] == approx(elastic_tensor.k_reuss, rel=1e-3) + + for i, prop in enumerate( + ( + "k_reuss", + "k_voigt", + "k_vrh", + "g_reuss", + "g_voigt", + "g_vrh", + "y_mod", + "universal_anisotropy", + "homogeneous_poisson", + ) + ): + units = 1.0 / 1e9 if prop == "y_mod" else 1.0 + assert elastic_tensor.property_dict[prop] * units == approx( + written_elasticity[i], rel=1e-3 + ) + + assert written_elasticity[9:] == approx(elastic_tensor.flatten(), rel=1e-3) diff --git a/tests/test_elasticity_cli.py b/tests/test_elasticity_cli.py new file mode 100644 index 000000000..6e1f9a5b0 --- /dev/null +++ b/tests/test_elasticity_cli.py @@ -0,0 +1,91 @@ +"""Test elasticity commandline interface.""" + +from __future__ import annotations + +from pathlib import Path + +from ase.io import read +import numpy as np +from pytest import approx +from typer.testing import CliRunner +import yaml + +from janus_core.cli.janus import app +from tests.utils import assert_log_contains, chdir, strip_ansi_codes + +DATA_PATH = Path(__file__).parent / "data" +MACE_PATH = Path(__file__).parent / "models" / "mace_mp_small.model" + +runner = CliRunner() + + +def test_help(): + """Test calling janus elasticity --help.""" + result = runner.invoke(app, ["elasticity", "--help"]) + assert result.exit_code == 0 + assert "Usage: janus elasticity [OPTIONS]" in strip_ansi_codes(result.stdout) + + +def test_elasticity(tmp_path): + """Test calculating the ElasticTensor from the command line.""" + with chdir(tmp_path): + results_dir = Path("./janus_results") + elasticity_path = results_dir / "NaCl-elastic_tensor.dat" + log_path = results_dir / "NaCl-elasticity-log.yml" + summary_path = results_dir / "NaCl-elasticity-summary.yml" + generated_path = results_dir / "NaCl-generated.extxyz" + + result = runner.invoke( + app, + [ + "elasticity", + "--struct", + DATA_PATH / "NaCl.cif", + "--arch", + "mace_mp", + "--n-strains", + "2", + "--minimize-all", + "--write-structures", + ], + ) + + assert result.exit_code == 0 + + assert elasticity_path.exists() + assert log_path.exists() + assert summary_path.exists() + assert generated_path.exists() + + written_elasticity = np.loadtxt(elasticity_path) + + assert written_elasticity[0] == approx(27.368617328271498) + assert written_elasticity[-1] == approx(2.1775175649257585) + + generated = read(generated_path, index=":") + assert len(generated) == 13 + + assert_log_contains(log_path, includes=["Minimising initial structure"]) + + with open(log_path) as io: + logs = yaml.safe_load(io.read()) + optimization_logs = sum( + log["message"] == ["Starting geometry optimization"] for log in logs + ) + assert optimization_logs == 13 + + with open(summary_path, encoding="utf8") as file: + elasticity_summary = yaml.safe_load(file) + + assert "command" in elasticity_summary + assert "janus elasticity" in elasticity_summary["command"] + assert "start_time" in elasticity_summary + assert "config" in elasticity_summary + assert "info" in elasticity_summary + assert "end_time" in elasticity_summary + + assert "emissions" in elasticity_summary + assert elasticity_summary["emissions"] > 0 + + assert "n_strains" in elasticity_summary["config"] + assert elasticity_summary["config"]["n_strains"] == 2