elastic_constants_static - Methodology and code

Python imports

[1]:
# Standard library imports
from pathlib import Path
import shutil
import datetime
from copy import deepcopy
from math import floor
from typing import Optional, Tuple

# http://www.numpy.org/
import numpy as np

# https://ipython.org/
from IPython.display import display, Markdown

# https://github.com/usnistgov/atomman
import atomman as am
import atomman.lammps as lmp
import atomman.unitconvert as uc
from atomman.tools import filltemplate

# https://github.com/usnistgov/iprPy
import iprPy
from iprPy.tools import read_calc_file

print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)
Notebook last executed on 2022-03-10 using iprPy version 0.11.2

1. Load calculation and view description

1.1. Load the calculation

[2]:
# Load the calculation being demoed
calculation = iprPy.load_calculation('elastic_constants_static')

1.2. Display calculation description and theory

[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))

elastic_constants_static calculation style

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

Introduction

The elastic_constants_static calculation style computes the elastic constants, \(C_{ij}\), for a system by applying small strains and performing static energy minimizations of the initial and strained configurations. Three estimates of the elastic constants are returned: one for applying positive strains, one for applying negative strains, and a normalized estimate that averages the ± strains and the symmetric components of the \(C_{ij}\) tensor.

Version notes
  • 2018-07-09: Notebook added.

  • 2019-07-30: Description updated and small changes due to iprPy version.

  • 2020-05-22: Version 0.10 update - potentials now loaded from database.

  • 2020-09-22: Setup and parameter definition streamlined.

  • 2022-03-11: Notebook updated to reflect version 0.11.

Additional dependencies
Disclaimers
  • NIST disclaimers

  • Unlike the previous LAMMPS_ELASTIC calculation, this calculation does not perform a box relaxation on the system prior to evaluating the elastic constants. This allows for the static elastic constants to be evaluated for systems that are relaxed to different pressures.

  • The elastic constants are estimated using small strains. Depending on the potential, the values for the elastic constants may vary with the size of the strain. This can come about either if the strain exceeds the linear elastic regime.

  • Some classical interatomic potentials have discontinuities in the fourth derivative of the energy function with respect to position. If the strained states straddle one of these discontinuities the resulting static elastic constants values will be nonsense.

Method and Theory

The calculation method used here for computing elastic constants is based on the method used in the ELASTIC demonstration script created by Steve Plimpton and distributed with LAMMPS.

The math in this section uses Voigt notation, where indicies i,j correspond to 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, and 6=xy, and x, y and z are orthogonal box vectors.

A LAMMPS simulation performs thirteen energy/force minimizations

  • One for relaxing the initial system.

  • Twelve for relaxing systems in which a small strain of magnitude \(\Delta \epsilon\) is applied to the system in both the positive and negative directions of the six Voigt strain components, \(\pm \Delta \epsilon_{i}\).

The system virial pressures, \(P_{i}\), are recorded for each of the thirteen relaxed configurations. Two estimates for the \(C_{ij}\) matrix for the system are obtained as

\[C_{ij}^+ = - \frac{P_i(\Delta \epsilon_j) - P_i(0)}{\Delta \epsilon},\]
\[C_{ij}^- = - \frac{P_i(0) - P_i(-\Delta \epsilon_j)}{\Delta \epsilon}.\]

The negative out front comes from the fact that the system-wide stress state is \(\sigma_i = -P_i\). A normalized, average estimate is also obtained by averaging the positive and negative strain estimates, as well as the symmetric components of the tensor

\[C_{ij} = \frac{C_{ij}^+ + C_{ij}^- + C_{ji}^+ + C_{ji}^-}{4}.\]

2. Define calculation functions and generate files

This section defines the calculation functions and associated resource files exactly as they exist inside the iprPy package. This allows for the code used to be directly visible and modifiable by anyone looking to see how it works.

2.1. elastic_constants_static()

This is the primary function for the calculation. The version of this function built in iprPy can be accessed by calling the calc() method of an object of the associated calculation class.

[4]:
def elastic_constants_static(lammps_command: str,
                             system: am.System,
                             potential: lmp.Potential,
                             mpi_command: Optional[str] = None,
                             strainrange: float = 1e-6,
                             etol: float = 0.0,
                             ftol: float = 0.0,
                             maxiter: int = 10000,
                             maxeval: int = 100000,
                             dmax: float = uc.set_in_units(0.01, 'angstrom')) -> dict:
    """
    Repeatedly runs the ELASTIC example distributed with LAMMPS until box
    dimensions converge within a tolerance.

    Parameters
    ----------
    lammps_command :str
        Command for running LAMMPS.
    system : atomman.System
        The system to perform the calculation on.
    potential : atomman.lammps.Potential
        The LAMMPS implemented potential to use.
    mpi_command : str, optional
        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS
        will run serially.
    strainrange : float, optional
        The small strain value to apply when calculating the elastic
        constants (default is 1e-6).
    etol : float, optional
        The energy tolerance for the structure minimization. This value is
        unitless. (Default is 0.0).
    ftol : float, optional
        The force tolerance for the structure minimization. This value is in
        units of force. (Default is 0.0).
    maxiter : int, optional
        The maximum number of minimization iterations to use (default is 10000).
    maxeval : int, optional
        The maximum number of minimization evaluations to use (default is
        100000).
    dmax : float, optional
        The maximum distance in length units that any atom is allowed to relax
        in any direction during a single minimization iteration (default is
        0.01 Angstroms).

    Returns
    -------
    dict
        Dictionary of results consisting of keys:

        - **'raw_Cij_negative'** (*numpy.ndarray*) - The values of Cij obtained
          from only the negative strains.
        - **'raw_Cij_positive'** (*numpy.ndarray*) - The values of Cij obtained
          from only the positive strains.
        - **'C'** (*atomman.ElasticConstants*) - The computed elastic constants
          obtained from averaging the negative and positive strain values.
    """

    # Convert hexagonal cells to orthorhombic to avoid LAMMPS tilt issues
    if am.tools.ishexagonal(system.box):
        system = system.rotate([[2,-1,-1,0], [0, 1, -1, 0], [0,0,0,1]])

    # Get lammps units
    lammps_units = lmp.style.unit(potential.units)

    # Get lammps version date
    lammps_date = lmp.checkversion(lammps_command)['date']

    # Define lammps variables
    lammps_variables = {}
    system_info = system.dump('atom_data', f='init.dat',
                              potential=potential)
    lammps_variables['atomman_system_pair_info'] = system_info
    lammps_variables['restart_commands'] = restart_commands(potential, system.symbols)
    lammps_variables['strainrange'] = strainrange
    lammps_variables['etol'] = etol
    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])
    lammps_variables['maxiter'] = maxiter
    lammps_variables['maxeval'] = maxeval
    lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])

    # Fill in template files
    lammps_script = 'cij.in'
    template = read_calc_file('iprPy.calculation.elastic_constants_static',
                              'cij.template')
    with open(lammps_script, 'w') as f:
        f.write(filltemplate(template, lammps_variables, '<', '>'))

    # Run LAMMPS
    output = lmp.run(lammps_command, script_name=lammps_script,
                     mpi_command=mpi_command)

    # Pull out initial state
    thermo = output.simulations[0]['thermo']
    pxx0 = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
    pyy0 = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
    pzz0 = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
    pyz0 = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
    pxz0 = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
    pxy0 = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

    # Negative strains
    cij_n = np.empty((6,6))
    for i in range(6):
        j = 1 + i * 2
        # Pull out strained state
        thermo = output.simulations[j]['thermo']
        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

        # Calculate cij_n using stress changes
        cij_n[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,
                             pyz - pyz0, pxz - pxz0, pxy - pxy0]) / strainrange

    # Positive strains
    cij_p = np.empty((6,6))
    for i in range(6):
        j = 2 + i * 2
        # Pull out strained state
        thermo = output.simulations[j]['thermo']
        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

        # Calculate cij_p using stress changes
        cij_p[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,
                              pyz - pyz0, pxz - pxz0, pxy - pxy0]) / -strainrange

    # Average symmetric values
    cij = (cij_n + cij_p) / 2
    for i in range(6):
        for j in range(i):
            cij[i,j] = cij[j,i] = (cij[i,j] + cij[j,i]) / 2

    # Define results_dict
    results_dict = {}
    results_dict['raw_Cij_negative'] = cij_n
    results_dict['raw_Cij_positive'] = cij_p
    results_dict['C'] = am.ElasticConstants(Cij=cij)

    return results_dict

