Introduction to atomman: ElasticConstants class

Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.

Disclaimers

1. Introduction

The ElasticConstants class represents the elastic constants of a crystal. The class methods focus on:

  • Allowing values to be set/retrieved in a number of different formats.

  • Correctly handling transformations to other Cartesian orientations.

Library Imports

[1]:

# Standard Python libraries
import datetime

# http://www.numpy.org/
import numpy as np
np.set_printoptions(precision=4, suppress=True)

# https://github.com/usnistgov/atomman
import atomman as am
import atomman.unitconvert as uc

# Show atomman version
print('atomman version =', am.__version__)

# Show date of Notebook execution
print('Notebook executed on', datetime.date.today())

atomman version = 1.4.10
Notebook executed on 2023-07-28

2. Elastic constants representations

The ElasticConstants object allows for various representations of the elastic constants to be retrieved:

  • Cij (6, 6) array of Voigt representation of elastic stiffness.

  • Sij (6, 6) array of Voigt representation of elastic compliance.

  • Cij9 (9, 9) array representation of elastic stiffness.

  • Cijkl (3, 3, 3, 3) array representation of elastic stiffness.

  • Sijkl (3, 3, 3, 3) array representation of elastic compliance.

2.1. Build demonstration (see Section #3 for more setting options)

[2]:

# 0 K values for Al
C11 = uc.set_in_units(1.143e12, 'dyn/cm^2')
C12 = uc.set_in_units(0.619e12, 'dyn/cm^2')
C44 = uc.set_in_units(0.316e12, 'dyn/cm^2')

C = am.ElasticConstants(C11=C11, C12=C12, C44=C44)

2.2. Show different representations

Cij is the 6x6 Voigt representation.

[3]:

print('Cij (GPa) ->')
print(uc.get_in_units(C.Cij, 'GPa'))

Cij (GPa) ->
[[114.3  61.9  61.9   0.    0.    0. ]
 [ 61.9 114.3  61.9   0.    0.    0. ]
 [ 61.9  61.9 114.3   0.    0.    0. ]
 [  0.    0.    0.   31.6   0.    0. ]
 [  0.    0.    0.    0.   31.6   0. ]
 [  0.    0.    0.    0.    0.   31.6]]

Cij9 is the full 9x9 representation.

[4]:

print('Cij9 (GPa) ->')
print(uc.get_in_units(C.Cij9, 'GPa'))

Cij9 (GPa) ->
[[114.3  61.9  61.9   0.    0.    0.    0.    0.    0. ]
 [ 61.9 114.3  61.9   0.    0.    0.    0.    0.    0. ]
 [ 61.9  61.9 114.3   0.    0.    0.    0.    0.    0. ]
 [  0.    0.    0.   31.6   0.    0.   31.6   0.    0. ]
 [  0.    0.    0.    0.   31.6   0.    0.   31.6   0. ]
 [  0.    0.    0.    0.    0.   31.6   0.    0.   31.6]
 [  0.    0.    0.   31.6   0.    0.   31.6   0.    0. ]
 [  0.    0.    0.    0.   31.6   0.    0.   31.6   0. ]
 [  0.    0.    0.    0.    0.   31.6   0.    0.   31.6]]

Cijkl is the full 3x3x3x3 representation.

[5]:

print('Cijkl (GPa) ->')
print(uc.get_in_units(C.Cijkl, 'GPa'))

Cijkl (GPa) ->
[[[[114.3   0.    0. ]
   [  0.   61.9   0. ]
   [  0.    0.   61.9]]

  [[  0.   31.6   0. ]
   [ 31.6   0.    0. ]
   [  0.    0.    0. ]]

  [[  0.    0.   31.6]
   [  0.    0.    0. ]
   [ 31.6   0.    0. ]]]


 [[[  0.   31.6   0. ]
   [ 31.6   0.    0. ]
   [  0.    0.    0. ]]

  [[ 61.9   0.    0. ]
   [  0.  114.3   0. ]
   [  0.    0.   61.9]]

  [[  0.    0.    0. ]
   [  0.    0.   31.6]
   [  0.   31.6   0. ]]]


 [[[  0.    0.   31.6]
   [  0.    0.    0. ]
   [ 31.6   0.    0. ]]

  [[  0.    0.    0. ]
   [  0.    0.   31.6]
   [  0.   31.6   0. ]]

  [[ 61.9   0.    0. ]
   [  0.   61.9   0. ]
   [  0.    0.  114.3]]]]

Sij is the Voigt 6x6 elastic compliances.

[6]:

print('Sij (1/GPa) ->')
print(uc.get_in_units(C.Sij, '1/GPa'))

Sij (1/GPa) ->
[[ 0.0141 -0.005  -0.005   0.      0.      0.    ]
 [-0.005   0.0141 -0.005   0.      0.      0.    ]
 [-0.005  -0.005   0.0141  0.      0.      0.    ]
 [ 0.      0.      0.      0.0316  0.      0.    ]
 [ 0.      0.      0.      0.      0.0316  0.    ]
 [ 0.      0.      0.      0.      0.      0.0316]]

Sijkl is the full 3x3x3x3 elastic compliances.

[7]:

print('Sijkl (1/GPa) ->')
print(uc.get_in_units(C.Sijkl, '1/GPa'))

Sijkl (1/GPa) ->
[[[[ 0.0141  0.      0.    ]
   [ 0.     -0.005   0.    ]
   [ 0.      0.     -0.005 ]]

  [[ 0.      0.0079  0.    ]
   [ 0.0079  0.      0.    ]
   [ 0.      0.      0.    ]]

  [[ 0.      0.      0.0079]
   [ 0.      0.      0.    ]
   [ 0.0079  0.      0.    ]]]


 [[[ 0.      0.0079  0.    ]
   [ 0.0079  0.      0.    ]
   [ 0.      0.      0.    ]]

  [[-0.005   0.      0.    ]
   [ 0.      0.0141  0.    ]
   [ 0.      0.     -0.005 ]]

  [[ 0.      0.      0.    ]
   [ 0.      0.      0.0079]
   [ 0.      0.0079  0.    ]]]


 [[[ 0.      0.      0.0079]
   [ 0.      0.      0.    ]
   [ 0.0079  0.      0.    ]]

  [[ 0.      0.      0.    ]
   [ 0.      0.      0.0079]
   [ 0.      0.0079  0.    ]]

  [[-0.005   0.      0.    ]
   [ 0.     -0.005   0.    ]
   [ 0.      0.      0.0141]]]]

3. Setting values

3.1. Setting using attributes

The values of the elastic constants can be initialized or set using any of the above representations.

  • During initialization, pass one of the following as the only parameter: Cij, Sij, Cij9, Cijkl, or Sijkl.

  • For an already initialized object, set any of the representations directly.

[8]:

Sijkl = C.Sijkl

# Initialize a new C with Sijkl values
newC = am.ElasticConstants(Sijkl = Sijkl)

# Show values to be the same
print(np.allclose(newC.Cij, C.Cij))

True

[9]:

# Initialize all zeros elastic constants array
newC = am.ElasticConstants()
print(newC)

# Set values using Cij9
newC.Cij9 = C.Cij9

# Show values to be the same
print(np.allclose(newC.Cij, C.Cij))

[[0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0.]]
True

3.2. Setting using crystal-specific constants

Alternatively, the elastic constants can be defined by supplying a full set of unique constants for a given crystal structure in the standard reference frame. During initialization, the specific form is inferred from the given parameters. After initialization, the functions corresponding to a specific crystal family can be called.

  • isotropic (two unique values required)

    • C11

    • C12

    • C44 (2*C44 = C11-C12)

    • M P-wave modulus (M = C11)

    • lambda Lame’s first parameter (lambda = C12)

    • mu Shear modulus (mu = C44)

    • E Young’s modulus

    • nu Poisson’s ratio

    • K Bulk modulus

  • cubic (all three values required)

    • C11

    • C12

    • C44

  • hexagonal (five unique values required)

    • C11

    • C12

    • C13

    • C33

    • C44

    • C66 (2*C66 = C11-C12)

  • tetragonal (six values required, one optional)

    • C11

    • C12

    • C13

    • C16 optional (C16=0 for some space groups)

    • C33

    • C44

    • C66

  • rhombohedral (six unique values required, one optional)

    • C11

    • C12

    • C13

    • C14

    • C15 optional (C15=0 for some space groups)

    • C33

    • C44

    • C66 (2*C66 = C11-C12)

  • orthorhombic (all nine values required)

    • C11

    • C12

    • C13

    • C22

    • C23

    • C33

    • C44

    • C55

    • C66

  • monoclinic (all thirteen values required)

    • C11

    • C12

    • C13

    • C15

    • C22

    • C23

    • C25

    • C33

    • C35

    • C44

    • C46

    • C55

    • C66

  • triclinic (all twenty-one values required)

    • all Cij where i <= j

[10]:

# Define dict with Zr HCP constants
Cdict = {}
Cdict['C11'] = uc.set_in_units(144.0, 'GPa')
Cdict['C12'] = uc.set_in_units(74.0, 'GPa')
Cdict['C13'] = uc.set_in_units(67.0, 'GPa')
Cdict['C33'] = uc.set_in_units(166.0, 'GPa')
Cdict['C44'] = uc.set_in_units(33.0, 'GPa')

# Initialize by passing dict key-values as parameters
C = am.ElasticConstants(**Cdict)

# Show that Cij array is properly constructed
print('Cij (GPa) ->')
print(uc.get_in_units(C.Cij, 'GPa'))
print()

# Show that 2 * C66 = C11 - C12
print('Does 2*C66 = C11-C12?')
print('C11 - C12 (GPa) =', uc.get_in_units(C.Cij[0,0] - C.Cij[0,1], 'GPa'))
print('2 * C66 (GPa) =  ', uc.get_in_units(2 * C.Cij[5,5], 'GPa'))

Cij (GPa) ->
[[144.  74.  67.   0.   0.   0.]
 [ 74. 144.  67.   0.   0.   0.]
 [ 67.  67. 166.   0.   0.   0.]
 [  0.   0.   0.  33.   0.   0.]
 [  0.   0.   0.   0.  33.   0.]
 [  0.   0.   0.   0.   0.  35.]]

Does 2*C66 = C11-C12?
C11 - C12 (GPa) = 70.0
2 * C66 (GPa) =   70.0

[11]:

# Define dict with Bi2Te3 rhombohedral constants
Cdict = {}
Cdict['C11'] = uc.set_in_units(7.436, '10^11 * dyn/cm^2')
Cdict['C12'] = uc.set_in_units(2.619, '10^11 * dyn/cm^2')
Cdict['C33'] = uc.set_in_units(5.160, '10^11 * dyn/cm^2')
Cdict['C44'] = uc.set_in_units(3.135, '10^11 * dyn/cm^2')
Cdict['C13'] = uc.set_in_units(2.917, '10^11 * dyn/cm^2')
Cdict['C14'] = uc.set_in_units(1.541, '10^11 * dyn/cm^2')

# Set values to existing C using the rhombohedral function
C.rhombohedral(**Cdict)

# Show that Cij array is properly constructed
print('Cij (GPa) ->')
print(uc.get_in_units(C.Cij, 'GPa'))
print()

# Show that 2 * C66 = C11 - C12
print('Does 2*C66 = C11-C12?')
print('C11 - C12 (GPa) =', uc.get_in_units(C.Cij[0,0] - C.Cij[0,1], 'GPa'))
print('2 * C66 (GPa) =  ', uc.get_in_units(2 * C.Cij[5,5], 'GPa'))

Cij (GPa) ->
[[ 74.36   26.19   29.17   15.41    0.      0.   ]
 [ 26.19   74.36   29.17  -15.41    0.      0.   ]
 [ 29.17   29.17   51.6     0.      0.      0.   ]
 [ 15.41  -15.41    0.     31.35    0.      0.   ]
 [  0.      0.      0.      0.     31.35   15.41 ]
 [  0.      0.      0.      0.     15.41   24.085]]

Does 2*C66 = C11-C12?
C11 - C12 (GPa) = 48.170000000000016
2 * C66 (GPa) =   48.170000000000016

4. ElasticConstants.tranform()

The transform() method is included for convenience allowing for the elastic constants to be transformed to a different set of axes.

Parameters

  • axes (3, 3) array giving three right-handed orthogonal vectors to use for transforming.

  • tol optional relative tolerance to use in identifying near-zero terms.

Returns a new ElasticConstants object with constants transformed to the new axes.

[12]:

# Set C back to a cubic system (Cu this time)
Cdict = {}
Cdict['C11'] = uc.set_in_units(1.762, '10^12 * dyn/cm^2')
Cdict['C12'] = uc.set_in_units(1.249, '10^12 * dyn/cm^2')
Cdict['C44'] = uc.set_in_units(0.818, '10^12 * dyn/cm^2')

C.cubic(**Cdict)
print('Cij (GPa) ->')
print(uc.get_in_units(C.Cij, 'GPa'))
print()


# Transform by a 45 degree rotation around z-axis
axes = [[ 1, 1, 0],
        [-1, 1, 0],
        [ 0, 0, 1]]
newC = C.transform(axes)

print('After transforming, Cij (GPa) ->')
print(uc.get_in_units(newC.Cij, 'GPa'))

Cij (GPa) ->
[[176.2 124.9 124.9   0.    0.    0. ]
 [124.9 176.2 124.9   0.    0.    0. ]
 [124.9 124.9 176.2   0.    0.    0. ]
 [  0.    0.    0.   81.8   0.    0. ]
 [  0.    0.    0.    0.   81.8   0. ]
 [  0.    0.    0.    0.    0.   81.8]]

After transforming, Cij (GPa) ->
[[232.35  68.75 124.9    0.     0.     0.  ]
 [ 68.75 232.35 124.9    0.     0.     0.  ]
 [124.9  124.9  176.2    0.     0.     0.  ]
 [  0.     0.     0.    81.8    0.     0.  ]
 [  0.     0.     0.     0.    81.8    0.  ]
 [  0.     0.     0.     0.     0.    25.65]]

5. Shear and bulk modulus estimates

5.1. ElasticConstants.bulk()

Three style options of estimates are available:

  • ‘Voigt’ Voigt estimate.

\[K_{Voigt} = \frac{ \left(C_{11} + C_{22} + C_{33} \right) + 2 \left(C_{12} + C_{13} + C_{23} \right) }{9}\]

  • ‘Reuss’ Reuss estimate.

\[K_{Reuss} = \frac{1}{ \left( S_{11} + S_{22} + S_{33} \right) + 2 \left(S_{12} + S_{13} + S_{23} \right) }\]

  • ‘Hill’ Hill estimate (default).

\[K_{Hill} = \frac{K_{Reuss} + K_{Voigt}}{2}\]
[13]:

print('Voigt bulk modulus estimate =', uc.get_in_units(C.bulk('Voigt'), 'GPa'), 'GPa')
print('Reuss bulk modulus estimate =', uc.get_in_units(C.bulk('Reuss'), 'GPa'), 'GPa')
print('Hill bulk modulus estimate = ', uc.get_in_units(C.bulk(), 'GPa'), 'GPa')

Voigt bulk modulus estimate = 142.0 GPa
Reuss bulk modulus estimate = 141.9999999999999 GPa
Hill bulk modulus estimate =  141.99999999999994 GPa

5.2. ElasticConstants.shear()

Three style options of estimates are available:

  • ‘Voigt’ Voigt estimate.

\[\mu_{Voigt} = \frac{ \left( C_{11} + C_{22} + C_{33} \right) - \left( C_{12} + C_{23} + C_{13} \right) + 3 \left( C_{44} + C_{55} + C_{66} \right) }{15}\]

  • ‘Reuss’ Reuss estimate.

\[\mu_{Reuss} = \frac{15}{4 \left( S_{11} + S_{22} + S_{33} \right) - 4 \left( S_{12} + S_{23} + S_{13} \right) + 3 \left( S_{44} + S_{55} + S_{66} \right)}\]

  • ‘Hill’ Hill estimate (default).

\[\mu_{Hill} = \frac{\mu_{Reuss} + \mu_{Voigt}}{2}\]
[14]:

print('Voigt shear modulus estimate =', uc.get_in_units(C.shear('Voigt'), 'GPa'), 'GPa')
print('Reuss shear modulus estimate =', uc.get_in_units(C.shear('Reuss'), 'GPa'), 'GPa')
print('Hill shear modulus estimate = ', uc.get_in_units(C.shear(), 'GPa'), 'GPa')

Voigt shear modulus estimate = 59.34000000000001 GPa
Reuss shear modulus estimate = 43.611930991477855 GPa
Hill shear modulus estimate =  51.47596549573894 GPa

6. Crystal system normalized values

The elastic constants tensor should have certain components that are equal or dependent based on the crystal symmetry of the system. However, calculations of the elastic constants for an atomic system may show some variability across the measured values of these symmetrically-dependent components. The methods listed here help to handle the symmetry components.

6.1. ElasticConstants.normalized_as()

Returns a new ElasticConstants object where values of the current are averaged or zeroed out according to a standard crystal system setting.

NOTE: no validation checks are made to evaluate whether such normalizations should be done! That is left up to you (compare values before and after normalization).

Parameters

  • crystal_system (str) Indicates the crystal system representation to use when building a data model.

Returns

  • (atomman.ElasticConstants) The elastic constants normalized according to the crystal system symmetries.

[15]:

# Create a new ElasticConstants object with some variability
scale = uc.set_in_units(1e-2, 'GPa')
Cij_with_v = C.Cij + scale * np.random.rand(6,6) - scale/2
for i in range(6):
    for j in range(i, 6):
        Cij_with_v[i,j] = Cij_with_v[j, i] = (Cij_with_v[i,j]+Cij_with_v[j, i])/2
C_with_v = am.ElasticConstants(Cij=Cij_with_v)

print('C_with_v:')
print(uc.get_in_units(C_with_v.Cij, 'GPa'))
print()

# Normalize to cubic
C_norm = C_with_v.normalized_as('cubic')
print('C_norm')
print(uc.get_in_units(C_norm.Cij, 'GPa'))

C_with_v:
[[176.2003 124.898  124.9001   0.0021   0.0018   0.0014]
 [124.898  176.204  124.9001   0.0037   0.0005  -0.0002]
 [124.9001 124.9001 176.2046   0.0019  -0.0016   0.0012]
 [  0.0021   0.0037   0.0019  81.7951  -0.0008  -0.0001]
 [  0.0018   0.0005  -0.0016  -0.0008  81.7992  -0.002 ]
 [  0.0014  -0.0002   0.0012  -0.0001  -0.002   81.7966]]

C_norm
[[176.203  124.8994 124.8994   0.       0.       0.    ]
 [124.8994 176.203  124.8994   0.       0.       0.    ]
 [124.8994 124.8994 176.203    0.       0.       0.    ]
 [  0.       0.       0.      81.797    0.       0.    ]
 [  0.       0.       0.       0.      81.797    0.    ]
 [  0.       0.       0.       0.       0.      81.797 ]]

6.2. is_normal()

Checks if current elastic constants agree with values normalized to a specified crystal family (within tolerances).

Parameters

  • crystal_system (str) Indicates the crystal system representation to use when building a data model.

  • atol (float, optional) Absolute tolerance to use. Default value is 1e-4.

  • rtol (float, optional) Relative tolerance to use. Default value is 1e-4.

Returns

  • (bool) True if all Cij match within the tolerances, false otherwise.

[16]:

print("C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4) ->", C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4))
print("C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7) ->", C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7))
print("C_with_v.is_normal('hexagonal') ->", C_with_v.is_normal('hexagonal'))
print()

