Elastic Constants from Stress-Strain Simulations¶

This notebook provides the theoretical background and a practical example for calculating the elastic constants of a glass system using the Stress-Strain (Finite Difference) method.

Elastic Constants via Finite Differences¶

1. Theoretical Background¶

In the framework of linear elasticity, the relationship between the stress tensor ($\sigma_{ij}$) and the strain tensor ($\epsilon_{ij}$) is governed by Hooke's Law:

$$\sigma_{ij} = C_{ijkl}\epsilon_{kl}$$

Using Voigt notation, the $3\times3\times3$ stiffness tensor $C_{ijkl}$ is reduced to a $6\times6$ matrix ($C_{ij}$), where indices 1-3 represent normal components and 4-6 represent shear components.

The Stress-Strain Method¶

To determine a specific $C_{ij}$, we apply a small deformation (strain) to the simulation cell and measure the resulting change in the internal stress tensor. We use a central difference scheme to improve accuracy and cancel out first-order errors:

$$C_{ij} = \frac{\sigma_i(+\epsilon) - \sigma_i(-\epsilon)}{2\epsilon}$$

Isotropic Averaging¶

For glass systems, which are macroscopically isotropic, the 21 independent components of $C_{ij}$ reduce to just two independent constants (typically the Bulk modulus $B$ and Shear modulus $G$). We use the Voigt-Reuss-Hill (VRH) average to estimate these:

• Voigt Bound ($G_V$): Assumes uniform strain.
• Reuss Bound ($G_R$): Assumes uniform stress.
• Hill Average ($G$): The arithmetic mean of Voigt and Reuss bounds, providing the best estimate for polycrystals and glasses.

---

2. Practical Implementation¶

Below is an example of how to run the elastic_simulation workflow.

Prerequisites¶

structure of glass pre_quenched.

3. Workflow Notes¶

• Convergence: Elastic constants in MD are highly sensitive to statistical noise. Ensure your production_steps are long enough to get a converged average of the pressure tensor.

• Strain Magnitude: A strain of $10^{-3}$ (0.1%) is usually small enough to stay in the linear elastic regime but large enough to overcome thermal noise.

• Symmetry: Glasses are suposed to be isotropic and any significant deviation between $C_{11}$, $C_{22}$ and $C_{33}$ may indicate the system size is too small or not well-equilibrated. The code checks for cubic symmetry and return a warning if the respective elastic constants are not symmetric.

Important: In the example below a previously generated silica glass using PMMCS potential is used as an example of how to use the workflow. For production simulations, system size, cooling rate, equilibration time, and strain magnitude should be tested to ensure the robustness of the results.

In [1]:

from ase.io import read
from executorlib import SingleNodeExecutor


from amorphouspy import (
    generate_potential,
    get_structure_dict,
    elastic_simulation,
)

from ase.io import read from executorlib import SingleNodeExecutor from amorphouspy import ( generate_potential, get_structure_dict, elastic_simulation, )

In [2]:

# Project setup,
# generating an initial random structure
# and setting up the potential information
with SingleNodeExecutor(plot_dependency_graph=True) as exe:
    atoms_dict_future = exe.submit(
        get_structure_dict,  # this line is needed to generate the atoms_dict so the potential can be generated
        composition={"SiO2": 100},  #
        target_atoms=30,
    )

    # Reading in the glass structure
    structure = read("data/SiO2_glass_300_atoms.xyz")

    # Generating the potential
    generated_potential_future = exe.submit(
        generate_potential,
        atoms_dict=atoms_dict_future,
        potential_type="pmmcs",
    )
    # Specification of the cpu parameters
    ncpus = 2
    server_kwargs = {"cores": ncpus}

    # Specification of the melt-quenching parameters
    elastic_sim = exe.submit(
        elastic_simulation,
        structure=structure,
        potential=generated_potential_future,
        temperature_sim=300.0,
        pressure=None,
        timestep=1.0,
        equilibration_steps=10_000,  # in production settings I would increase this to 1_000_000
        production_steps=1_000,  # in production settings I would increase this to at least 10_000
        n_print=1,
        strain=1e-3,  # strain magnitude should be tested for convergence this values is just a good starting point
        server_kwargs=server_kwargs,
        langevin=False,
        seed=12345,
    )

