Elasticity¶
Set up environment (optional)¶
These steps are required to run this tutorial with Google Colab. To do so, uncomment and run the cell below.
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.
[1]:
# ! uv pip install git+https://github.com/WillBaldwin0/graph_electrostatics@v0.4.0
# ! uv pip install janus-core[mace,orb,visualise] data-tutorials --system # Install janus-core with MACE, Orb, and WeasWidget, and data-tutorials
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/
[2]:
from janus_core import __version__
print(__version__)
0.9.9
Prepare modules¶
[3]:
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from ase.build import nanotube, bulk
from ase.lattice.cubic import Diamond
from ase.io import write, read
from weas_widget import WeasWidget
Path("../data").mkdir(exist_ok=True)
Generation of samples¶
In janus_core we can calculate the elasticity tensor via pymatgen. As an example we will do this for three samples: Aluminium, Diamond, and a Carbon-nanotube.
Using the ASE we can build each sample using the utility functions in ase.build.
Firstly Aluminium can be generated using ase.build.bulk.
[4]:
al = bulk("Al", crystalstructure="fcc")*(3,3,3)
write("../data/Aluminium.xyz", al)
v=WeasWidget()
v.from_ase(al)
v
[4]:
Next we can build Diamond,
[5]:
diamond = Diamond('C')*(2,2,2)
write("../data/Diamond.xyz", diamond)
v=WeasWidget()
v.from_ase(diamond)
v
[5]:
And finally the Carbon-nanotube,
[6]:
nt = nanotube(6, 6, length=4)
# Place into a box
nt.cell = [nt.cell[2][2], nt.cell[2][2], nt.cell[2][2]]
nt.pbc = [True, True, True]
write("../data/Carbon-nanotube.xyz", nt)
v=WeasWidget()
v.from_ase(nt)
v
[6]:
Elasticity calculations¶
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:
\[\sigma^{ij} = C^{ijkl} E^{kl}\]
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:
[7]:
! janus elasticity --help
Usage: janus elasticity [OPTIONS]
Calculate elasticity tensors.
╭─ Options ────────────────────────────────────────────────────────────────────╮
│ --config <str> Path to configuration file. │
│ --help Show this message and exit. │
╰──────────────────────────────────────────────────────────────────────────────╯
╭─ MLIP calculator ────────────────────────────────────────────────────────────╮
│ * --arch <mace|mace_mp|mace_off|c MLIP architecture to use │
│ hgnet|sevennet|nequip|de for calculations. │
│ epmd|orb|mattersim|grace [required] │
│ |upet|fairchem|mace_omol │
│ |mace_polar|dpa3> │
│ --device <cpu|cuda|mps|xpu> Device to run calculations │
│ on. │
│ [default: cpu] │
│ --model <str> MLIP model name, or path │
│ to model. │
│ --calc-kwargs DICT Keyword arguments to pass │
│ to selected calculator. │
│ Must be passed as a │
│ dictionary wrapped in │
│ quotes, e.g. "{'key': │
│ value}". │
╰──────────────────────────────────────────────────────────────────────────────╯
╭─ Calculation ────────────────────────────────────────────────────────────────╮
│ * --struct <path> Path of structure │
│ to simulate. │
│ [required] │
│ --shear-magnitude <float> Magnitude of shear │
│ strain for deformed │
│ structures. │
│ [default: 0.06] │
│ --normal-magnitude <float> Magnitude of normal │
│ strain for deformed │
│ structures. │
│ [default: 0.01] │
│ --n-strains <int> Number of normal │
│ and shear strains │
│ to use for deformed │
│ structures. │
│ [default: 4] │
│ --write-voigt --no-write-voigt Write the │
│ ElasticityTensor in │
│ Voigt notation. │
│ [default: │
│ write-voigt] │
│ --minimize --no-minimize Whether to minimize │
│ initial structure │
│ before │
│ calculations. │
│ [default: minimize] │
│ --minimize-all --no-minimize-all Whether to minimize │
│ all generated │
│ structures for │
│ calculations. │
│ [default: │
│ no-minimize-all] │
│ --fmax <float> Maximum force for │
│ optimization │
│ convergence. │
│ [default: 0.1] │
│ --minimize-kwargs DICT Keyword arguments │
│ to pass to geometry │
│ optimizer, │
│ including │
│ "opt_kwargs", │
│ "filter_kwargs", │
│ and "traj_kwargs". │
│ Must be passed as a │
│ dictionary wrapped │
│ in quotes, e.g. │
│ "{'key': value}". │
│ --write-structures --no-write-structu… Whether to write │
│ out all genereated │
│ structures. │
│ [default: │
│ no-write-structure… │
╰──────────────────────────────────────────────────────────────────────────────╯
╭─ Structure I/O ──────────────────────────────────────────────────────────────╮
│ --file-prefix <path> Prefix for output files, including │
│ directories. Default directory is │
│ ./janus_results, and default filename prefix │
│ is inferred from the input stucture filename. │
│ --read-kwargs DICT Keyword arguments to pass to ase.io.read. Must │
│ be passed as a dictionary wrapped in quotes, │
│ e.g. "{'key': value}". By default, │
│ read_kwargs['index'] = -1, so only the last │
│ structure is read. │
│ --write-kwargs DICT Keyword arguments to pass to ase.io.write when │
│ saving any structures. Must be passed as a │
│ dictionary wrapped in quotes, e.g. "{'key': │
│ value}". │
╰──────────────────────────────────────────────────────────────────────────────╯
╭─ Logging/summary ────────────────────────────────────────────────────────────╮
│ --log <path> Path to save logs to. Default is │
│ inferred from `file_prefix` │
│ --tracker --no-tracker Whether to save carbon emissions of │
│ calculation │
│ [default: tracker] │
│ --summary <path> Path to save summary of inputs, │
│ start/end time, and carbon emissions. │
│ Default is inferred from `file_prefix`. │
╰──────────────────────────────────────────────────────────────────────────────╯
To calculate the elasticity tensor we can use the following CLI command for our Aluminium sample.
[8]:
! janus elasticity --arch mace_mp --struct ../data/Aluminium.xyz --no-tracker
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/e3nn/o3/_wigner.py:10: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
_Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))
cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.
Downloading MACE model from 'https://github.com/ACEsuit/mace-mp/releases/download/mace_mp_0/2023-12-10-mace-128-L0_energy_epoch-249.model'
Downloading: 100.0% (31.1 MB / 31.1 MB)
Cached MACE model to /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:274: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
torch.load(f=model_path, map_location=device)
Step Time Energy fmax
LBFGS: 0 12:22:59 -100.158373 0.019203
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.
[9]:
! cat janus_results/Aluminium-elastic_tensor.dat
# 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]
80.34637735135469 80.3463773513547 80.3463773513547 22.54802608157281 23.34462105588001 22.94632356872641 62.855302809644954 0.1766440599778729 0.36961597838074944 104.69105725175157 68.17403740115658 68.17403740115653 0.0 0.0 0.0 68.17403740115637 104.69105725175139 68.17403740115644 0.0 0.0 0.0 68.17403740115614 68.17403740115614 104.69105725175105 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601745 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601745 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601745
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.
For elasticity the main options for user control are the magnitudes of the applied shear and normal strains as well as their number.
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.
For example we can increase the shear and normal magnitudes, and apply 16 strains of each type as follows:
[10]:
! janus elasticity --arch mace_mp --struct ../data/Aluminium.xyz --no-tracker --n-strains 16 --shear-magnitude 0.3 --normal-magnitude 0.2
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/e3nn/o3/_wigner.py:10: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
_Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))
cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:274: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
torch.load(f=model_path, map_location=device)
Step Time Energy fmax
LBFGS: 0 12:23:09 -100.158373 0.019203
With the results now being different.
[11]:
! cat janus_results/Aluminium-elastic_tensor.dat
# 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]
96.0598935747555 96.05989357475555 96.05989357475553 28.115301864608714 94.23212116428371 61.17371151444621 151.38556869290915 11.758155686548445 0.2373417023843598 112.97498688614613 87.6023469190603 87.6023469190603 0.0 0.0 0.0 87.60234691906024 112.9749868861461 87.60234691906024 0.0 0.0 0.0 87.60234691906027 87.60234691906027 112.97498688614607 0.0 0.0 0.0 0.0 0.0 0.0 148.5959886181109 0.0 0.0 0.0 0.0 0.0 0.0 148.5959886181109 0.0 0.0 0.0 0.0 0.0 0.0 148.59598861811094
In the same way we can simply calculate the data for our Diamond and Carbon nano-tube samples,
[12]:
! janus elasticity --arch mace_mp --struct ../data/Diamond.xyz --no-tracker --write-structures
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/e3nn/o3/_wigner.py:10: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
_Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))
cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:274: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
torch.load(f=model_path, map_location=device)
Step Time Energy fmax
LBFGS: 0 12:23:31 -581.199344 0.165171
LBFGS: 1 12:23:32 -581.200507 0.163391
LBFGS: 2 12:23:33 -581.253268 0.012436
[13]:
! cat janus_results/Diamond-elastic_tensor.dat
# 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]
450.1943212514214 450.1943212514214 450.1943212514214 274.3844055824093 326.05315787990963 300.2187817311594 736.8608128396435 0.9415395198540546 0.22720638693999043 689.5017274807417 330.5406181366512 330.54061813665203 0.0 0.0 0.0 330.54061813670734 689.5017274810456 330.54061813662037 0.0 0.0 0.0 330.54061813670756 330.54061813662094 689.5017274810461 0.0 0.0 0.0 -1.0136478390474652e-10 1.1472701608614945e-10 -1.3813023198768173e-10 423.7682266850985 0.0 0.0 0.0 0.0 0.0 0.0 423.768226685073 0.0 0.0 0.0 0.0 0.0 0.0 423.76822668507316
Noticed that we supplied also the option --write-structures. This generates two new files in janus_results: Diamond-elasticity-generated.extxyz and Diamond-elasticity-opt.extxyz.
The latter is simply our initial diamond structure, but geometry optimized. The former is the set of strained structures used to calculate the elasticity. We can read the file and observe the impact of the strains like so
[14]:
v=WeasWidget()
diamond_strained = read("janus_results/Diamond-elasticity-generated.extxyz", index=":")
v.from_ase(diamond_strained)
v
[14]:
[15]:
! janus elasticity --arch mace_mp --struct ../data/Carbon-nanotube.xyz --no-tracker
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/e3nn/o3/_wigner.py:10: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
_Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))
cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:274: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
torch.load(f=model_path, map_location=device)
Step Time Energy fmax
LBFGS: 0 12:24:04 -849.799763 6.268367
LBFGS: 1 12:24:06 -856.222492 2.349892
LBFGS: 2 12:24:07 -858.847051 1.109023
LBFGS: 3 12:24:08 -859.823886 1.067070
LBFGS: 4 12:24:08 -860.237243 1.028980
LBFGS: 5 12:24:09 -860.538898 1.005038
LBFGS: 6 12:24:10 -860.681007 0.988611
LBFGS: 7 12:24:11 -860.934222 0.974843
LBFGS: 8 12:24:12 -861.275268 0.962106
LBFGS: 9 12:24:13 -861.676469 0.952216
LBFGS: 10 12:24:14 -862.066407 0.940885
LBFGS: 11 12:24:15 -862.447240 0.925222
LBFGS: 12 12:24:16 -862.814264 0.903153
LBFGS: 13 12:24:17 -863.170981 0.878710
LBFGS: 14 12:24:18 -863.530212 0.913888
LBFGS: 15 12:24:18 -863.905406 0.925973
LBFGS: 16 12:24:19 -864.294730 0.864437
LBFGS: 17 12:24:20 -864.669249 0.695976
LBFGS: 18 12:24:21 -864.983970 0.648393
LBFGS: 19 12:24:22 -865.207944 0.605894
LBFGS: 20 12:24:23 -865.354798 0.572142
LBFGS: 21 12:24:24 -865.489388 0.554320
LBFGS: 22 12:24:24 -865.728634 0.538098
LBFGS: 23 12:24:25 -865.993684 0.645311
LBFGS: 24 12:24:26 -866.282437 0.749929
LBFGS: 25 12:24:27 -866.577762 0.756553
LBFGS: 26 12:24:28 -866.863855 0.709103
LBFGS: 27 12:24:29 -867.129347 0.667180
LBFGS: 28 12:24:30 -867.368495 0.711989
LBFGS: 29 12:24:30 -867.579368 0.760721
LBFGS: 30 12:24:31 -867.763200 0.812266
LBFGS: 31 12:24:32 -867.922977 0.865492
LBFGS: 32 12:24:33 -868.062776 0.919399
LBFGS: 33 12:24:34 -868.187302 0.972989
LBFGS: 34 12:24:34 -868.301853 1.025241
LBFGS: 35 12:24:35 -868.412367 1.074994
LBFGS: 36 12:24:36 -868.525549 1.120851
LBFGS: 37 12:24:37 -868.649003 1.160999
LBFGS: 38 12:24:38 -868.791688 1.192842
LBFGS: 39 12:24:38 -868.965389 1.212000
LBFGS: 40 12:24:39 -869.189376 1.208879
LBFGS: 41 12:24:40 -869.464759 1.159925
LBFGS: 42 12:24:41 -869.655324 1.086800
LBFGS: 43 12:24:42 -869.820087 0.999967
LBFGS: 44 12:24:43 -869.966942 0.904172
LBFGS: 45 12:24:43 -870.094124 0.801747
LBFGS: 46 12:24:44 -870.202430 0.694082
LBFGS: 47 12:24:45 -870.302706 0.582461
LBFGS: 48 12:24:46 -870.416439 0.468606
LBFGS: 49 12:24:47 -870.562834 0.414222
LBFGS: 50 12:24:48 -870.748578 0.403267
LBFGS: 51 12:24:48 -870.969965 0.338580
LBFGS: 52 12:24:49 -871.133062 0.273087
LBFGS: 53 12:24:50 -871.246834 0.190993
LBFGS: 54 12:24:51 -871.310110 0.102759
LBFGS: 55 12:24:52 -871.321534 0.082324
[16]:
! cat janus_results/Carbon-nanotube-elastic_tensor.dat
# 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]
193.9503009215358 218.54266486979637 206.24648289566608 77.9064140507454 139.5610079179821 108.73371098436374 277.4446128784308 4.083762446676342 0.27579850980312903 235.08150639637853 178.51925017896437 113.66022222756877 1.492422534759681e-10 -1.3045736181331495e-09 18.447300674106685 178.51925017897062 235.08150639627286 113.66022222748538 -3.7253708983672235e-10 -2.764491875272417e-09 -18.44730067325802 116.15034813601955 116.15034813592992 680.0613299505777 -1.6191960643265413e-09 4.624288414024965e-10 0.0 0.0 -3.1239908199942987e-10 -1.1933931504160093e-09 135.05386607246464 -1.185311510799643 1.1726378343734233e-09 3.253917258159012e-09 -8.701226262964039e-10 -1.6737687911431308e-09 1.1853115107937473 135.0538660724815 1.1539382576114167e-09 20.608899695476488 -20.608899695023567 -1.0362459038068705e-10 -2.763023995162113e-10 7.394001549643116e-10 179.5690914085609
We can examine the results to asses the comparative stiffness of the materials, do the results line up with experiments?
[17]:
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, layout='constrained')
for material, marker in zip(('Aluminium', 'Diamond', 'Carbon-nanotube'), "o^s"):
tensor = np.loadtxt(f"janus_results/{material}-elastic_tensor.dat")
bulk, shear, youngs = tensor[[0, 3, 6]]
theta = [i*2.0*np.pi/3 for i in range(3)]
ax.scatter(theta, [bulk, shear,youngs], marker=marker, label=material)
ax.set_xticks(theta, ("Bulk", "Shear", "Youngs"))
fig.legend(loc='outside upper right')
fig.suptitle("Elastic modulii for our samples")
[17]:
Text(0.5, 0.98, 'Elastic modulii for our samples')
