Elasticity Routines#

Code to calculate the elastic modulus tensor for athermal systems.

The elastic modulus tensor describes a material’s response to different boundary deformations. Specifically, for a small deformation given by a symmetric strain tensor e, the change in energy is

\[U / V^0 = U^0/V^0 + s^0_{ij} e_{ji} + (1/2) C_{ijkl} e_{ij} e_{kl} + \ldots{}\]

where \(V^0\) is the volume, \(U^0\) is the initial energy, \(s^0\) is the residual stress tensor of the undeformed system, and \(C\) is the elastic modulus tensor. \(C\) is a fourth-rank tensor of shape (dimension,dimension,dimension,dimension), with the following symmetries.

Minor symmetries:

\[C_{ijkl} = C_{jikl} = C_{ijlk}\]

Major symmetries:

\[C_{ijkl} = C_{lkij}\]

The minor symmetries are also reflected in the symmetric nature of stress and strain tensors:

\[s_{ij} = s_{ji} e_{ij} = e_{ji}\]

In general, there are 21 independent elastic constants in 3 dimension (6 in 2 dimensions). While systems with additional symmetries (e.g. isotropic, orthotropic, etc.) can be expressed with fewer constants, we do not assume any such additional symmetries.

At zero temperature, the response of every particle to a deformation can be calculated explicitly to linear order, enabling the exact calculation of the elastic modulus tensor without the need of any finite differences or any approximations.

Primary Functions#

jax_md.elasticity.athermal_moduli(energy_fn, tether_strength=1e-10, gradient_check=None, cg_tol=1e-07, check_convergence=False)[source]#

Setup calculation of elastic modulus tensor.

Parameters

  • energy_fn (Callable[…, Array]) – A function that computes the energy of the system. This function must take as an argument perturbation which perturbs the box shape. Any energy function constructed using smap or in energy.py with a standard space will satisfy this property.

  • tether_strength (float) – Scalar. Strength of the “tether” applied to each particle, which can be necessary to make the Hessian matrix non-singular. Solving for the non-affine response of each particle requires that the Hessian is positive definite. However, there can often be zero modes (eigenvectors of the Hessian with zero eigenvalue) that do not couple to the boundary, and therefore do not affect the elastic constants despite the zero eigenvalue. The most common example is the global translational modes. To solve for the non-affine response, we consider the “tethered Hessian” H + tether_strength * I, where I is the identity matrix and tether_strength is a small constant.

  • gradient_check (Optional[Array]) – None or scalar. If not None, a check will be performed to guarantee that the maximum component of the gradient is less than gradient_check. In other words, that jnp.amax(jnp.abs(grad(energy_fn)(R, box=box))) < gradient_check == True NOTE: JAX currently does not support proper runtime error handling. Therefore, if this check fails, the calculation will return an array of jnp.nan’s. It is the users responsibility, if they want to use this check, to then ensure that the returned array is not full of nans.

  • cg_tol (float) – scalar. Tolerance used when solving for the non-affine response.

  • check_convergence (bool) – bool. If true, calculate_EMT will return a boolean flag specifying if the cg solve routine converged to the desired tolerance. The default is False, but convergence checking is highly recommended especially when using 32-bit precision data.

Return: A function to calculate the elastic modulus tensor

Return type

Callable[…, Array]

Helper Functions#

jax_md.elasticity.tensor_to_mandel(T)[source]#

Convert a tensor to Mandel notation.

Mandel notation is a way to represent symmetric second-rank tensors and fourth-rank tensors with minor symmetries in a reduced form. Pairs of indices are combined as follows:

For tensors of shape (2,2) or (2,2,2,2): tensor indices Mandel indices 0,0 ————–> 0 1,1 ————–> 1 0,1 or 1,0 ——-> 2

For tensors of shape (3,3) or (3,3,3,3): tensor indices Mandel indices 0,0 ————–> 0 1,1 ————–> 1 2,2 ————–> 2 1,2 or 2,1 ——-> 3 0,2 or 2,0 ——-> 4 0,1 or 1,0 ——-> 5

If mandel_index(i,j) performs the above index mapping, then the input T and output M satisfy

M[mandel_index(i,j)] = T[i,j] * w

or

M[mandel_index(i,j), mandel_index(k,l)] = T[i,j,k,l] * w(i,j) * w(k,l)

where

w(i,j) = 1 if i==j

= sqrt(2) if i!=j

