5
$\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]]
Improve this question
$\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$
    – Nike Dattani
    Aug 31, 2022 at 20:23
  • $\begingroup$ @NikeDattani I have updated question $\endgroup$
    – Pranav kumar
    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$
    – Vandan Revanur
    Apr 25 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$
    – Pranav kumar
    Apr 25 at 12:20
  • $\begingroup$ I will answer as soon as there is a bounty setup :) $\endgroup$
    – Vandan Revanur
    Apr 25 at 12:29

1 Answer 1

Reset to default
1
+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]])
Improve this answer
$\endgroup$

You must log in to answer this question.

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