{ "cells": [ { "cell_type": "markdown", "id": "1aad7bd0-ad93-4037-90be-39a90e0609fc", "metadata": {}, "source": [ "# Elastic Constants\n", "\n", "Solids respond to small external loads through a reversible elastic response. The strength of the response is characterized by the elastic moduli. {py:mod}`matscipy.elasticity` implements functions for computing elastic moduli from small deformation of atomistic systems that consider potential symmetries of the underlying atomic system, in particular for crystals. {py:mod}`matscipy.elasticity` also implements analytic calculation of elastic moduli for some interatomic potentials, described in more detail below. The computation of elastic moduli is a basic prerequisite for multi-scale modelling of materials, as they are the most basic parameters of continuum material models.\n", "\n", "In this tutorial, we will go over different ways that `matscipy` can compute elastic constants of an atomic configuration.\n", "\n", "## Problem setup\n", "\n", "Let's first create an FCC system and view it interactively:" ] }, { "cell_type": "code", "execution_count": null, "id": "f7edef8b-a391-44fb-8591-2e80db68f9fc", "metadata": {}, "outputs": [], "source": [ "from ase.lattice.cubic import FaceCenteredCubic\n", "\n", "def interactive_view(system):\n", " from ase.visualize import view\n", " # Interactive view of the lattice\n", " v = view(system, viewer='ngl')\n", "\n", " # Resize widget\n", " v.view._remote_call(\"setSize\", target=\"Widget\", args=[\"300px\", \"300px\"])\n", " v.view.center()\n", " return v\n", "\n", "# Define FCC crystal\n", "system = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], size=(3,3,3), symbol='Cu', pbc=(1,1,1))\n", "interactive_view(system)" ] }, { "cell_type": "markdown", "id": "5eb45b2e-ec03-4695-bc4b-3ce832c61964", "metadata": {}, "source": [ "We'll assign a force-field to this system based on the Lennard-Jones potential ({py:class}`LennardJonesCut `) implemented in ``matscipy``." ] }, { "cell_type": "code", "execution_count": null, "id": "3fb38057-8e69-479b-bb17-d5060d05e05c", "metadata": {}, "outputs": [], "source": [ "from matscipy.calculators.pair_potential import PairPotential, LennardJonesCut\n", "from ase.data import reference_states, atomic_numbers\n", "\n", "import numpy as np\n", "\n", "Cu_num = atomic_numbers[\"Cu\"]\n", "lattice = reference_states[Cu_num][\"a\"]\n", "sigma = lattice / (2**(2/3))\n", "\n", "system.calc = PairPotential({\n", " # Type map: define Copper-Copper pair interaction\n", " (Cu_num, Cu_num): LennardJonesCut(epsilon=1, sigma=lattice * 2**(-1/6) / np.sqrt(2), cutoff=3 * lattice)\n", "})" ] }, { "cell_type": "markdown", "id": "d95c366a-9180-426e-944a-1d4745e0f827", "metadata": {}, "source": [ "We test we have a sensible potential energy and that we have equilibrium." ] }, { "cell_type": "code", "execution_count": null, "id": "b9bc1a2c-10a5-49bf-985f-d98d54e77bca", "metadata": {}, "outputs": [], "source": [ "from numpy.testing import assert_allclose\n", "from IPython.display import display, Latex\n", "\n", "display(Latex(f\"$E_\\\\mathrm{{pot}} = {system.get_potential_energy():.1f}$\"))\n", "\n", "# Testing negative potential energy\n", "assert system.get_potential_energy() < 0\n", "\n", "# Testing equilibrium\n", "assert_allclose(system.get_forces(), 0, atol=3e-14)" ] }, { "cell_type": "markdown", "id": "3a3405e6-4efc-49a3-a453-d3d01b73f5a9", "metadata": {}, "source": [ "## Crystalline elastic constants\n", "\n", "Let us first compute elastic constants with the {py:mod}`matscipy.elasticity` module, with two different functions:\n", "\n", "- {py:func}`measure_triclinic_elastic_constants ` computes the full Hooke's tensor by finite differences\n", "- {py:func}`fit_elastic_constants ` computes a range of deformed configurations and fits a Hooke's tensor with least-squares" ] }, { "cell_type": "code", "execution_count": null, "id": "5b96c21b-fe28-4e93-8482-58ec73092b28", "metadata": {}, "outputs": [], "source": [ "from matscipy.elasticity import measure_triclinic_elastic_constants, fit_elastic_constants\n", "from matscipy.elasticity import full_3x3x3x3_to_Voigt_6x6\n", "\n", "C_finite_differences = full_3x3x3x3_to_Voigt_6x6(measure_triclinic_elastic_constants(system))\n", "C_least_squares, _ = fit_elastic_constants(system, verbose=False)" ] }, { "cell_type": "markdown", "id": "12f6c608-7ca3-4782-b054-a18b399c4a8b", "metadata": {}, "source": [ "Let's plot the Hooke tensor:" ] }, { "cell_type": "code", "execution_count": null, "id": "5b2f88bb-e86e-4cb3-a356-985a654b4f6a", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "def spy_constants(ax, constants):\n", " ax.imshow(constants, cmap='RdPu', interpolation='none')\n", " labels = np.full_like(constants, \"\", dtype=object)\n", " labels[:3, :3] = \"$\\\\lambda$\\n\"\n", " labels[(np.arange(3), np.arange(3))] = \"$\\\\lambda + 2\\\\mu$\\n\"\n", " labels[(np.arange(3, 6), np.arange(3, 6))] = \"$\\\\mu$\\n\"\n", "\n", " max_C = constants.max()\n", " \n", " for i in range(constants.shape[0]):\n", " for j in range(constants.shape[1]):\n", " color = \"white\" if constants[i, j] / max_C > 0.7 else \"black\"\n", " numeric = f\"${constants[i, j]:.2f}$\" if np.abs(constants[i, j]) / max_C > 1e-3 else \"$\\\\sim 0$\"\n", "\n", " ax.annotate(labels[i, j] + numeric,\n", " (i, j),\n", " horizontalalignment='center',\n", " verticalalignment='center', color=color)\n", " \n", " ax.set_xticks(np.arange(constants.shape[1]))\n", " ax.set_yticks(np.arange(constants.shape[0]))\n", " ax.set_xticklabels([f\"C$_{{i{j+1}}}$\" for j in range(constants.shape[1])])\n", " ax.set_yticklabels([f\"C$_{{{i+1}j}}$\" for i in range(constants.shape[0])])" ] }, { "cell_type": "code", "execution_count": null, "id": "0cfb76fb-ee48-4b23-96df-4ad78bc85df2", "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(9, 5))\n", "\n", "spy_constants(axs[0], C_finite_differences)\n", "spy_constants(axs[1], C_least_squares)\n", "\n", "axs[0].set_title(f\"Finite differences\")\n", "axs[1].set_title(f\"Least squares\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "11599564-f976-4813-8ebd-408792c9c6df", "metadata": {}, "source": [ "### With second-order derivatives\n", "\n", "Most calculators in ``matscipy`` actually implement second-order derivatives, and can therefore directly compute elastic constants:" ] }, { "cell_type": "code", "execution_count": null, "id": "db51ca03-b229-463d-b6b8-43c19548928d", "metadata": {}, "outputs": [], "source": [ "C = full_3x3x3x3_to_Voigt_6x6(system.calc.get_property(\"elastic_constants\", system))" ] }, { "cell_type": "code", "execution_count": null, "id": "c664eedf-8890-4960-8a59-12a429223eff", "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(9, 5))\n", "\n", "spy_constants(axs[0], C_least_squares)\n", "spy_constants(axs[1], C)\n", "\n", "axs[0].set_title(f\"Least squares\")\n", "axs[1].set_title(f\"Direct evaluation\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "81338c89-9e2b-4146-a2f1-36118a2e5af7", "metadata": {}, "source": [ "## Amorphous elastic constants\n", "\n", "In amorphous solids, non-affine relaxation modes play an important role in elastic deformation.\n", "\n", "### Problem setup\n", "\n", "Let's randomize our atoms to mimic a glassy structure:" ] }, { "cell_type": "code", "execution_count": null, "id": "5643f2cb-c81d-4d7a-bd48-2555492bd101", "metadata": {}, "outputs": [], "source": [ "from ase.io import read\n", "from pathlib import Path\n", "\n", "data_path = Path(\"..\") / \"..\" / \"tests\"\n", "\n", "system = read(data_path / \"CuZr_glass_460_atoms.gz\")\n", "interactive_view(system)" ] }, { "cell_type": "markdown", "id": "872bbcc8-be7a-4433-8d3d-573b45a37b5d", "metadata": {}, "source": [ "Setting up the calculator:" ] }, { "cell_type": "code", "execution_count": null, "id": "b683cc4d-832c-4679-a74d-ae891a873f79", "metadata": {}, "outputs": [], "source": [ "from matscipy.calculators.eam import EAM\n", "\n", "system.calc = EAM(data_path / 'ZrCu.onecolumn.eam.alloy')" ] }, { "cell_type": "code", "execution_count": null, "id": "213de06c-b987-41e9-af09-11998076bc37", "metadata": {}, "outputs": [], "source": [ "display(Latex(f\"$E_\\\\mathrm{{pot}} = {system.get_potential_energy():.1f}$\"))" ] }, { "cell_type": "markdown", "id": "fab46693-0ab7-4fb6-acff-3ee602236587", "metadata": {}, "source": [ "### Fitting elastic constants\n", "\n", "The {py:func}`fit_elastic_constants ` function accepts a minimizer procedure as argument to account for non-affine relaxation modes." ] }, { "cell_type": "code", "execution_count": null, "id": "d5019537-e39d-4ee8-be92-22336b046c11", "metadata": {}, "outputs": [], "source": [ "from ase.optimize import FIRE\n", "\n", "delta = 1e-4 # Configuration change increment\n", "\n", "C_affine, _ = fit_elastic_constants(system, verbose=False, delta=delta)\n", "C_relaxed, _ = fit_elastic_constants(system, verbose=False, delta=delta,\n", " # adjust fmax to desired precision\n", " optimizer=FIRE, fmax=5 * delta)" ] }, { "cell_type": "code", "execution_count": null, "id": "56b55c89-980b-4264-8e55-188e4aa60259", "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 2, figsize=(9, 5))\n", "\n", "spy_constants(axs[0], C_affine)\n", "spy_constants(axs[1], C_relaxed)\n", "\n", "axs[0].set_title(f\"Affine only\")\n", "axs[1].set_title(f\"Affine + Non-affine\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "56b085ef-06ab-464f-8eaf-bceba08fd3ba", "metadata": {}, "source": [ "One can see that elastic constants are significantly reduced when the internal relaxation is included. However, mind that the reduction is *very* dependent on the optimizer's stopping criterion ``fmax``, which should ideally be lower than the deformation increment (we set it higher in the example above for demonstration purposes only)." ] }, { "cell_type": "code", "execution_count": null, "id": "12744ced-a9e2-4adb-aa5f-748ed6cc22b5", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.2" } }, "nbformat": 4, "nbformat_minor": 5 }