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.
This will replace pre-installed versions of numpy and torch in Colab with versions that are known to be compatible with janus-core.
It may be possible to skip the steps that uninstall and reinstall torch, which will save a considerable amount of time.
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]:
# import locale
# locale.getpreferredencoding = lambda: "UTF-8"
# ! pip uninstall numpy -y # Uninstall pre-installed numpy
# ! pip uninstall torch torchaudio torchvision transformers -y # Uninstall pre-installed torch
# ! uv pip install torch==2.5.1 # Install pinned version of torch
# ! uv pip install janus-core[mace,visualise] data-tutorials --system # Install janus-core with MACE and WeasWidget, and data-tutorials
# get_ipython().kernel.do_shutdown(restart=True) # Restart kernel to update libraries. This may warn that your session has crashed.
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.8.7
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 TEXT Path to configuration file. │
│ --help Show this message and exit. │
╰──────────────────────────────────────────────────────────────────────────────╯
╭─ MLIP calculator ────────────────────────────────────────────────────────────╮
│ * --arch [mace|mace_mp|mace_off|m MLIP architecture to use │
│ 3gnet|chgnet|alignn|seve for calculations. │
│ nnet|nequip|dpa3|orb|mat [required] │
│ tersim|grace|esen|equifo │
│ rmer|pet_mad|uma|mace_om │
│ ol] │
│ --device [cpu|cuda|mps|xpu] Device to run calculations │
│ on. │
│ [default: cpu] │
│ --model TEXT 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 INTEGER 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'
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:197: 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)
Using head Default out of ['Default']
Step Time Energy fmax
LBFGS: 0 13:21:04 -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:197: 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)
Using head Default out of ['Default']
Step Time Energy fmax
LBFGS: 0 13:21:13 -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:197: 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)
Using head Default out of ['Default']
Step Time Energy fmax
LBFGS: 0 13:21:32 -581.199344 0.165171
LBFGS: 1 13:21:33 -581.200507 0.163391
LBFGS: 2 13:21: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.1943212512929 450.1943212512928 450.1943212512929 274.38440558240325 326.05315787984875 300.218781731126 736.8608128395381 0.9415395198530749 0.2272063869399515 689.5017274808512 330.5406181365162 330.54061813651265 0.0 0.0 0.0 330.54061813651333 689.5017274808511 330.5406181365137 0.0 0.0 0.0 330.5406181365131 330.5406181365129 689.5017274808517 0.0 0.0 0.0 1.0783876082080196e-10 -1.1332034214307304e-10 1.375591364101266e-10 423.76822668496874 0.0 0.0 0.0 0.0 0.0 0.0 423.7682266849697 0.0 0.0 0.0 0.0 0.0 0.0 423.76822668496845
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:197: 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)
Using head Default out of ['Default']
Step Time Energy fmax
LBFGS: 0 13:21:54 -849.799763 6.268367
LBFGS: 1 13:21:54 -856.222492 2.349892
LBFGS: 2 13:21:55 -858.847051 1.109023
LBFGS: 3 13:21:56 -859.823886 1.067070
LBFGS: 4 13:21:56 -860.237243 1.028980
LBFGS: 5 13:21:57 -860.538898 1.005038
LBFGS: 6 13:21:57 -860.681007 0.988611
LBFGS: 7 13:21:58 -860.934222 0.974843
LBFGS: 8 13:21:59 -861.275268 0.962106
LBFGS: 9 13:21:59 -861.676469 0.952216
LBFGS: 10 13:22:00 -862.066407 0.940885
LBFGS: 11 13:22:00 -862.447240 0.925222
LBFGS: 12 13:22:01 -862.814264 0.903153
LBFGS: 13 13:22:02 -863.170981 0.878710
LBFGS: 14 13:22:02 -863.530212 0.913888
LBFGS: 15 13:22:03 -863.905406 0.925973
LBFGS: 16 13:22:03 -864.294730 0.864437
LBFGS: 17 13:22:04 -864.669249 0.695976
LBFGS: 18 13:22:05 -864.983970 0.648393
LBFGS: 19 13:22:05 -865.207944 0.605894
LBFGS: 20 13:22:06 -865.354798 0.572142
LBFGS: 21 13:22:06 -865.489388 0.554320
LBFGS: 22 13:22:07 -865.728634 0.538098
LBFGS: 23 13:22:07 -865.993684 0.645311
LBFGS: 24 13:22:08 -866.282437 0.749929
LBFGS: 25 13:22:09 -866.577762 0.756553
LBFGS: 26 13:22:09 -866.863855 0.709103
LBFGS: 27 13:22:10 -867.129347 0.667180
LBFGS: 28 13:22:10 -867.368495 0.711989
LBFGS: 29 13:22:11 -867.579368 0.760721
LBFGS: 30 13:22:11 -867.763200 0.812266
LBFGS: 31 13:22:12 -867.922977 0.865492
LBFGS: 32 13:22:13 -868.062776 0.919399
LBFGS: 33 13:22:13 -868.187302 0.972989
LBFGS: 34 13:22:14 -868.301853 1.025241
LBFGS: 35 13:22:14 -868.412367 1.074994
LBFGS: 36 13:22:15 -868.525549 1.120851
LBFGS: 37 13:22:15 -868.649003 1.160999
LBFGS: 38 13:22:16 -868.791688 1.192842
LBFGS: 39 13:22:16 -868.965389 1.212000
LBFGS: 40 13:22:17 -869.189376 1.208879
LBFGS: 41 13:22:17 -869.464759 1.159925
LBFGS: 42 13:22:18 -869.655324 1.086800
LBFGS: 43 13:22:18 -869.820087 0.999967
LBFGS: 44 13:22:19 -869.966942 0.904172
LBFGS: 45 13:22:20 -870.094124 0.801747
LBFGS: 46 13:22:20 -870.202430 0.694082
LBFGS: 47 13:22:21 -870.302706 0.582461
LBFGS: 48 13:22:21 -870.416439 0.468606
LBFGS: 49 13:22:22 -870.562834 0.414222
LBFGS: 50 13:22:22 -870.748578 0.403267
LBFGS: 51 13:22:23 -870.969965 0.338580
LBFGS: 52 13:22:23 -871.133062 0.273087
LBFGS: 53 13:22:24 -871.246834 0.190993
LBFGS: 54 13:22:24 -871.310110 0.102759
LBFGS: 55 13:22:25 -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.95030092167832 218.54266486996087 206.2464828958196 77.90641405101773 139.56100791800858 108.73371098451315 277.44461287878585 4.0837624466467535 0.2757985098030089 235.08150639466425 178.5192501790822 113.6602222286648 0.0 -1.323096452128578e-09 18.447300697318852 178.51925017916716 235.08150639821588 113.66022222669963 -3.7230073771373626e-10 -2.6791053464482114e-09 -18.447300650233707 116.15034813714018 116.15034813508179 680.0613299509321 -1.5752902402338083e-09 4.763653570662971e-10 -1.9996008173953012e-09 -2.0380856396216894e-10 0.0 -1.195651786614127e-09 135.0538660731713 -1.185311510077846 1.186370459485761e-09 3.1845303157638846e-09 -8.095139378505415e-10 -1.6490127340783737e-09 1.1853115115172188 135.05386607176266 1.2362876360267125e-09 20.608899705385046 -20.60889968096377 -4.224279299948755e-09 -4.912058798658525e-10 8.248188509099604e-10 179.56909140865375
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')