is a weight that is used to ensure proper contraction rules. Here (and only here) we do not assume major symmetries in fourth-rank tensors.

Parameters

T (Array) –

Array with 4 possible shapes:

  1. T.shape == (2,2) Convert a symmetric array of shape (2,2) to an array of shape (3,)

  2. T.shape == (3,3) Convert a symmetric array of shape (3,3) to an array of shape (6,)

  3. T.shape == (2,2,2,2) Convert a tensor of shape (2,2,2,2) with minor symmetries to an array of shape (3,3)

  4. T.shape == (3,3,3,3) Convert a tensor of shape (3,3,3,3) with minor symmetries to an array of shape (6,6)

Output: Array of shape (3,), (6,), (3,3), or (6,6)

see: https://sbrisard.github.io/janus/mandel.html (accessed 21 April, 2021)

Return type

Array

jax_md.elasticity.mandel_to_tensor(M)[source]#

Perform the inverse of M = tensor_to_mandel(T).

Parameters

M (Array) – Array of shape (3,), (6,), (3,3), or (6,6)

Output: Array of shape (2,2), (3,3), (2,2,2,2), or (3,3,3,3)

Return type

Array

jax_md.elasticity.extract_elements(C)[source]#

Convert an elastic modulus tensor into a list of unique elements.

In 2d, these are: cxxxx,cyyyy,cxyxy,cxxyy,cxxxy,cyyxy

In 3d, these are: cxxxx,cyyyy,czzzz,cyzyz,cxzxz,cxyxy,cyyzz,cxxzz,cxxyy,cxxyz,cxxxz,cxxxy, cyyyz,cyyxz,cyyxy,czzyz,czzxz,czzxy,cyzxz,cyzxy,cxzxy

Parameters

C (Array) – A previously calculated elastic modulus tensor represented as an array of shape (spatial_dimension,spatial_dimension,spatial_dimension, spatial_dimension), where spatial_dimension is either 2 or 3. C must satisfy both the major and minor symmetries, but this is not checked.

Return: a dict of the 6 (21) unique elastic constants in 2 (3) dimensions.

Return type

Dict

jax_md.elasticity.extract_isotropic_moduli(C)[source]#

Extract commonly used isotropic constants.

There are a number of important constants used to describe the linear elastic behavior of isotropic systems, including the bulk modulus, B, the shear modulus, G, the longitudinal modulus, M, the Young’s modulus, E, and the Poisson’s ratio, nu. This convenience function extracts them from an elastic modulus tensor C.

Angle averaged quantities: While these quantities are defined for isotropic systems, one can still define an “angle-averaged shear modulus”, for example, that averages over all possible shear deformations. This can be useful for systems that are statistically isotropic but where a particular realization is slightly anisotropic. The precise definitions are as follows:

First, we define the “response”, R, to a certain strain tensor e to be

R = 2 * (U / V^0 - U^0/V^0 - s^0_{ij} e_{ji}) = C_{ijkl} e_{ij} e_{kl}

Bulk modulus, B: This is the response to the rotationally invariant strain tensor:

e = (1/2) * ( 1 0 ) or e = (1/3) * ( 1 0 0 )
( 0 1 ) ( 0 1 0 )

( 0 0 1 )

Shear modulus, G: This is the response to the strain tensor:

e = (1/2) * ( 0 1 ) or e = (1/2) * ( 0 1 0 )
( 1 0 ) ( 1 0 0 )

( 0 0 0 )

averaged over all possible orientations. For perfectly isotropic systems, it should be equal to C[0,1,0,1].

Longitudinal modulus, M: This is the response to the strain tensor:

e = ( 1 0 ) or e = ( 1 0 0 )
( 0 0 ) ( 0 0 0 )

( 0 0 0 )

averaged over all possible orientations. For perfectly isotropic systems, it should be equal to C[0,0,0,0].

Young’s modulus, E: This is a measure of tensile stiffness and is calculated a using well- known expression in terms of B and G.

Poisson’s ratio: This is a measure of deformations in directions perpendicular to an applied load and is calculated a using well-known expression in terms of B and G.

Parameters

C (Array) – A previously calculated elastic modulus tensor represented as an array of shape (spatial_dimension,spatial_dimension,spatial_dimension, spatial_dimension), where spatial_dimension is either 2 or 3. C must satisfy both the major and minor symmetries, but this is not checked.

Return: a dictionary containing the elastic constants.

Return type

Dict