5
5
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
6
# SPDX - License - Identifier: GPL - 3.0 +
7
7
from mantid .api import AlgorithmFactory , FileProperty , FileAction , WorkspaceProperty , PythonAlgorithm
8
- from mantid .kernel import Direction
8
+ from mantid .kernel import Direction , V3D
9
9
from mantid .geometry import SymmetryOperationFactory , SpaceGroupFactory
10
10
from os import path , makedirs
11
+ import numpy as np
11
12
12
13
13
14
class SaveINS (PythonAlgorithm ):
14
15
LATT_TYPE_MAP = {type : itype + 1 for itype , type in enumerate (["P" , "I" , "R" , "F" , "A" , "B" , "C" ])}
16
+ IDENTIY_OP = SymmetryOperationFactory .createSymOp ("x,y,z" )
15
17
INVERSION_OP = SymmetryOperationFactory .createSymOp ("-x,-y,-z" )
18
+ ROTATION_OPS = {1 : [IDENTIY_OP , INVERSION_OP ], - 1 : [IDENTIY_OP ]}
19
+ A_CENTERING_OP = SymmetryOperationFactory .createSymOp ("x,y+1/2,z+1/2" )
20
+ B_CENTERING_OP = (SymmetryOperationFactory .createSymOp ("x+1/2,y,z+1/2" ),)
21
+ C_CENTERING_OP = SymmetryOperationFactory .createSymOp ("x+1/2,y+1/2,z" )
22
+ CENTERING_OPS = {
23
+ 1 : [IDENTIY_OP ],
24
+ 2 : [
25
+ IDENTIY_OP ,
26
+ SymmetryOperationFactory .createSymOp ("x+1/2,y+1/2,z+1/2" ),
27
+ ],
28
+ 3 : [
29
+ IDENTIY_OP ,
30
+ SymmetryOperationFactory .createSymOp ("x+1/3,y+2/3,z+2/3" ),
31
+ SymmetryOperationFactory .createSymOp ("x+2/3,y+1/3,z+1/3" ),
32
+ ],
33
+ 4 : [
34
+ IDENTIY_OP ,
35
+ A_CENTERING_OP ,
36
+ B_CENTERING_OP ,
37
+ C_CENTERING_OP ,
38
+ ],
39
+ 5 : [
40
+ IDENTIY_OP ,
41
+ A_CENTERING_OP ,
42
+ ],
43
+ 6 : [IDENTIY_OP , B_CENTERING_OP ],
44
+ 7 : [
45
+ IDENTIY_OP ,
46
+ C_CENTERING_OP ,
47
+ ],
48
+ }
16
49
DUMMY_WAVELENGTH = 1.0
17
50
18
51
def category (self ):
@@ -113,7 +146,7 @@ def PyExec(self):
113
146
f_handle .write (f"LATT { latt_type } \n " )
114
147
115
148
# print sym operations
116
- for sym_str in spgr . getSymmetryOperationStrings ( ):
149
+ for sym_str in self . _get_shelx_symmetry_operators ( spgr , latt_type ):
117
150
f_handle .write (f"SYMM { sym_str } \n " )
118
151
119
152
# print atom info
@@ -140,5 +173,84 @@ def PyExec(self):
140
173
f_handle .write ("HKLF 2\n " ) # tells SHELX the columns saved in the reflection file
141
174
f_handle .write ("END" )
142
175
176
+ def _symmetry_operation_key (self , rot1_mat , trans1_vec , rot2_mat = np .eye (3 ), trans2_vec = np .zeros (3 )):
177
+ """
178
+ Generate a key for symmetry operation comparison.
179
+ Combines rotation and translation into a unique tuple representation.
180
+ Ex: "x,y,z+1/2" is equivalent to "x,y,z+0.5"
181
+ """
182
+ rot_mat = rot1_mat @ rot2_mat
183
+ trans_vec = rot1_mat @ trans2_vec + trans1_vec
184
+ trans_vec = np .mod (trans_vec , 1 ) # Ensure trans_vec is within [0, 1)
185
+ return tuple (np .round (rot_mat , 0 ).astype (int ).flatten ().tolist () + np .round (trans_vec , 3 ).tolist ())
186
+
187
+ def _symmetry_matrix_vector (self , symop ):
188
+ """
189
+ Extract the rotation matrix (rot_mat) and translation vector (trans_vec) from a symmetry element.
190
+ This symmetry operation transform any point via a matrix/translation pair.
191
+ """
192
+ rot_mat = np .linalg .inv (np .vstack ([symop .transformHKL (V3D (* vec )) for vec in np .eye (3 )]))
193
+ trans_vec = np .array (symop .transformCoordinates (V3D (0 , 0 , 0 )))
194
+ return rot_mat , trans_vec
195
+
196
+ def _generate_equivalent_operators (self , rotation_ops , centering_ops ):
197
+ """
198
+ Generate all equivalent symmetry operators for the given lattice rotation and centering operations.
199
+ """
200
+ equivalent_ops = set ()
201
+ for rot in rotation_ops :
202
+ rot2_mat , _ = self ._symmetry_matrix_vector (rot )
203
+ for cent in centering_ops :
204
+ _ , trans2_vec = self ._symmetry_matrix_vector (cent )
205
+ key = self ._symmetry_operation_key (np .eye (3 ), np .zeros (3 ), rot2_mat , trans2_vec )
206
+ equivalent_ops .add (key )
207
+ return equivalent_ops
208
+
209
+ def _update_symmetry_dict (self , rot1_mat , trans1_vec , sym3 , sym_key , sym_ops_dict , rot_mat_dict , trans_vec_dict ):
210
+ """
211
+ Update the symmetry operations dictionary with priority for closeness to identity/origin.
212
+ This bias improves readability.
213
+ # Ex: lattice type 3; "-x+y,-x,z" is simpler than "-x+y+1/3,-x+2/3,z+2/3"
214
+ """
215
+ if sym3 not in sym_ops_dict or (
216
+ np .linalg .det (rot1_mat ) > np .linalg .det (rot_mat_dict [sym3 ]) # identity preferred
217
+ or np .linalg .norm (trans1_vec ) < np .linalg .norm (trans_vec_dict [sym3 ]) # origin preferred
218
+ ):
219
+ sym_ops_dict [sym3 ] = sym_key
220
+ rot_mat_dict [sym3 ] = rot1_mat
221
+ trans_vec_dict [sym3 ] = trans1_vec
222
+
223
+ def _get_shelx_symmetry_operators (self , spgr , latt_type ):
224
+ """
225
+ Get SHELX symmetry operators for the given space group and lattice type.
226
+ Returns symmetry set.
227
+ """
228
+ latt_numb = abs (latt_type )
229
+ latt_sign = 1 if latt_type > 0 else - 1
230
+
231
+ # Generate equivalent lattice type operators common to lattice type.
232
+ latt_type_ops_set = self ._generate_equivalent_operators (self .ROTATION_OPS [latt_sign ], self .CENTERING_OPS [latt_numb ])
233
+
234
+ sym_ops = spgr .getSymmetryOperations ()
235
+ sym_ops_dict = {}
236
+ rot_mat_dict = {}
237
+ trans_vec_dict = {}
238
+
239
+ for sym_op in sym_ops :
240
+ rot1_mat , trans1_vec = self ._symmetry_matrix_vector (sym_op )
241
+ sym_key = sym_op .getIdentifier ()
242
+ sym1 = self ._symmetry_operation_key (rot1_mat , trans1_vec )
243
+
244
+ if sym1 not in latt_type_ops_set :
245
+ # re-iterate over lattice operators to map equivalently generated
246
+ for rot in self .ROTATION_OPS [latt_sign ]:
247
+ rot2_mat , _ = self ._symmetry_matrix_vector (rot )
248
+ for cent in self .CENTERING_OPS [latt_numb ]:
249
+ _ , trans2_vec = self ._symmetry_matrix_vector (cent )
250
+ sym3 = self ._symmetry_operation_key (rot1_mat , trans1_vec , rot2_mat , trans2_vec ) # sym3 = sym1 X sym2
251
+ self ._update_symmetry_dict (rot1_mat , trans1_vec , sym3 , sym_key , sym_ops_dict , rot_mat_dict , trans_vec_dict )
252
+
253
+ return set (sym_ops_dict .values ())
254
+
143
255
144
256
AlgorithmFactory .subscribe (SaveINS )
0 commit comments