matscipy.elasticity

Computing elastic moduli of atomistic systems under small deformations.

Functions

Voigt_6_to_full_3x3_strain(strain_vector)

Form a 3x3 strain matrix from a 6 component vector in Voigt notation

Voigt_6_to_full_3x3_stress(stress_vector)

Form a 3x3 stress matrix from a 6 component vector in Voigt notation

Voigt_6x6_to_cubic(C)

Convert the Voigt 6x6 representation into the cubic elastic constants C11, C12 and C44.

Voigt_6x6_to_full_3x3x3x3(C)

Convert from the Voigt representation of the stiffness matrix to the full 3x3x3x3 representation.

coalesce_elastic_constants([C11, C12, C44, ...])

Convert either all of C11, C12, C44, or 6x6 C, or 3x3x3x3 C to the form defined by convention.

cubic_to_Voigt_6x6(C11, C12, C44)

elastic_moduli(C[, l, R, tol])

Calculate elastic moduli from 6x6 elastic constant matrix C_{ij}.

fit_elastic_constants(a[, symmetry, ...])

Compute elastic constants by linear regression of stress vs.

full_3x3_to_Voigt_6_index(i, j)

full_3x3_to_Voigt_6_strain(strain_matrix)

Form a 6 component strain vector in Voigt notation from a 3x3 matrix

full_3x3_to_Voigt_6_stress(stress_matrix)

Form a 6 component stress vector in Voigt notation from a 3x3 matrix

full_3x3x3x3_to_Voigt_6x6(C[, tol, ...])

Convert from the full 3x3x3x3 representation of the stiffness matrix to the representation in Voigt notation.

generate_strained_configs(at0[, symmetry, ...])

Generate a sequence of strained configurations

invariants(s[, syy, szz, syz, sxz, sxy, ...])

measure_triclinic_elastic_constants(a[, ...])

Brute-force measurement of elastic constants for a triclinic (general) unit cell.

nonaffine_elastic_contribution(atoms[, ...])

Compute the correction of non-affine displacements to the elasticity tensor.

poisson_ratio(C, l, m)

Calculate approximate Poisson ratio

rotate_cubic_elastic_constants(C11, C12, C44, A)

Return rotated elastic moduli for a cubic crystal given the elastic constant in standard C11, C12, C44 notation.

rotate_elastic_constants(C, A[, tol])

Return rotated elastic moduli for a general crystal given the elastic constant in Voigt notation.

youngs_modulus(C, l)

Calculate approximate Youngs modulus E_l from 6x6 elastic constants matrix C_ij

Classes

CubicElasticModuli(C11, C12, C44)

matscipy.elasticity.full_3x3_to_Voigt_6_index(i, j)
matscipy.elasticity.Voigt_6_to_full_3x3_strain(strain_vector)

Form a 3x3 strain matrix from a 6 component vector in Voigt notation

matscipy.elasticity.Voigt_6_to_full_3x3_stress(stress_vector)

Form a 3x3 stress matrix from a 6 component vector in Voigt notation

matscipy.elasticity.full_3x3_to_Voigt_6_strain(strain_matrix)

Form a 6 component strain vector in Voigt notation from a 3x3 matrix

matscipy.elasticity.full_3x3_to_Voigt_6_stress(stress_matrix)

Form a 6 component stress vector in Voigt notation from a 3x3 matrix

matscipy.elasticity.Voigt_6x6_to_full_3x3x3x3(C)

Convert from the Voigt representation of the stiffness matrix to the full 3x3x3x3 representation.

Parameters:

C (array_like) – 6x6 stiffness matrix (Voigt notation).

Returns:

C – 3x3x3x3 stiffness matrix.

Return type:

array_like

matscipy.elasticity.full_3x3x3x3_to_Voigt_6x6(C, tol=0.001, check_symmetry=True)

Convert from the full 3x3x3x3 representation of the stiffness matrix to the representation in Voigt notation. Checks symmetry in that process.

matscipy.elasticity.Voigt_6x6_to_cubic(C)

Convert the Voigt 6x6 representation into the cubic elastic constants C11, C12 and C44.

matscipy.elasticity.cubic_to_Voigt_6x6(C11, C12, C44)
matscipy.elasticity.coalesce_elastic_constants(C11=None, C12=None, C44=None, C=None, convention='Cij')

Convert either all of C11, C12, C44, or 6x6 C, or 3x3x3x3 C to the form defined by convention.

convention=”cubic”: returns C11, C12, C44 convention=”Cij”: returns 6x6 C matrix convention=”Cijkl”: returns 3x3x3x3 C tensor

matscipy.elasticity.invariants(s, syy=None, szz=None, syz=None, sxz=None, sxy=None, full_3x3_to_Voigt_6=<function full_3x3_to_Voigt_6_stress>)
matscipy.elasticity.rotate_cubic_elastic_constants(C11, C12, C44, A, tol=1e-06)

Return rotated elastic moduli for a cubic crystal given the elastic constant in standard C11, C12, C44 notation.