# Project setup, # generating an initial random structure # and setting up the potential information with SingleNodeExecutor(plot_dependency_graph=True) as exe: atoms_dict_future = exe.submit( get_structure_dict, # this line is needed to generate the atoms_dict so the potential can be generated composition={"SiO2": 100}, # target_atoms=30, ) # Reading in the glass structure structure = read("data/SiO2_glass_300_atoms.xyz") # Generating the potential generated_potential_future = exe.submit( generate_potential, atoms_dict=atoms_dict_future, potential_type="pmmcs", ) # Specification of the cpu parameters ncpus = 2 server_kwargs = {"cores": ncpus} # Specification of the melt-quenching parameters elastic_sim = exe.submit( elastic_simulation, structure=structure, potential=generated_potential_future, temperature_sim=300.0, pressure=None, timestep=1.0, equilibration_steps=10_000, # in production settings I would increase this to 1_000_000 production_steps=1_000, # in production settings I would increase this to at least 10_000 n_print=1, strain=1e-3, # strain magnitude should be tested for convergence this values is just a good starting point server_kwargs=server_kwargs, langevin=False, seed=12345, )

In [3]:

# Project setup,
# generating an initial random structure
# and setting up the potential information
with SingleNodeExecutor() as exe:
    atoms_dict_future = exe.submit(
        get_structure_dict,  # this line is needed to generate the atoms_dict so the potential can be generated
        composition={"SiO2": 100},  #
        target_atoms=30,
    )

    # Reading in the glass structure
    structure = read("data/SiO2_glass_300_atoms.xyz")

    # Generating the potential
    generated_potential_future = exe.submit(
        generate_potential,
        atoms_dict=atoms_dict_future,
        potential_type="pmmcs",
    )
    # Specification of the cpu parameters
    ncpus = 2
    server_kwargs = {"cores": ncpus}

    # Specification of the melt-quenching parameters
    elastic_sim_future = exe.submit(
        elastic_simulation,
        structure=structure,
        potential=generated_potential_future,
        temperature_sim=300.0,
        pressure=None,
        timestep=1.0,
        equilibration_steps=10_000,  # in production settings I would increase this to 1_000_000
        production_steps=1_000,  # in production settings I would increase this to at least 10_000
        n_print=1,
        strain=1e-3,  # strain magnitude should be tested for convergence this values is just a good starting point
        server_kwargs=server_kwargs,
        langevin=False,
        seed=12345,
    )
    results = elastic_sim_future.result()  # Running the workflow

# Project setup, # generating an initial random structure # and setting up the potential information with SingleNodeExecutor() as exe: atoms_dict_future = exe.submit( get_structure_dict, # this line is needed to generate the atoms_dict so the potential can be generated composition={"SiO2": 100}, # target_atoms=30, ) # Reading in the glass structure structure = read("data/SiO2_glass_300_atoms.xyz") # Generating the potential generated_potential_future = exe.submit( generate_potential, atoms_dict=atoms_dict_future, potential_type="pmmcs", ) # Specification of the cpu parameters ncpus = 2 server_kwargs = {"cores": ncpus} # Specification of the melt-quenching parameters elastic_sim_future = exe.submit( elastic_simulation, structure=structure, potential=generated_potential_future, temperature_sim=300.0, pressure=None, timestep=1.0, equilibration_steps=10_000, # in production settings I would increase this to 1_000_000 production_steps=1_000, # in production settings I would increase this to at least 10_000 n_print=1, strain=1e-3, # strain magnitude should be tested for convergence this values is just a good starting point server_kwargs=server_kwargs, langevin=False, seed=12345, ) results = elastic_sim_future.result() # Running the workflow

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning: Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity
  warnings.warn(
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/executorlib/standalone/interactive/backend.py:22: UserWarning: System may not be cubic: C11 != C22 != C33
  return args[0].__call__(*args[1:], **kwargs)
/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/executorlib/standalone/interactive/backend.py:22: UserWarning: System may not be cubic: C44 != C55 != C66
  return args[0].__call__(*args[1:], **kwargs)

In [4]:

moduli = results["moduli"]
E = moduli["E"]
G = moduli["G"]
B = moduli["B"]
nu = moduli["nu"]

print(f"Young's modulus: {E} GPa")
print(f"Shear modulus: {G} GPa")
print(f"Bulk modulus: {B} GPa")
print(f"Poisson's ratio: {nu}")

moduli = results["moduli"] E = moduli["E"] G = moduli["G"] B = moduli["B"] nu = moduli["nu"] print(f"Young's modulus: {E} GPa") print(f"Shear modulus: {G} GPa") print(f"Bulk modulus: {B} GPa") print(f"Poisson's ratio: {nu}")

Young's modulus: 80.06022877384835 GPa
Shear modulus: 33.214140713552425 GPa
Bulk modulus: 45.26445113073877 GPa
Poisson's ratio: 0.2052130064768052