print("C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4) ->", C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4))
print("C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7) ->", C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7))
print("C_norm.is_normal('hexagonal') ->", C_norm.is_normal('hexagonal'))

C_with_v.is_normal('cubic', rtol=1e-4, atol=1e-4) -> True
C_with_v.is_normal('cubic', rtol=1e-7, atol=1e-7) -> False
C_with_v.is_normal('hexagonal') -> False

C_norm.is_normal('cubic', rtol=1e-4, atol=1e-4) -> True
C_norm.is_normal('cubic', rtol=1e-7, atol=1e-7) -> True
C_norm.is_normal('hexagonal') -> False

7. ElasticConstants.model()

The elastic constants can also be saved/retrieved as a JSON/XML data model. This is useful as it captures not only the elastic constant values but also the associated units.

Parameters

  • model (DataModelDict, string, or file-like object, optional) Data model containing exactly one ‘elastic-constants’ branch to read.

  • unit (str, optional) Units or pressure to save values in when building a data model. Default value is None (no conversion).

  • crystal_system (str, optional) Indicates the crystal system representation to use when building a data model. Default value is ‘triclinic’ (save all values in Cij).

Returns

  • (DataModelDict) If model is not given as a parameter.

[17]:

# Generate data model based on C
model = C.model(unit='GPa', crystal_system='cubic')