Parameters:

  • C11 (float) – Cubic elastic moduli.

  • C12 (float) – Cubic elastic moduli.

  • C44 (float) – Cubic elastic moduli.

  • A (array_like) – 3x3 rotation matrix.

Returns:

C – 6x6 matrix of rotated elastic constants (Voigt notation).

Return type:

array

matscipy.elasticity.rotate_elastic_constants(C, A, tol=1e-06)

Return rotated elastic moduli for a general crystal given the elastic constant in Voigt notation.

Parameters:

  • C (array_like) – 6x6 matrix of elastic constants (Voigt notation).

  • A (array_like) – 3x3 rotation matrix.

Returns:

C – 6x6 matrix of rotated elastic constants (Voigt notation).

Return type:

array

class matscipy.elasticity.CubicElasticModuli(C11, C12, C44)

Bases: object

Methods

compliance()

Return the compliance coefficients

rotate(A)

Compute the rotated stiffness matrix

stiffness()

Return the elastic constants
tol = 1e-06
__init__(C11, C12, C44)

Initialize a cubic system with elastic constants C11, C12, C44

rotate(A)

Compute the rotated stiffness matrix

stiffness()

Return the elastic constants

compliance()

Return the compliance coefficients

matscipy.elasticity.measure_triclinic_elastic_constants(a, delta=0.001, optimizer=None, logfile=None, **kwargs)

Brute-force measurement of elastic constants for a triclinic (general) unit cell.

Parameters:

  • a (ase.Atoms) – Atomic configuration.

  • optimizer (ase.optimizer.*) – Optimizer to use for atomic position. Does not optimize atomic position if set to None.

  • delta (float) – Strain increment for analytical derivatives of stresses.

Returns:

C – 6x6 matrix of the elastic constants in Voigt notation.

Return type:

array_like

matscipy.elasticity.generate_strained_configs(at0, symmetry='triclinic', N_steps=5, delta=0.01)

Generate a sequence of strained configurations

Parameters:

  • a (ase.Atoms) – Bulk crystal configuration - both unit cell and atomic positions should be relaxed before calling this routine.

  • symmetry (string) – Symmetry to use to determine which strain patterns are necessary. Default is ‘triclininc’, i.e. no symmetry.

  • N_steps (int) – Number of atomic configurations to generate for each strain pattern. Default is 5. Absolute strain values range from -delta*N_steps/2 to +delta*N_steps/2.

  • delta (float) – Strain increment for analytical derivatives of stresses. Default 1e-2

Returns:

  • Generator which yields a sequence of ase.Atoms instances corresponding

  • to the minima set of strained conifugurations required to evaluate the

  • full 6x6 C_ij matrix under the assumed symmetry.

matscipy.elasticity.fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=0.01, optimizer=None, verbose=True, graphics=False, logfile=None, stress_err=None, grad_scale=0.6241509125883258, intercept_scale=0.03120754562941629, **kwargs)

Compute elastic constants by linear regression of stress vs. strain

Parameters:

  • a (ase.Atoms or list of ase.Atoms) – Either a single atomic configuration or a list of strained configurations. If a single configuration is given, it is passed generate_strained_configs() along with symmetry, N_steps, and delta to generate the set of strained configurations.

  • symmetry (string) – Symmetry to use to determine which strain patterns are necessary. Default is ‘triclininc’, i.e. no symmetry.

  • N_steps (int) – Number of atomic configurations to generate for each strain pattern. Default is 5. Absolute strain values range from -delta*N_steps/2 to +delta*N_steps/2.

  • delta (float) – Strain increment for analytical derivatives of stresses. Default is 1e-2.

  • optimizer (ase.optimizer.*) – Optimizer to use for atomic positions (internal degrees of freedom) for each applied strain. Initial config a is not optimised, and should already have relaxed cell and atomic positions. Does not optimize atomic positions if set to None.

  • verbose (bool) – If True, print additional infomation about the quality of fit and summarise results of C_ij and estimated errors. Default True.

  • graphics (bool) – If True, use matplotlib.pyplot to plot the stress vs. strain curve for each C_ij component fitted. Default True.

  • logfile (bool) – Log file to write optimizer output to. Default None (i.e. suppress output).

  • stress_err (float) – Assumed error on each stress measurement, in eV/A**2

  • grad_scale (float) – Prior assumption of scale of elastic constants in eV/A**2 (default 100GPa)

  • intercept_scale (float) – Prior assumption of scale of intercept in eV/A**2 (default 5GPa)

  • **kwargs (dict) – Additional arguments to pass to optimizer.run() method e.g. fmax.

Returns:

  • C (array_like) – 6x6 matrix of the elastic constants in Voigt notation.

  • C_err (array_like) – If scipy.stats module is available then error estimates for each C_ij component are obtained from the accuracy of the linear regression. Otherwise an array of np.zeros((6,6)) is returned.

Notes