2.2. restart_commands

[5]:
def restart_commands(potential: lmp.Potential,
                     symbols: list) -> str:
    """
    Command lines to restart calculation from the initial relaxation

    Parameters
    ----------
    potential : lmp.Potential
        The interatomic potential.
    symbols : list
        The list of symbol models associated with the interatomic potential.
    """

    if potential.pair_style == 'kim':
        pair_info = potential.pair_info(symbols)
        commands = '\n'.join([
            pair_info.split('\n')[0],
            'read_restart initial.restart',
        ])
        commands += '\n' + '\n'.join(pair_info.split('\n')[1:])

    else:
        commands = '\n'.join([
            'read_restart initial.restart',
            potential.pair_info(symbols),
        ])

    commands += '\n'.join([
        '',
        '# Setup minimization style',
        'min_modify dmax ${dmax}',
        '',
        '# Setup output',
        'thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe',
        'thermo_modify format float %.13e'])

    return commands

2.3. cij.template file

[6]:
with open('cij.template', 'w') as f:
    f.write("""# Performs simulations to statically evaluate elastic constants using small strains
# Based on the LAMMPS_ELASTIC script by Aidan Thompson (Sandia, athomps@sandia.gov)

box tilt large

<atomman_system_pair_info>

change_box all triclinic

# Specify strain
variable strain equal <strainrange>

# Define minimization parameters
variable etol equal <etol>
variable ftol equal <ftol>
variable maxiter equal <maxiter>
variable maxeval equal <maxeval>
variable dmax equal <dmax>

# Specify variables of the initial configuration's dimensions
variable lx0 equal $(lx)
variable ly0 equal $(ly)
variable lz0 equal $(lz)

# Specify the thermo properties to calculate
variable peatom equal pe/atoms

# Read in potential and thermo information
# Setup minimization style
min_modify dmax ${dmax}

# Setup output
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e

# Relax initial configuration and save as restart
minimize ${etol} ${ftol} ${maxiter} ${maxeval}
write_restart initial.restart

# Apply -xx strain
clear
<restart_commands>

variable delta equal -${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xx strain
clear
<restart_commands>

variable delta equal ${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -yy strain
clear
<restart_commands>

variable delta equal -${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +yy strain
clear
<restart_commands>

variable delta equal ${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -zz strain
clear
<restart_commands>

variable delta equal -${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +zz strain
clear
<restart_commands>

variable delta equal ${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -yz strain
clear
<restart_commands>

variable delta equal -${strain}*${lz0}
change_box all yz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +yz strain
clear
<restart_commands>

variable delta equal ${strain}*${lz0}
change_box all yz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -xz strain
clear
<restart_commands>

variable delta equal -${strain}*${lz0}
change_box all xz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xz strain
clear
<restart_commands>

variable delta equal ${strain}*${lz0}
change_box all xz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -xy strain
clear
<restart_commands>

variable delta equal -${strain}*${ly0}
change_box all xy delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xy strain
clear
<restart_commands>

variable delta equal ${strain}*${ly0}
change_box all xy delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}""")

3. Specify input parameters

3.1. System-specific paths

  • lammps_command is the LAMMPS command to use (required).

  • mpi_command MPI command for running LAMMPS in parallel. A value of None will run simulations serially.

[7]:
lammps_command = 'lmp_serial'
mpi_command = None

3.2. Interatomic potential

  • potential_name gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.

  • potential is an atomman.lammps.Potential object (required).

[8]:
potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'

# Retrieve potential and parameter file(s) using atomman
potential = am.load_lammps_potential(id=potential_name, getfiles=True)

3.3. Initial unit cell system

  • ucell is an atomman.System representing a fundamental unit cell of the system (required). Here, this is generated using the load parameters and symbols.

