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 data and modules¶
[3]:
from weas_widget import WeasWidget
from ase.build import bulk, nanotube
from ase.lattice.cubic import Diamond
from ase.io import read
import numpy as np
import matplotlib.pyplot as plt
from janus_core.calculations.elasticity import Elasticity
Use data_tutorials to get the data required for this tutorial:
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. Then compare their relative stiffnesses.
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)
v=WeasWidget()
v.from_ase(al)
v
[4]:
Next we can build Diamond,
[5]:
diamond = Diamond('C')*(2,2,2)
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]
v=WeasWidget()
v.from_ase(nt)
v
[6]:
Elasticity calculations¶
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:
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.
[7]:
elasticity_aluminium = Elasticity(
struct=al.copy(),
arch="mace_mp",
device="cpu",
model="small",
calc_kwargs={"default_dtype": "float64"}
).run()
/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.
Using head Default out of ['Default']
/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)
Step Time Energy fmax
LBFGS: 0 13:38:39 -100.158373 0.019203
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.
[8]:
! cat janus_results/Al27-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.34637735135443 80.34637735135442 80.34637735135442 22.54802608157287 23.344621055880058 22.946323568726463 62.85530280964506 0.17664405997786936 0.3696159783807488 104.6910572517515 68.17403740115623 68.17403740115626 0.0 0.0 0.0 68.17403740115591 104.69105725175115 68.17403740115606 0.0 0.0 0.0 68.17403740115594 68.17403740115589 104.69105725175092 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.735361809601763 0.0 0.0 0.0 0.0 0.0 0.0 26.73536180960177
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:
[9]:
calc_elasticity_aluminium = Elasticity(
struct=al.copy(),
arch="mace_mp",
device="cpu",
model="small",
calc_kwargs={"default_dtype": "float64"},
shear_magnitude=0.3,
normal_magnitude=0.2,
n_strains=16
)
elasticity_aluminium = calc_elasticity_aluminium.run()
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.
Using head Default out of ['Default']
/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)
Step Time Energy fmax
LBFGS: 0 13:38:43 -100.158373 0.019203
Using Weas we can get a feel for how the applied deformation are impacting our structure by playing through them like any other trajectory.
In this case, since we saved the calculation object we can directly access the deformed_structures
[10]:
v=WeasWidget()
v.from_ase(calc_elasticity_aluminium.deformed_structures)
v
[10]:
Also note that the deformed structure contain the applied strain within their info attributes
[11]:
s = calc_elasticity_aluminium.deformed_structures[-1].info["strain"]
print(f"(Voigt) Strain recorded in structure: {calc_elasticity_aluminium.strains[-1].voigt}")
print(f"Original strain {calc_elasticity_aluminium.strains[-1]}")
(Voigt) Strain recorded in structure: [0. 0. 0. 0.6 0. 0. ]
Original strain [[0. 0. 0. ]
[0. 0. 0.3]
[0. 0.3 0. ]]
[12]:
elasticity_diamond = Elasticity(
struct=diamond.copy(),
arch="mace_mp",
device="cpu",
model="small",
write_structures=True,
calc_kwargs={"default_dtype": "float64"}
).run()
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.
Using head Default out of ['Default']
/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)
Step Time Energy fmax
LBFGS: 0 13:38:56 -581.199344 0.165171
LBFGS: 1 13:38:56 -581.200507 0.163391
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/scipy/_lib/_util.py:1181: RuntimeWarning: logm result may be inaccurate, approximate err = 7.169467197482553e-13
return f(*arrays, *other_args, **kwargs)
LBFGS: 2 13:38:57 -581.253268 0.012436
Notice for diamond we used the write_structures=True key word argument, in this case our results directory will include the file C64-elasticity-generated.extxyz and C64-elasticity-opt.extxyz.
The latter is simply our initial diamond structure, but geometry optimized. The former is again the set of strained structures used to calculate the elasticity. We can read the file and observe the impact of the strains like so
[13]:
v=WeasWidget()
diamond_strained = read("janus_results/C64-elasticity-generated.extxyz", index=":")
v.from_ase(diamond_strained)
v
[13]:
[14]:
elasticity_nt = Elasticity(
struct=nt.copy(),
arch="mace_mp",
device="cpu",
model="small",
calc_kwargs={"default_dtype": "float64"}
).run()
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.
Using head Default out of ['Default']
/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)
Step Time Energy fmax
LBFGS: 0 13:39:12 -849.799763 6.268367
LBFGS: 1 13:39:13 -856.222492 2.349892
LBFGS: 2 13:39:13 -858.847051 1.109023
LBFGS: 3 13:39:14 -859.823886 1.067070
LBFGS: 4 13:39:14 -860.237243 1.028980
LBFGS: 5 13:39:15 -860.538898 1.005038
LBFGS: 6 13:39:16 -860.681007 0.988611
LBFGS: 7 13:39:16 -860.934222 0.974843
LBFGS: 8 13:39:17 -861.275268 0.962106
LBFGS: 9 13:39:17 -861.676469 0.952216
LBFGS: 10 13:39:18 -862.066407 0.940885
LBFGS: 11 13:39:19 -862.447240 0.925222
LBFGS: 12 13:39:19 -862.814264 0.903153
LBFGS: 13 13:39:20 -863.170981 0.878710
LBFGS: 14 13:39:20 -863.530212 0.913888
LBFGS: 15 13:39:21 -863.905406 0.925973
LBFGS: 16 13:39:22 -864.294730 0.864437
LBFGS: 17 13:39:22 -864.669249 0.695976
LBFGS: 18 13:39:23 -864.983970 0.648393
LBFGS: 19 13:39:23 -865.207944 0.605894
LBFGS: 20 13:39:24 -865.354798 0.572142
LBFGS: 21 13:39:25 -865.489388 0.554320
LBFGS: 22 13:39:25 -865.728634 0.538098
LBFGS: 23 13:39:26 -865.993684 0.645311
LBFGS: 24 13:39:26 -866.282437 0.749929
LBFGS: 25 13:39:27 -866.577762 0.756553
LBFGS: 26 13:39:27 -866.863856 0.709103
LBFGS: 27 13:39:28 -867.129347 0.667180
LBFGS: 28 13:39:28 -867.368495 0.711989
LBFGS: 29 13:39:29 -867.579368 0.760721
LBFGS: 30 13:39:30 -867.763200 0.812266
LBFGS: 31 13:39:30 -867.922977 0.865492
LBFGS: 32 13:39:31 -868.062776 0.919399
LBFGS: 33 13:39:31 -868.187302 0.972989
LBFGS: 34 13:39:32 -868.301853 1.025241
LBFGS: 35 13:39:32 -868.412367 1.074994
LBFGS: 36 13:39:33 -868.525549 1.120851
LBFGS: 37 13:39:33 -868.649003 1.160999
LBFGS: 38 13:39:34 -868.791688 1.192842
LBFGS: 39 13:39:34 -868.965389 1.212000
LBFGS: 40 13:39:35 -869.189376 1.208879
LBFGS: 41 13:39:36 -869.464759 1.159925
LBFGS: 42 13:39:36 -869.655324 1.086800
LBFGS: 43 13:39:37 -869.820087 0.999967
LBFGS: 44 13:39:37 -869.966942 0.904172
LBFGS: 45 13:39:38 -870.094124 0.801747
LBFGS: 46 13:39:38 -870.202430 0.694082
LBFGS: 47 13:39:39 -870.302706 0.582461
LBFGS: 48 13:39:39 -870.416439 0.468606
LBFGS: 49 13:39:40 -870.562834 0.414222
LBFGS: 50 13:39:40 -870.748578 0.403267
LBFGS: 51 13:39:41 -870.969965 0.338580
LBFGS: 52 13:39:41 -871.133062 0.273087
LBFGS: 53 13:39:42 -871.246834 0.190993
LBFGS: 54 13:39:43 -871.310110 0.102759
LBFGS: 55 13:39:43 -871.321534 0.082324
[15]:
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, layout='constrained')
for tensor, material, marker in zip((elasticity_aluminium, elasticity_diamond, elasticity_nt), ("Aluminium", "Diamond", "Carbon-nanotube"), "o^s"):
ecs = dict()
ecs["Bulk modulus"] = tensor.k_voigt
ecs["Shear modulus"] = tensor.g_voigt
# Note https://github.com/materialsproject/pymatgen/issues/4435
ecs["Young's modulus"] = tensor.y_mod/1e9
theta = [i*2.0*np.pi/(len(ecs)) for i in range(len(ecs))]
ax.scatter(theta, ecs.values(), marker=marker, label=material)
ax.set_xticks(theta, ecs.keys())
fig.legend(loc='outside upper right')
fig.suptitle("Elastic modulii for our samples")
[15]:
Text(0.5, 0.98, 'Elastic modulii for our samples')
