Elastic constants¶

mdapy.lattice_dynamics.elastic_constants(system, calculator, *, strain=1e-3, relax=True, max_force=1e-4, max_steps=5000) gives the elastic constants of a periodic crystal from its stress under small strains.

from mdapy.lattice_dynamics import elastic_constants
from mdapy.simulation import minimize

relaxed = minimize(crystal, potential, cell="full", max_stress=1e-4).system
result = elastic_constants(relaxed, potential)
result.constants  # (6, 6) GPa
result.bulk_modulus(), result.shear_modulus("voigt")

Relax the box first

The constants are those of the box as it is given. Relax it to zero stress first, as above, or the constants belong to a crystal under stress.

Name	Type	Default	Unit	Meaning
`system`	`System`			the crystal, periodic in every direction, or a sheet, with a `mass` or an `element` column when its positions relax
`calculator`	`Calculator`			computes the virial, and energies and forces when the positions relax
`strain`	float64	1e-3		the magnitude \(\delta\) of the strains
`relax`	bool	`True`		whether the positions relax in each strained box (relaxed-ion constants) or stay put (clamped-ion constants)
`max_force`	float64	1e-4	eV/Å	convergence of each relaxation, as in `minimize`
`max_steps`	int32	5000		largest number of steps of each relaxation
`symmetry`	bool or `None`	`None`		take the stress under a strain that an operation makes from another as the image of its stress, with the symmetry operations of the crystal, found by spglib: `True`, `False`, or by default whenever spglib is installed

A cubic crystal takes 3 strains instead of 12: a stretch both ways and a shear.

A box periodic along two vectors in the x-y plane only makes the system a sheet, whose constants are in N/m.

Definition¶

The box and the atoms are strained by \(\pm\delta\) in each of the six components \(k\) of the strain in Voigt order xx, yy, zz, yz, xz, xy, the shear components as engineering strains (half on each off-diagonal entry), and the constants are central differences of the stress \(\sigma\) (positive under tension) in the same order:

\[ C_{jk} = \frac{\sigma_j(+\delta\, \mathbf e_k) - \sigma_j(-\delta\, \mathbf e_k)}{2\delta}. \]

With \(S = C^{-1}\), the moduli of a polycrystal are the Voigt and Reuss bounds and their mean, the Hill average (Hill (1952)):

\[ \begin{aligned} K_\mathrm V &= \tfrac19\left[C_{11} + C_{22} + C_{33} + 2(C_{12} + C_{13} + C_{23})\right], & K_\mathrm R &= \left[S_{11} + S_{22} + S_{33} + 2(S_{12} + S_{13} + S_{23})\right]^{-1}, \\ G_\mathrm V &= \tfrac1{15}\left[C_{11} + C_{22} + C_{33} - (C_{12} + C_{13} + C_{23}) + 3(C_{44} + C_{55} + C_{66})\right], & G_\mathrm R &= 15\left[4(S_{11} + S_{22} + S_{33}) - 4(S_{12} + S_{13} + S_{23}) + 3(S_{44} + S_{55} + S_{66})\right]^{-1}, \\ K_\mathrm H &= \tfrac12 (K_\mathrm V + K_\mathrm R), & G_\mathrm H &= \tfrac12 (G_\mathrm V + G_\mathrm R). \end{aligned} \]

Results¶

It returns an ElasticResult:

Attribute or method	Type	Meaning
`constants`	(6, 6) float64	elastic constants \(C\) in GPa, in Voigt order; read-only
`bulk_modulus(average="hill")`	float64	the bulk modulus in GPa: the `"voigt"`, `"reuss"` or `"hill"` average
`shear_modulus(average="hill")`	float64	the shear modulus in GPa, with the same averages
`errors`	(6, 6) float64 or `None`	standard errors in GPa, from molecular dynamics; None here

elastic_constants raises InvalidParameter for a box that is not periodic in every direction and is not that of a sheet, a non-positive strain, a calculator without the properties it needs, or a relaxation that does not converge within max_steps; bulk_modulus and shear_modulus raise it for an unknown average.