[9]:
# Create ucell by loading prototype record
ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc')

print(ucell)
avect =  [ 3.520,  0.000,  0.000]
bvect =  [ 0.000,  3.520,  0.000]
cvect =  [ 0.000,  0.000,  3.520]
origin = [ 0.000,  0.000,  0.000]
natoms = 4
natypes = 1
symbols = ('Ni',)
pbc = [ True  True  True]
per-atom properties = ['atype', 'pos']
     id   atype  pos[0]  pos[1]  pos[2]
      0       1   0.000   0.000   0.000
      1       1   0.000   1.760   1.760
      2       1   1.760   0.000   1.760
      3       1   1.760   1.760   0.000

3.4. System modifications

  • sizemults list of three integers specifying how many times the ucell vectors of \(a\), \(b\) and \(c\) are replicated in creating system.

  • system is an atomman.System to perform the scan on (required).

[10]:
sizemults = [3, 3, 3]

# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 108

3.5. Calculation-specific parameters

  • strainrange specifies the \(\Delta \epsilon\) strain range to use in estimating \(C_{ij}\).

  • energytolerance is the energy tolerance to use during the minimizations. This is unitless.

  • forcetolerance is the force tolerance to use during the minimizations. This is in energy/length units.

  • maxiterations is the maximum number of minimization iterations to use.

  • maxevaluations is the maximum number of minimization evaluations to use.

  • maxatommotion is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units.

[11]:
strainrange = 1e-7
energytolerance = 1e-8
forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')
maxiterations = 10000
maxevaluations = 100000
maxatommotion = uc.set_in_units(0.01, 'angstrom')

4. Run calculation and view results

4.1. Run calculation

All primary calculation method functions take a series of inputs and return a dictionary of outputs.

[12]:
results_dict = elastic_constants_static(lammps_command, system, potential,
                                        mpi_command = mpi_command,
                                        strainrange = strainrange,
                                        etol = energytolerance,
                                        ftol = forcetolerance,
                                        maxiter = maxiterations,
                                        maxeval = maxevaluations,
                                        dmax = maxatommotion)
print(results_dict.keys())
dict_keys(['raw_Cij_negative', 'raw_Cij_positive', 'C'])

4.2. Report results

Values returned in the results_dict:

  • ‘raw_Cij_negative’ (numpy.ndarray) - The values of Cij obtained from only the negative strains.

  • ‘raw_Cij_positive’ (numpy.ndarray) - The values of Cij obtained from only the positive strains.

  • ‘C’ (atomman.ElasticConstants) - The computed elastic constants obtained from averaging the negative and positive strain values.

[13]:
pressure_unit = 'GPa'