Code originally adapted from elastics.py script, available from http://github.com/djw/elastic-constants

matscipy.elasticity.youngs_modulus(C, l)

Calculate approximate Youngs modulus E_l from 6x6 elastic constants matrix C_ij

This is the modulus for loading in the l direction. For the exact answer, taking into account elastic anisotropuy, rotate the C_ij matrix to the correct frame, compute the compliance matrix:

C = ...  # 6x6 C_ij matrix in crystal frame
A = ...  # rotation matrix
Cr = rotate_elastic_constants(C, A)
S = np.inv(Cr)
E_x = 1/S[0, 0]  # Young's modulus for a pull in x direction
E_y = 1/S[1, 1]  # Young's modulus for a pull in y direction
E_z = 1/S[0, 0]  # Young's modulus for a pull in z direction

Notes

Formula is from W. Brantley, Calculated elastic constants for stress problems associated with semiconductor devices. J. Appl. Phys., 44, 534 (1973).

matscipy.elasticity.poisson_ratio(C, l, m)

Calculate approximate Poisson ratio

u_{lm} from 6x6 elastic constant matrix C_{ij}

This is the response in m direction to pulling in l direction. Result is dimensionless.

Notes

Formula is from W. Brantley, Calculated elastic constants for stress problems associated with semiconductor devices. J. Appl. Phys., 44, 534 (1973).

matscipy.elasticity.elastic_moduli(C, l=array([1, 0, 0]), R=None, tol=1e-06)

Calculate elastic moduli from 6x6 elastic constant matrix C_{ij}.

The elastic moduli calculated are: Young’s muduli, Poisson’s ratios, shear moduli, bulk mudulus and bulk mudulus tensor.

If a direction l is specified, the system is rotated to have it as its x direction (see Notes for details). If R is specified the system is rotated according to it.

Parameters:

  • C (array_like) – 6x6 stiffness matrix (Voigt notation).

  • l (array_like, optional) – 3D direction vector for pull (the default is the x direction of the original system)

  • R (array_like, optional) – 3x3 rotation matrix.

  • tol (float, optional) – tolerance for checking validity of rotation and comparison of vectors.

Returns:

  • E (array_like) – Young’s modulus for a stress in each of the three directions of the rotated system.

  • nu (array_like) – 3x3 matrix with Poisson’s ratios.

  • Gm (array_like) – 3x3 matrix with shear moduli.

  • B (float) – Bulk modulus.

  • K (array_like) – 3x3 matrix with bulk modulus tensor.

Notes

It works by rotating the elastic constant tensor to the desired direction, so it should be valid for arbitrary crystal structures. If only l is specified there is an infinite number of possible rotations. The chosen one is a rotation along the axis orthogonal to the plane defined by the vectors (1, 0, 0) and l.

Bulk modulus tensor as defined in O. Rand and V. Rovenski, “Analytical Methods in Anisotropic Elasticity”, Birkh”auser (2005), pp. 71.

matscipy.elasticity.nonaffine_elastic_contribution(atoms, eigenvalues=None, eigenvectors=None, pc_parameters=None, cg_parameters={'M': None, 'atol': 1e-05, 'callback': None, 'maxiter': None, 'tol': 1e-05, 'x0': None})

Compute the correction of non-affine displacements to the elasticity tensor. The computation of the occuring inverse of the Hessian matrix is bypassed by using a cg solver.

If eigenvalues and and eigenvectors are given the inverse of the Hessian can be easily computed.

Parameters:

  • atoms (ase.Atoms) – Atomic configuration in a local or global minima.

  • eigenvalues (array) – Eigenvalues in ascending order obtained by diagonalization of Hessian matrix. If given, use eigenvalues and eigenvectors to compute non-affine contribution.

  • eigenvectors (array) – Eigenvectors corresponding to eigenvalues.

  • cg_parameters (dict) –

    Dictonary for the conjugate-gradient solver.

    x0: {array, matrix}

    Starting guess for the solution.

    tol/atol: float, optional

    Tolerances for convergence, norm(residual) <= max(tol*norm(b), atol).

    maxiter: int

    Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved.

    M: {sparse matrix, dense matrix, LinearOperator}

    Preconditioner for A.

    callback: function

    User-supplied function to call after each iteration.

  • pc_parameters (dict) –

    Dictonary for the incomplete LU decomposition of the Hessian.

    A: array_like

    Sparse matrix to factorize.

    drop_tol: float

    Drop tolerance for an incomplete LU decomposition.

    fill_factor: float

    Specifies the fill ratio upper bound.

    drop_rule: str

    Comma-separated string of drop rules to use.

    permc_spec: str

    How to permute the columns of the matrix for sparsity.

    diag_pivot_thresh: float

    Threshold used for a diagonal entry to be an acceptable pivot.

    relax: int

    Expert option for customizing the degree of relaxing supernodes.

    panel_size: int

    Expert option for customizing the panel size.

    options: dict

    Dictionary containing additional expert options to SuperLU.