Elastic constants#
The theory in this section relates to the functions under the elastic
module.
Introduction#
Elastic constants measure the proportionality between strain and stress in a crystal, provided that the strain is not so large as to violate Hooke’s law.
—D. A. Papaconstantopolous
As noted by Papconstantopolous, the elastic constants may be considered
a three-dimensional generalisation of Hooke’s law, where put simply, the
restoring force
Discussion#
For the following we will be using the Voigt notation to simplify higher
order tensors to lower order ones. In particular reduction of the
fourth-order elastic stiffness tensor and second order stress/strain
tensors allows us to reduce the entire elasticity problem to a linear
equation. Derivation of the elastic tensor is quite lengthy and will not
be discussed here, instead we focus specifically on the stress and
strain matrices
For these works it is important to define what the directions
Resuming our discussion, regardless of orientation chosen, under Voigt notation these matrices are assumed symmetric. Therefore in the mapping, we make sure to take only the symmetric part, for stress we use:
or explicitly:
this conversion may be found in convert_stress_in_matrix_form_to_voigt_form()
,
whereas for strain (paying special attention the last three components) the mapping is:
or explicitly:
this conversion may be found in convert_strain_in_matrix_form_to_voigt_form()
.
We may then reverse the mapping to give our second-order tensors:
where the asymmetric components have been explicitly discarded, these functions are found under
convert_stress_in_voigt_form_to_matrix_form()
and
convert_strain_in_voigt_form_to_matrix_form()
respectively.
After mapping to Voigt notation, the most general stress-strain relation is thence:
where
Similarly, the total energy of the system is computed using a Taylor
expansion in
where
Computing elastic constants#
There are two main strategies for computing the elastic constants computationally. The first, known as the “Energy approach”, is based off of the energy changes one observes due to deforming the crystal using (2). The second, known as the “Stress theorem” is based on varying the strain of the crystal and monitoring the changes in the stress via (1). Both have their strengths and weaknesses as we will discuss below. To make the discussion more accessible, we demonstrate both methods using a model cubic system. Under the symmetries of a cubic system, the elastic tensor reduces to:
containing only 3 independent elastic constants.
Aside: Symmetry#
Certain methods such as ULICS [ULICS] use a set of dense vectors to sample
as much of phase space as possible by coupling as many degrees of
freedom as possible. This method though sensible in its goal of coupling
degrees of freedom falls short when the symmetry of the problem is taken
into account. Consider for instance a square lattice, the Irreducible
Brillouin Zone (IBZ) is reduced to 1/8th the size of the full Brillouin
Zone (BZ). If the symmetry is broken, and that all other properties
converge just as quickly, this would be at least an
Energy-based#
To compute properties efficiently using the energy approach we use the
basis sets detailed in the papers by Jamal et al. [Jamal2014], [Jamal2018],
where the authors optimized their vector choices to preserve the symmetry as
much as possible and (interestingly) simultaneously constraining them to
conserve volume where possible, where this second correction is most
important for high-pressure systems due to the nature of the energy expansion
in P and V. By preserving symmetry we can use fewer (
As energy converges much more quickly than other properties of Ab Initio simulations this method is numerically well-behaved but suffers from one major drawback; only one degree of freedom may be sampled as the energy method is a many-to-one operation. Nevertheless, it is both a useful method especially in the context of systems where the stress tensor (and the forces on atoms) cannot be computed. In such a case only the so-called “clamped” elastic constants can be computed as one cannot relax the atomic positions, or if say a method like GW, they are relaxed before applying the correction.
We break this method up into three parts, volumetric expansion, orthorhombic deformation and monoclinic deformation.
Volumetric expansion#
For volumetric expansion, the strain tensor is given by:
where
where (for cubic crystals)
where
we also use this method in our workflow to attain accurate estimates for the equilibrium lattice constant. To see why this is so helpful, first we take the Taylor series of (4) we arrive at:
which is the simplest expression that gives
and can be seen to be equivalent to (3), thus fitting the Birch-Murnaghan equation of state (which supports high/low pressure solids without or without nonlinear corrections) with least squares gives us a useful estimate where our energy scales act perturbatively, obtain a Bulk modulus value and get a reasonable lattice parameter estimate all in one go.
Orthorhombic deformation#
Orthorhombic deformation, as the name suggests, contracts one dimension while expanding another, it is given by:
where the (roughly) quadratic correction
where (for cubic crystals)
Monoclinic deformation#
The last deformation type is monoclinic, which unfortunately breaks most of the symmetry in a cubic lattice. It is given by:
where once again the
One then computes this using finite differences or linear regression.
which is unity for isotropic materials.
Stress-based#
As before, the eponymous stress-based method calculates the elastic constants by exploiting the many-to-many nature of the stress/strain relations. In order to efficiently and expediently calculate the properties of these crystals we use the Optimized High-Efficiency Strain-matrix Set (OHESS) [OHESS2020] strain sets which have been chosen to maximize both the symmetry and the number of degrees of freedom sampled in one evaluation.
In this example, the OHESS vector is able to sample all three elastic constants of a cubic crystal can be computed using only a single deformation type:
therefore by scaling
Data fitting method#
In this section we will be covering how we extract parameters from these fits with high confidence. We use an over-determined, least-squares approach using the Levenberg-Marquardt method of the lmfit python package to the entire dataset. This method provides several advantages over a simple least squares approach for each line.
Let’s say we deform a cubic crystal using a strain vector
where
where
Similarly, this extends to energy-based methods in a analogous fashion,
though we would need more functions to find the value of all three
elastic constants, we see that we have removed a degree of freedom when
it comes to the zero-strain solution as well as the zero-strain
components in
and similarly we gain energy information of the form:
Computing elastic constants using stress, stress, and energy least-squares fits respectively returns:
which we can see do not agree with each other. However, by constraining the reference strains to be equal, the elastic constants to be equal, and fitting using simultaneous least squares returns answers of:
where we can see that the elastic constants have changed minimally but
now their sum is consistent, the fig:diamondfits
,
we see that the reason for the discrepancy was likely due to the poor
fit of the energy-strain curve.

Stress-strain curve for

Stress-strain curve for

Energy-strain curve for given strain.#
Jamal, M., et al. "Elastic constants of cubic crystals." Computational Materials Science 95 (2014): 592-599.
Jamal, M., et al. "IRelast package." Journal of Alloys and Compounds 735 (2018): 569-579.
Liu, Zhong-Li. “High-efficiency calculation of elastic constants enhanced by the optimized strain-matrix sets." arXiv preprint: arXiv:2002.00005 (2020).
Yu, R., et al. "Calculations of single-crystal elastic constants made simple." Comput. Phys. Commun. 181 (2010): 671-675.