print('Raw Cij for negative strains ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['raw_Cij_negative'], pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
print()

print('Raw Cij for positive strains ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['raw_Cij_positive'], pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
Raw Cij for negative strains (GPa) =
[ 247.8622  147.8285  147.8285    0.0000    0.0000   -0.0000]
[ 147.8285  247.8622  147.8285    0.0000    0.0000   -0.0000]
[ 147.8285  147.8285  247.8622   -0.0000    0.0000   -0.0000]
[  -0.0000    0.0001    0.0001  124.8381   -0.0000   -0.0000]
[   0.0001   -0.0000    0.0001   -0.0000  124.8381   -0.0000]
[   0.0001    0.0001   -0.0000   -0.0000   -0.0000  124.8381]

Raw Cij for positive strains (GPa) =
[ 247.8625  147.8283  147.8283   -0.0000    0.0000    0.0000]
[ 147.8283  247.8625  147.8283    0.0000   -0.0000    0.0000]
[ 147.8283  147.8283  247.8625    0.0000   -0.0000    0.0000]
[   0.0000   -0.0001   -0.0001  124.8381    0.0000    0.0000]
[  -0.0001    0.0000   -0.0001    0.0000  124.8381    0.0000]
[  -0.0001   -0.0001    0.0000    0.0000    0.0000  124.8381]
[14]:
print('Cij ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['C'].Cij, pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
Cij (GPa) =
[ 247.8624  147.8284  147.8284    0.0000    0.0000    0.0000]
[ 147.8284  247.8624  147.8284    0.0000    0.0000    0.0000]
[ 147.8284  147.8284  247.8624    0.0000    0.0000    0.0000]
[   0.0000    0.0000    0.0000  124.8381    0.0000    0.0000]
[   0.0000    0.0000    0.0000    0.0000  124.8381    0.0000]
[   0.0000    0.0000    0.0000    0.0000    0.0000  124.8381]
[15]:
family = am.tools.identifyfamily(ucell.box)
C = results_dict['C']

if not C.is_normal(family):
    print("Cij not consistent with ucell's box")

else:
    norm_C = C.normalized_as(family)

    if family == 'cubic':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
    elif family == 'hexagonal':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C13 =', uc.get_in_units(norm_C.Cij[0,2], 'GPa'), 'GPa')
        print('C33 =', uc.get_in_units(norm_C.Cij[2,2], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
    elif family == 'tetragonal':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C13 =', uc.get_in_units(norm_C.Cij[0,2], 'GPa'), 'GPa')
        print('C16 =', uc.get_in_units(norm_C.Cij[0,5], 'GPa'), 'GPa')
        print('C33 =', uc.get_in_units(norm_C.Cij[2,2], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
        print('C66 =', uc.get_in_units(norm_C.Cij[5,5], 'GPa'), 'GPa')
    elif family == 'rhombohedral':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C13 =', uc.get_in_units(norm_C.Cij[0,2], 'GPa'), 'GPa')
        print('C14 =', uc.get_in_units(norm_C.Cij[0,3], 'GPa'), 'GPa')
        print('C15 =', uc.get_in_units(norm_C.Cij[0,4], 'GPa'), 'GPa')
        print('C33 =', uc.get_in_units(norm_C.Cij[2,2], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
    elif family == 'orthorhombic':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C13 =', uc.get_in_units(norm_C.Cij[0,2], 'GPa'), 'GPa')
        print('C22 =', uc.get_in_units(norm_C.Cij[1,1], 'GPa'), 'GPa')
        print('C23 =', uc.get_in_units(norm_C.Cij[1,2], 'GPa'), 'GPa')
        print('C33 =', uc.get_in_units(norm_C.Cij[2,2], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
        print('C55 =', uc.get_in_units(norm_C.Cij[4,4], 'GPa'), 'GPa')
        print('C66 =', uc.get_in_units(norm_C.Cij[5,5], 'GPa'), 'GPa')
    elif family == 'monoclinic':
        print('C11 =', uc.get_in_units(norm_C.Cij[0,0], 'GPa'), 'GPa')
        print('C12 =', uc.get_in_units(norm_C.Cij[0,1], 'GPa'), 'GPa')
        print('C13 =', uc.get_in_units(norm_C.Cij[0,2], 'GPa'), 'GPa')
        print('C15 =', uc.get_in_units(norm_C.Cij[0,4], 'GPa'), 'GPa')
        print('C22 =', uc.get_in_units(norm_C.Cij[1,1], 'GPa'), 'GPa')
        print('C23 =', uc.get_in_units(norm_C.Cij[1,2], 'GPa'), 'GPa')
        print('C25 =', uc.get_in_units(norm_C.Cij[1,4], 'GPa'), 'GPa')
        print('C33 =', uc.get_in_units(norm_C.Cij[2,2], 'GPa'), 'GPa')
        print('C35 =', uc.get_in_units(norm_C.Cij[2,4], 'GPa'), 'GPa')
        print('C44 =', uc.get_in_units(norm_C.Cij[3,3], 'GPa'), 'GPa')
        print('C46 =', uc.get_in_units(norm_C.Cij[3,5], 'GPa'), 'GPa')
        print('C55 =', uc.get_in_units(norm_C.Cij[4,4], 'GPa'), 'GPa')
        print('C66 =', uc.get_in_units(norm_C.Cij[5,5], 'GPa'), 'GPa')
    else:
        print('system is triclinic: just look at Cij above')
C11 = 247.86236444973335 GPa
C12 = 147.8284132524967 GPa
C44 = 124.83811766566167 GPa