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 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.4
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 <function> MLIP architecture to use for │
│ calculations. │
│ [required] │
│ --device <function> 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:226: 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 19:58:52 -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.34637735135448 80.34637735135449 80.34637735135448 22.548026081572864 23.344621055880047 22.946323568726456 62.85530280964503 0.17664405997786936 0.3696159783807489 104.69105725175152 68.17403740115631 68.1740374011563 0.0 0.0 0.0 68.17403740115608 104.69105725175119 68.17403740115606 0.0 0.0 0.0 68.17403740115597 68.17403740115597 104.69105725175099 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601734 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601756 0.0 0.0 0.0 0.0 0.0 0.0 26.735361809601756
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:226: 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 19:59:02 -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.05989357475553 96.0598935747555 96.05989357475552 28.11530186460873 94.23212116428365 61.17371151444619 151.38556869290912 11.758155686548424 0.23734170238435986 112.9749868861461 87.60234691906024 87.60234691906025 0.0 0.0 0.0 87.60234691906021 112.97498688614604 87.60234691906024 0.0 0.0 0.0 87.60234691906021 87.60234691906018 112.97498688614604 0.0 0.0 0.0 0.0 0.0 0.0 148.59598861811085 0.0 0.0 0.0 0.0 0.0 0.0 148.5959886181108 0.0 0.0 0.0 0.0 0.0 0.0 148.59598861811085
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:226: 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 19:59:22 -581.199344 0.165171
LBFGS: 1 19:59:23 -581.200507 0.163391
LBFGS: 2 19:59:24 -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.19432125130555 450.1943212513056 450.1943212513056 274.38440558240325 326.05315787985336 300.2187817311283 736.8608128395465 0.9415395198531602 0.22720638693995612 689.5017274808616 330.5406181365285 330.540618136529 0.0 0.0 0.0 330.5406181365284 689.5017274808603 330.5406181365268 0.0 0.0 0.0 330.5406181365283 330.54061813652737 689.5017274808603 0.0 0.0 0.0 0.0 0.0 0.0 423.76822668497834 0.0 0.0 0.0 0.0 0.0 0.0 423.7682266849781 0.0 0.0 0.0 0.0 0.0 0.0 423.7682266849775
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:226: 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 19:59:45 -849.799763 6.268367
LBFGS: 1 19:59:46 -856.222492 2.349892
LBFGS: 2 19:59:47 -858.847051 1.109023
LBFGS: 3 19:59:48 -859.823886 1.067070
LBFGS: 4 19:59:48 -860.237243 1.028980
LBFGS: 5 19:59:49 -860.538898 1.005038
LBFGS: 6 19:59:49 -860.681007 0.988611
LBFGS: 7 19:59:50 -860.934222 0.974843
LBFGS: 8 19:59:51 -861.275268 0.962106
LBFGS: 9 19:59:51 -861.676469 0.952216
LBFGS: 10 19:59:52 -862.066407 0.940885
LBFGS: 11 19:59:52 -862.447240 0.925222
LBFGS: 12 19:59:53 -862.814264 0.903153
LBFGS: 13 19:59:54 -863.170981 0.878710
LBFGS: 14 19:59:54 -863.530212 0.913888
LBFGS: 15 19:59:55 -863.905406 0.925973
LBFGS: 16 19:59:55 -864.294730 0.864437
LBFGS: 17 19:59:56 -864.669249 0.695976
LBFGS: 18 19:59:56 -864.983970 0.648393
LBFGS: 19 19:59:57 -865.207944 0.605894
LBFGS: 20 19:59:58 -865.354798 0.572142
LBFGS: 21 19:59:58 -865.489388 0.554320
LBFGS: 22 19:59:59 -865.728634 0.538098
LBFGS: 23 19:59:59 -865.993684 0.645311
LBFGS: 24 20:00:00 -866.282437 0.749929
LBFGS: 25 20:00:00 -866.577762 0.756553
LBFGS: 26 20:00:01 -866.863855 0.709103
LBFGS: 27 20:00:02 -867.129347 0.667180
LBFGS: 28 20:00:02 -867.368495 0.711989
LBFGS: 29 20:00:03 -867.579368 0.760721
LBFGS: 30 20:00:03 -867.763200 0.812266
LBFGS: 31 20:00:04 -867.922977 0.865492
LBFGS: 32 20:00:04 -868.062776 0.919399
LBFGS: 33 20:00:05 -868.187302 0.972989
LBFGS: 34 20:00:05 -868.301853 1.025241
LBFGS: 35 20:00:06 -868.412367 1.074994
LBFGS: 36 20:00:07 -868.525549 1.120851
LBFGS: 37 20:00:07 -868.649003 1.160999
LBFGS: 38 20:00:08 -868.791688 1.192842
LBFGS: 39 20:00:08 -868.965389 1.212000
LBFGS: 40 20:00:09 -869.189376 1.208879
LBFGS: 41 20:00:09 -869.464759 1.159925
LBFGS: 42 20:00:10 -869.655324 1.086800
LBFGS: 43 20:00:10 -869.820087 0.999967
LBFGS: 44 20:00:11 -869.966942 0.904172
LBFGS: 45 20:00:11 -870.094124 0.801747
LBFGS: 46 20:00:12 -870.202430 0.694082
LBFGS: 47 20:00:13 -870.302706 0.582461
LBFGS: 48 20:00:13 -870.416439 0.468606
LBFGS: 49 20:00:14 -870.562834 0.414222
LBFGS: 50 20:00:14 -870.748578 0.403267
LBFGS: 51 20:00:15 -870.969965 0.338580
LBFGS: 52 20:00:15 -871.133062 0.273087
LBFGS: 53 20:00:16 -871.246834 0.190993
LBFGS: 54 20:00:16 -871.310110 0.102759
LBFGS: 55 20:00:17 -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.95030092156225 218.5426648698607 206.24648289571147 77.90641405075749 139.56100791795427 108.73371098435588 277.4446128784228 4.0837624466733455 0.2757985098031848 235.08150639601578 178.51925017896647 113.66022222779782 1.1485431457546693e-10 -1.308654257039186e-09 18.44730067703937 178.51925017898978 235.08150639659797 113.66022222754263 -3.793191767913997e-10 -2.691046251396711e-09 -18.44730067033899 116.15034813619116 116.15034813591083 680.0613299507339 -1.6193395659099362e-09 4.589392098610714e-10 -2.6270744589619876e-10 -1.0185161947905582e-10 -2.6864473475838133e-10 -1.2020836335285596e-09 135.05386607255127 -1.1853115106897139 1.156772524213653e-09 3.177430386497328e-09 -8.295537752891772e-10 -1.6627447740518e-09 1.1853115109030181 135.05386607230332 1.17904686087031e-09 20.608899696695698 -20.608899693218863 -5.83365940527819e-10 -2.575670967107553e-10 7.642173668023703e-10 179.5690914085698
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')
