symmetrizing-the-elastic-tensor-with-crystal-symmetry-operations.md

symmetrizing-the-elastic-tensor-with-crystal-symmetry-operations.md

Introduction

Density functional theory (DFT) has proven to be an effective way to obtain the physical properties of material systems. Using the linear stress-strain relationship, one is able to obtain the elastic constants from a set of DFT calculations. One particular challenge is that the elastic tensor thus obtained is usually not numerically symmetric with respect to the crystal symmetry. Neumann's principle, or principle of symmetry, states that ¹:

If a crystal is invariant with respect to certain symmetry operations, any of its physical properties must also be invariant with respect to the same symmetry operations, or otherwise stated, the symmetry operations of any physical property of a crystal must include the symmetry operations of the point group of the crystal.

Therefore, it is worthwhile working out a routine to symmetrize the elastic tensor using the symmetry operations of the crystal.

Theoretical formulation

A $3\times 3 \times 3\times 3$ elastic tensor ${c_{ijkl}}$ has the intrinsic symmetry $c_{ijkl} = c_{jikl}$ and $c_{ijkl} = c_{ijlk}$, which simplifies the 81 elements to 36. With Voigt's notation, one can use a $6\times 6$ matrix ${C_{ij}}$ to represent the 36 elements. One is able to set $C_{ij} = C_{ji}$, so that 36 elements become 21. Further simplification depends on the symmetry operations of the crystal structure ².

Suppose the space group of the structure is $\mathbf G$, and the point group is $\mathbf R$, it can be stated that the rotational parts of the operations in $\mathbf G$ form $\mathbf R$. For each rotational operation $\mathbf{Q}$ that operates on a column vector $\mathbf v$ ³, $$\mathbf{v}' = \mathbf{Q}⋅\mathbf{v}.$$ Therefore, for the 4th order tensor $\mathbf c$, $$\mathbf{c}' = \mathbf{Q}⋅\mathbf{Q}⋅\mathbf{c}⋅\mathbf{Q}^T⋅\mathbf{Q}^T,$$ or in Einstein's summation, $$ \begin{align} c'{ijkl} & = Q{im} Q_{jn} c_{mnop} (Q^T){ok} (Q^T){pl} \ & = Q_{im} Q_{jn} Q_{ko} Q_{lp} c_{mnop}. \end{align} $$

So the idea is let all the rotations in $\mathbf R$ operate on $\mathbf c$ to get $\mathbf c'$, and average the resultant tensors.

Python implementation

Luckily, our effort doesn't have to go from scratch. There are excellent packages in Python to deal with basic math and linear algebra, such as NumPy⁴, and to deal with crystal structures, pymatgen⁵, leveraging the spglib⁶ to handle symmetry detection and obtain space group operations.

Here we assume users have installed Python 3.x and packages NumPy and pymatgen (v.3.4.0). We define a new function that takes in a $6\times 6$ elastic tensor in Voigt's notation, the structure used in the calculation to symmetrize the tensor.

def get_symmetrized_elastic_tensor(C, structure, align=False, tol=1e-4):
    """
    Parameters:
        C: a 6 x 6 elasic tensor.
        structure: the structure used in the calculation, here given to detect symmetry.
        align: if the tensor is provided in the standard orientation, but the structure is not,
                switch to True to align the structure.
        tol: tolerance for pymatgen.analysis.elasticity.elastic.ElasticTensor initializer,
            here default to 1e-4 rather than 1e-5.

    Returns:
        C_prime: the symmetrized elastic tensor.
    """
    import numpy as np
    from pymatgen.analysis.elasticity.elastic import ElasticTensor
    from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
    # get the space group analyzer of the primitive cell.
    sga = SpacegroupAnalyzer(SpacegroupAnalyzer(structure).find_primitive())
    if align:
        sga = SpacegroupAnalyzer(sga.get_primitive_standard_structure())
	# get the rotational parts of symmetry operations of the space group in the cartesian coordinates.
    ops_list = sga.get_symmetry_operations(cartesian=True)
    ops_rotation_list = [op.rotation_matrix for op in ops_list]
	# symmetrize C with $C_{ij} = C_{ji}$.
    C = ElasticTensor(ElasticTensor(C).symmetrized)
    # average the rotated tensors.
    C_prime_sum = 0
    for idx, op in enumerate(ops_rotation_list):
	    # perform the 4th order tensor rotation.
        C_prime = ElasticTensor.from_full_tensor(np.einsum('im,jn,mnop,ok,pl', op, op, C.full_tensor, op.T, op.T), tol=tol)
        C_prime_sum += C_prime
    C_prime = C_prime_sum / (idx + 1)
    # print out some useful information.
    print(sga.get_spacegroup_symbol(), sga.get_point_group(), len(ops_rotation_list))
    
    return C_prime

Here is a test, getting the structures and elastic tensors from the Materials Project⁷:

import pymatgen as mg
import pandas as pd

def get_MP_data(chemsys_formula_id, data_type='vasp', prop=''):
    m = mg.MPRester()
    df = pd.io.json.json_normalize(m.get_data(chemsys_formula_id, data_type, prop)).set_index('material_id')\
            .drop(['spacegroup.source', 'spacegroup.point_group', 'spacegroup.hall'], axis=1)
    for col in df.columns:
        if 'unit_cell_formula.' in col:
            df.drop(col, axis=1, inplace=True)
    return df

# test
for i in ['Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn']:
    print(i)
    for idx, data in get_MP_data(i + '-N').iterrows():
        if ('elasticity.elastic_tensor' in data) and (not np.isnan(data['elasticity.elastic_tensor']).all()):
            st = mg.Structure.from_str(data['cif'], fmt='cif')
            C = pd.DataFrame(data['elasticity.elastic_tensor'], index=range(1, 7), columns=range(1, 7))
            C_prime = pd.DataFrame(get_symmetrized_elastic_tensor(C, st, align=True), index=range(1, 7), columns=range(1, 7))
            print(C_prime)
            print(C_prime - C)
            print('=====')

Footnotes

1. http://reference.iucr.org/dictionary/Neumann's_principle ↩

2. J.F. Nye, Physical Properties of Crystals: Their Representation by Tensors and Matrices (Clarendon Press, 1985). ↩

3. http://www.continuummechanics.org/cm/coordxforms.html ↩

4. http://www.numpy.org/ ↩

5. http://pymatgen.org/ ↩

6. http://atztogo.github.io/spglib/ ↩

7. http://materialsproject.org/ ↩

@MilMou

From your error traceback, it looks like you have not imported the pymatgen library and aliased it to mg.
It can be done like this:
import pymatgen as mg

And that should solve your error.