# Show json version of model
print(model.json(indent=4))

{
    "elastic-constants": {
        "Cij": {
            "value": [
                176.20000000000002,
                124.9,
                124.9,
                0.0,
                0.0,
                0.0,
                124.9,
                176.20000000000002,
                124.9,
                0.0,
                0.0,
                0.0,
                124.9,
                124.9,
                176.20000000000002,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                81.8,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                81.8,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                0.0,
                81.8
            ],
            "shape": [
                6,
                6
            ],
            "unit": "GPa"
        }
    }
}

[18]:

# If crystal_system is not given, model will contain all 21 triclinic constants
model = C.model(unit='GPa')
print(model.json())

{"elastic-constants": {"Cij": {"value": [176.20000000000002, 124.9, 124.9, 0.0, 0.0, 0.0, 124.9, 176.20000000000002, 124.9, 0.0, 0.0, 0.0, 124.9, 124.9, 176.20000000000002, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 81.8], "shape": [6, 6], "unit": "GPa"}}}

The model can also be read in to an ElasticConstants object either during initialization or using the model() method.

[19]:

C = am.ElasticConstants(model=model)
print(uc.get_in_units(C.Cij, 'GPa'))

[[176.2 124.9 124.9   0.    0.    0. ]
 [124.9 176.2 124.9   0.    0.    0. ]
 [124.9 124.9 176.2   0.    0.    0. ]
 [  0.    0.    0.   81.8   0.    0. ]
 [  0.    0.    0.    0.   81.8   0. ]
 [  0.    0.    0.    0.    0.   81.8]]

[20]:

newC.model(model=model)
print(uc.get_in_units(newC.Cij, 'GPa'))

[[176.2 124.9 124.9   0.    0.    0. ]
 [124.9 176.2 124.9   0.    0.    0. ]
 [124.9 124.9 176.2   0.    0.    0. ]
 [  0.    0.    0.   81.8   0.    0. ]
 [  0.    0.    0.    0.   81.8   0. ]
 [  0.    0.    0.    0.    0.   81.8]]