# Symmetrization of elastic constants based on rotational symmetry in Pymatgen

I am trying to write a small python script which will symmetrize an elastic tensor based on rotational symmetry of the space group. I started with this link and created a new python script as below.

I started with cubic elastic constants in Voigt form and expected to get same results as input. But the final result shows C12=C14 which is not correct. What might I be doing wrong?

from pymatgen.core import Structure, Lattice
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
import numpy as np
from pymatgen.analysis.elasticity.elastic import ElasticTensor

structure = Structure.from_spacegroup(
"Fm-3m", Lattice.cubic(1.0), ["element"], [[0, 0, 0]])
sga = SpacegroupAnalyzer(SpacegroupAnalyzer(structure).find_primitive())
ops_list = sga.get_symmetry_operations(cartesian=True)
ops_rotation_list = [op.rotation_matrix for op in ops_list]

C = np.array([[271.64, 141.6, 141.6,0,0,0],
[141.6,271.64,141.6,0,0,0],
[141.6,141.6,271.64,0,0,0],
[0,0,0,105.18,0,0],
[0,0,0,0,105.18,0],
[0,0,0,0,0,105.18]])
C_prime_sum = 0

CC = ElasticTensor(ElasticTensor.from_voigt(C).symmetrized)
for idx, op in enumerate(ops_rotation_list):
# perform the 4th order tensor rotation.
C_prime = ElasticTensor(
np.einsum('im,jn,mnop,ok,pl', op, op, CC, op.T, op.T), tol=1e-5)
C_prime_sum += C_prime
C_prime = C_prime_sum / (idx + 1)
A=C_prime.voigt


Expected output (Same as input matrix C)

C = np.array([[271.64, 141.6, 141.6,0,0,0],
[141.6,271.64,141.6,0,0,0],
[141.6,141.6,271.64,0,0,0],
[0,0,0,105.18,0,0],
[0,0,0,0,105.18,0],
[0,0,0,0,0,105.18]])



current output (Which is not correct)

[[271.64 117.32 117.32   0.     0.     0.  ]
[117.32 271.64 117.32   0.     0.     0.  ]
[117.32 117.32 271.64   0.     0.     0.  ]
[  0.     0.     0.   117.32   0.     0.  ]
[  0.     0.     0.     0.   117.32   0.  ]
[  0.     0.     0.     0.     0.   117.32]]


• +1. I was sorry to see that this Bounty didn't result in any answers. I did try to bump up the question to the top of the main page so that it could get more attention while the bounty was active. I've also posted a link to the question here but the author of that file seems to have moved away from working on science. Do you think you can show us the output that you get from running the script in your question, and say what you the expected output would be? Aug 31, 2022 at 20:23
• @NikeDattani I have updated question Sep 1, 2022 at 16:17
• @NikeDattani I am now ready to answer this question. Is it possible to issue a bounty on this question again ? :P Apr 25 at 11:14
• @VandanRevanur second bounty will start from 100. I think if you answer correctly, we may give our old bounty of 50.@NikeDattani can give more details. Apr 25 at 12:20
• I will answer as soon as there is a bounty setup :) Apr 25 at 12:29

### Solution

On Line 21 of your code, where you initialize the elastic tensor, if the .symmetrized part is removed, then you get the desired output.

### Modified Code:

from pymatgen.core import Structure, Lattice
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
import numpy as np
from pymatgen.analysis.elasticity.elastic import ElasticTensor

structure = Structure.from_spacegroup(
"Fm-3m", Lattice.cubic(1.0), ["element"], [[0, 0, 0]])
sga = SpacegroupAnalyzer(SpacegroupAnalyzer(structure).find_primitive())
ops_list = sga.get_symmetry_operations(cartesian=True)
ops_rotation_list = [op.rotation_matrix for op in ops_list]

C = np.array([[271.64, 141.6 , 141.6 ,   0.  ,   0.  ,   0.  ],
[141.6 , 271.64, 141.6 ,   0.  ,   0.  ,   0.  ],
[141.6 , 141.6 , 271.64,   0.  ,   0.  ,   0.  ],
[  0.  ,   0.  ,   0.  , 105.18,   0.  ,   0.  ],
[  0.  ,   0.  ,   0.  ,   0.  , 105.18,   0.  ],
[  0.  ,   0.  ,   0.  ,   0.  ,   0.  , 105.18]])
C_prime_sum = 0

CC = ElasticTensor(ElasticTensor.from_voigt(C))
for idx, op in enumerate(ops_rotation_list):
# perform the 4th order tensor rotation.
C_prime = ElasticTensor(
np.einsum('im,jn,mnop,ok,pl', op, op, CC, op.T, op.T), tol=1e-5)
C_prime_sum += C_prime
C_prime = C_prime_sum / (idx + 1)
A=C_prime.voigt


### Output

print(A) yields the below values which is equal to C:

array([[271.64, 141.6 , 141.6 ,   0.  ,   0.  ,   0.  ],
[141.6 , 271.64, 141.6 ,   0.  ,   0.  ,   0.  ],
[141.6 , 141.6 , 271.64,   0.  ,   0.  ,   0.  ],
[  0.  ,   0.  ,   0.  , 105.18,   0.  ,   0.  ],
[  0.  ,   0.  ,   0.  ,   0.  , 105.18,   0.  ],
[  0.  ,   0.  ,   0.  ,   0.  ,   0.  , 105.18]])