6
$\begingroup$

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]]


$\endgroup$
7
  • $\begingroup$ +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? $\endgroup$ Commented Aug 31, 2022 at 20:23
  • $\begingroup$ @NikeDattani I have updated question $\endgroup$ Commented Sep 1, 2022 at 16:17
  • $\begingroup$ @NikeDattani I am now ready to answer this question. Is it possible to issue a bounty on this question again ? :P $\endgroup$ Commented Apr 25, 2023 at 11:14
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Apr 25, 2023 at 12:20
  • $\begingroup$ I will answer as soon as there is a bounty setup :) $\endgroup$ Commented Apr 25, 2023 at 12:29

1 Answer 1

2
+100
$\begingroup$

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]])
$\endgroup$
1
  • $\begingroup$ I am not getting nuances hidden here. What is the use of symmetrization of elastic constants if it does not bring any changes to final constants? Is it only to check if they are numerically symmetric with respect to crystal symmetry? Apart from 6x6 tensor matrix, is there anything that gets modified and needs to be used further? $\endgroup$
    – AbPhys
    Commented May 15 at 12:21

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .