Context Navigation

source: trunk/GSASIIspc.py @ 3379

Last change on this file since 3379 was 3379, checked in by vondreele, 7 years ago
fix problem of bad space group symbol (no field spaces) found in cif file
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 153.2 KB

Line
1	# -- coding: utf-8 --
2	"""
3	GSASIIspc: Space group module
4	-------------------------------
5
6	Space group interpretation routines. Note that space group information is
7	stored in a :ref:`Space Group (SGData)<SGData_table>` object.
8
9	"""
10	########### SVN repository information ###################
11	# $Date: 2018-05-10 17:48:20 +0000 (Thu, 10 May 2018) $
12	# $Author: vondreele $
13	# $Revision: 3379 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 3379 2018-05-10 17:48:20Z vondreele $
16	########### SVN repository information ###################
17	from __future__ import division, print_function
18	import numpy as np
19	import numpy.linalg as nl
20	import scipy.optimize as so
21	import sys
22	import copy
23	import os.path as ospath
24
25	import GSASIIpath
26	GSASIIpath.SetVersionNumber("$Revision: 3379 $")
27
28	npsind = lambda x: np.sin(x*np.pi/180.)
29	npcosd = lambda x: np.cos(x*np.pi/180.)
30	DEBUG = False
31
32	################################################################################
33	#### Space group codes
34	################################################################################
35
36	def SpcGroup(SGSymbol):
37	"""
38	Determines cell and symmetry information from a short H-M space group name
39
40	:param SGSymbol: space group symbol (string) with spaces between axial fields
41	:returns: (SGError,SGData)
42
43	* SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
44	* SGData - is a dict (see :ref:`Space Group object<SGData_table>`) with entries:
45
46	* 'SpGrp': space group symbol, slightly cleaned up
47	* 'SGFixed': True if space group data can not be changed, e.g. from magnetic cif; otherwise False
48	* 'SGGray': True if 1' in symbol - gray group for mag. incommensurate phases
49	* 'SGLaue': one of '-1', '2/m', 'mmm', '4/m', '4/mmm', '3R',
50	'3mR', '3', '3m1', '31m', '6/m', '6/mmm', 'm3', 'm3m'
51	* 'SGInv': boolean; True if centrosymmetric, False if not
52	* 'SGLatt': one of 'P', 'A', 'B', 'C', 'I', 'F', 'R'
53	* 'SGUniq': one of 'a', 'b', 'c' if monoclinic, '' otherwise
54	* 'SGCen': cell centering vectors [0,0,0] at least
55	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
56	* 'SGSys': one of 'triclinic', 'monoclinic', 'orthorhombic',
57	'tetragonal', 'rhombohedral', 'trigonal', 'hexagonal', 'cubic'
58	* 'SGPolax': one of ' ', 'x', 'y', 'x y', 'z', 'x z', 'y z',
59	'xyz', '111' for arbitrary axes
60	* 'SGPtGrp': one of 32 point group symbols (with some permutations), which
61	is filled by SGPtGroup, is external (KE) part of supersymmetry point group
62	* 'SSGKl': default internal (Kl) part of supersymmetry point group; modified
63	in supersymmetry stuff depending on chosen modulation vector for Mono & Ortho
64	* 'BNSlattsym': BNS lattice symbol & cenering op - used for magnetic structures
65
66	"""
67	LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
68	LattSym = ('P','A','B','C','I','F','R')
69	UniqSym = ('','','a','b','c','',)
70	SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
71	SGData = {}
72	if len(SGSymbol.split()) < 2:
73	return SGErrors(0),SGData
74	if ':R' in SGSymbol:
75	SGSymbol = SGSymbol.replace(':',' ') #get rid of ':' in R space group symbols from some cif files
76	SGData['SGGray'] = False
77	if "1'" in SGSymbol: #set for incommensurate magnetic
78	SGData['SGGray'] = True
79	SGSymbol = SGSymbol.replace("1'",'')
80	SGSymbol = SGSymbol.split(':')[0] #remove :1/2 setting symbol from some cif files
81	if '-2' in SGSymbol: #replace bad but legal symbols with correct equivalents
82	SGSymbol = SGSymbol.replace('-2','m')
83	if SGSymbol.split()[1] =='3/m':
84	SGSymbol = SGSymbol.replace('3/m','-6')
85	import pyspg
86	SGInfo = pyspg.sgforpy(SGSymbol)
87	SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
88	SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
89	SGData['SGInv'] = bool(SGInfo[1])
90	SGData['SGLatt'] = LattSym[SGInfo[2]-1]
91	SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
92	SGData['SGFixed'] = False
93	SGData['SGOps'] = []
94	SGData['SGGen'] = []
95	for i in range(SGInfo[5]):
96	Mat = np.array(SGInfo[6][i])
97	Trns = np.array(SGInfo[7][i])
98	SGData['SGOps'].append([Mat,Trns])
99	if 'array' in str(type(SGInfo[8])): #patch for old fortran bin?
100	SGData['SGGen'].append(int(SGInfo[8][i]))
101	SGData['BNSlattsym'] = [LattSym[SGInfo[2]-1],[0,0,0]]
102	lattSpin = []
103	if SGData['SGLatt'] == 'P':
104	SGData['SGCen'] = np.array(([0,0,0],))
105	elif SGData['SGLatt'] == 'A':
106	SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
107	lattSpin += [1,]
108	elif SGData['SGLatt'] == 'B':
109	SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
110	lattSpin += [1,]
111	elif SGData['SGLatt'] == 'C':
112	SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
113	lattSpin += [1,]
114	elif SGData['SGLatt'] == 'I':
115	SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
116	lattSpin += [1,]
117	elif SGData['SGLatt'] == 'F':
118	SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
119	lattSpin += [1,1,1,1]
120	elif SGData['SGLatt'] == 'R':
121	SGData['SGCen'] = np.array(([0,0,0],[2./3,1./3,1./3],[1./3,2./3,2./3]))
122
123	if SGData['SGInv']:
124	if SGData['SGLaue'] in ['-1','2/m','mmm']:
125	Ibar = 7
126	elif SGData['SGLaue'] in ['4/m','4/mmm']:
127	Ibar = 1
128	elif SGData['SGLaue'] in ['3R','3mR','3','3m1','31m','6/m','6/mmm']:
129	Ibar = 15 #8+4+2+1
130	else:
131	Ibar = 4
132	Ibarx = Ibar&14
133	else:
134	Ibarx = 8
135	if SGData['SGLaue'] in ['-1','2/m','mmm','m3','m3m']:
136	Ibarx = 0
137	moregen = []
138	for i,gen in enumerate(SGData['SGGen']):
139	if SGData['SGLaue'] in ['m3','m3m']:
140	if gen in [1,2,4]:
141	SGData['SGGen'][i] = 4
142	elif gen < 7:
143	SGData['SGGen'][i] = 0
144	elif SGData['SGLaue'] in ['4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm']:
145	if gen == 2:
146	SGData['SGGen'][i] = 4
147	elif gen in [3,5]:
148	SGData['SGGen'][i] = 3
149	elif gen == 6:
150	if SGData['SGLaue'] in ['4/m','4/mmm']:
151	SGData['SGGen'][i] = 128
152	else:
153	SGData['SGGen'][i] = 16
154	elif not SGData['SGInv'] and gen == 12:
155	SGData['SGGen'][i] = 8
156	elif (not SGData['SGInv']) and (SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']) and (gen == 1):
157	SGData['SGGen'][i] = 24
158	gen = SGData['SGGen'][i]
159	if gen == 99:
160	gen = 8
161	if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
162	gen = 3
163	elif SGData['SGLaue'] == 'm3m':
164	gen = 12
165	SGData['SGGen'][i] = gen
166	elif gen == 98:
167	gen = 8
168	if SGData['SGLaue'] in ['3m1','31m','6/m','6/mmm']:
169	gen = 4
170	SGData['SGGen'][i] = gen
171	elif not SGData['SGInv'] and gen in [23,] and SGData['SGLaue'] in ['m3','m3m']:
172	SGData['SGGen'][i] = 24
173	elif gen >= 16 and gen != 128:
174	if not SGData['SGInv']:
175	gen = 31
176	else:
177	gen ^= Ibarx
178	SGData['SGGen'][i] = gen
179	if SGData['SGInv']:
180	if gen < 128:
181	moregen.append(SGData['SGGen'][i]^Ibar)
182	else:
183	moregen.append(1)
184	SGData['SGGen'] += moregen
185	if SGData['SGLaue'] in '-1':
186	SGData['SGSys'] = SysSym[0]
187	elif SGData['SGLaue'] in '2/m':
188	SGData['SGSys'] = SysSym[1]
189	elif SGData['SGLaue'] in 'mmm':
190	SGData['SGSys'] = SysSym[2]
191	elif SGData['SGLaue'] in ['4/m','4/mmm']:
192	SGData['SGSys'] = SysSym[3]
193	elif SGData['SGLaue'] in ['3R','3mR']:
194	SGData['SGSys'] = SysSym[4]
195	elif SGData['SGLaue'] in ['3','3m1','31m']:
196	SGData['SGSys'] = SysSym[5]
197	elif SGData['SGLaue'] in ['6/m','6/mmm']:
198	SGData['SGSys'] = SysSym[6]
199	elif SGData['SGLaue'] in ['m3','m3m']:
200	SGData['SGSys'] = SysSym[7]
201	SGData['SGPolax'] = SGpolar(SGData)
202	SGData['SGPtGrp'],SGData['SSGKl'] = SGPtGroup(SGData)
203
204	if SGData['SGLatt'] == 'R':
205	if SGData['SGPtGrp'] in ['3',]:
206	SGData['SGSpin'] = 4*[1,]
207	elif SGData['SGPtGrp'] in ['-3','32','3m']:
208	SGData['SGSpin'] = 5*[1,]
209	elif SGData['SGPtGrp'] in ['-3m',]:
210	SGData['SGSpin'] = 6*[1,]
211
212	else:
213	if SGData['SGPtGrp'] in ['1','3','23',]:
214	SGData['SGSpin'] = lattSpin+[1,]
215	elif SGData['SGPtGrp'] in ['-1','2','m','4','-4','-3','312','321','3m1','31m','6','-6','432','-43m']:
216	SGData['SGSpin'] = lattSpin+[1,1,]
217	elif SGData['SGPtGrp'] in ['2/m','4/m','422','4mm','-42m','-4m2','-3m1','-31m',
218	'6/m','622','6mm','-6m2','-62m','m3','m3m']:
219	SGData['SGSpin'] = lattSpin+[1,1,1,]
220	else: #'222'-'mmm','4/mmm','6/mmm'
221	SGData['SGSpin'] = lattSpin+[1,1,1,1,]
222	return SGInfo[-1],SGData
223
224	def SGErrors(IErr):
225	'''
226	Interprets the error message code from SpcGroup. Used in SpaceGroup.
227
228	:param IErr: see SGError in :func:`SpcGroup`
229	:returns:
230	ErrString - a string with the error message or "Unknown error"
231	'''
232
233	ErrString = [' ',
234	'Less than 2 operator fields were found',
235	'Illegal Lattice type, not P, A, B, C, I, F or R',
236	'Rhombohedral lattice requires a 3-axis',
237	'Minus sign does not preceed 1, 2, 3, 4 or 6',
238	'Either a 5-axis anywhere or a 3-axis in field not allowed',
239	' ',
240	'I for COMPUTED GO TO out of range.',
241	'An a-glide mirror normal to A not allowed',
242	'A b-glide mirror normal to B not allowed',
243	'A c-glide mirror normal to C not allowed',
244	'D-glide in a primitive lattice not allowed',
245	'A 4-axis not allowed in the 2nd operator field',
246	'A 6-axis not allowed in the 2nd operator field',
247	'More than 24 matrices needed to define group',
248	' ',
249	'Improper construction of a rotation operator',
250	'Mirror following a / not allowed',
251	'A translation conflict between operators',
252	'The 2bar operator is not allowed',
253	'3 fields are legal only in R & m3 cubic groups',
254	'Syntax error. Expected I -4 3 d at this point',
255	' ',
256	'A or B centered tetragonal not allowed',
257	' ','unknown error in sgroup',' ',' ',' ',
258	'Illegal character in the space group symbol',
259	]
260	try:
261	return ErrString[IErr]
262	except:
263	return "Unknown error"
264
265	def SGpolar(SGData):
266	'''
267	Determine identity of polar axes if any
268	'''
269	POL = ('','x','y','x y','z','x z','y z','xyz','111')
270	NP = [1,2,4]
271	NPZ = [0,1]
272	for M,T in SGData['SGOps']:
273	for i in range(3):
274	if M[i][i] <= 0.: NP[i] = 0
275	if M[0][2] > 0: NPZ[0] = 8
276	if M[1][2] > 0: NPZ[1] = 0
277	NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]NPZ[1])(1-int(SGData['SGInv']))
278	return POL[NPol]
279
280	def SGPtGroup(SGData):
281	'''
282	Determine point group of the space group - done after space group symbol has
283	been evaluated by SpcGroup. Only short symbols are allowed
284
285	:param SGData: from :func SpcGroup
286	:returns: SSGPtGrp & SSGKl (only defaults for Mono & Ortho)
287	'''
288	Flds = SGData['SpGrp'].split()
289	if len(Flds) < 2:
290	return '',[]
291	if SGData['SGLaue'] == '-1': #triclinic
292	if '-' in Flds[1]:
293	return '-1',[-1,]
294	else:
295	return '1',[1,]
296	elif SGData['SGLaue'] == '2/m': #monoclinic - default for 2D modulation vector
297	if '/' in SGData['SpGrp']:
298	return '2/m',[-1,1]
299	elif '2' in SGData['SpGrp']:
300	return '2',[-1,]
301	else:
302	return 'm',[1,]
303	elif SGData['SGLaue'] == 'mmm': #orthorhombic
304	if SGData['SpGrp'].count('2') == 3:
305	return '222',[-1,-1,-1]
306	elif SGData['SpGrp'].count('2') == 1:
307	if SGData['SGPolax'] == 'x':
308	return '2mm',[-1,1,1]
309	elif SGData['SGPolax'] == 'y':
310	return 'm2m',[1,-1,1]
311	elif SGData['SGPolax'] == 'z':
312	return 'mm2',[1,1,-1]
313	else:
314	return 'mmm',[1,1,1]
315	elif SGData['SGLaue'] == '4/m': #tetragonal
316	if '/' in SGData['SpGrp']:
317	return '4/m',[1,-1]
318	elif '-' in Flds[1]:
319	return '-4',[-1,]
320	else:
321	return '4',[1,]
322	elif SGData['SGLaue'] == '4/mmm':
323	if '/' in SGData['SpGrp']:
324	return '4/mmm',[1,-1,1,1]
325	elif '-' in Flds[1]:
326	if '2' in Flds[2]:
327	return '-42m',[-1,-1,1]
328	else:
329	return '-4m2',[-1,1,-1]
330	elif '2' in Flds[2:]:
331	return '422',[1,-1,-1]
332	else:
333	return '4mm',[1,1,1]
334	elif SGData['SGLaue'] in ['3','3R']: #trigonal/rhombohedral
335	if '-' in Flds[1]:
336	return '-3',[-1,]
337	else:
338	return '3',[1,]
339	elif SGData['SGLaue'] == '3mR' or 'R' in Flds[0]:
340	if '2' in Flds[2]:
341	return '32',[1,-1]
342	elif '-' in Flds[1]:
343	return '-3m',[-1,1]
344	else:
345	return '3m',[1,1]
346	elif SGData['SGLaue'] == '3m1':
347	if '2' in Flds[2]:
348	return '321',[1,-1,1]
349	elif '-' in Flds[1]:
350	return '-3m1',[-1,1,1]
351	else:
352	return '3m1',[1,1,1]
353	elif SGData['SGLaue'] == '31m':
354	if '2' in Flds[3]:
355	return '312',[1,1,-1]
356	elif '-' in Flds[1]:
357	return '-31m',[-1,1,1]
358	else:
359	return '31m',[1,1,1]
360	elif SGData['SGLaue'] == '6/m': #hexagonal
361	if '/' in SGData['SpGrp']:
362	return '6/m',[1,-1]
363	elif '-' in SGData['SpGrp']:
364	return '-6',[-1,]
365	else:
366	return '6',[1,]
367	elif SGData['SGLaue'] == '6/mmm':
368	if '/' in SGData['SpGrp']:
369	return '6/mmm',[1,-1,1,1]
370	elif '-' in Flds[1]:
371	if '2' in Flds[2]:
372	return '-62m',[-1,-1,1]
373	else:
374	return '-6m2',[-1,1,-1]
375	elif '2' in Flds[2:]:
376	return '622',[1,-1,-1]
377	else:
378	return '6mm',[1,1,1]
379	elif SGData['SGLaue'] == 'm3': #cubic - no (3+1) supersymmetry
380	if '2' in Flds[1]:
381	return '23',[]
382	else:
383	return 'm3',[]
384	elif SGData['SGLaue'] == 'm3m':
385	if '4' in Flds[1]:
386	if '-' in Flds[1]:
387	return '-43m',[]
388	else:
389	return '432',[]
390	else:
391	return 'm3m',[]
392
393	def SGPrint(SGData,AddInv=False):
394	'''
395	Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
396
397	:param SGData: from :func:`SpcGroup`
398	:returns:
399	SGText - list of strings with the space group details
400	SGTable - list of strings for each of the operations
401	'''
402	if SGData.get('SGFixed',False): #inverses included in ops for cif fixed
403	Mult = len(SGData['SGCen'])*len(SGData['SGOps'])
404	else:
405	Mult = len(SGData['SGCen'])len(SGData['SGOps'])(int(SGData['SGInv'])+1)
406	SGText = []
407	SGText.append(' Space Group: '+SGData['SpGrp'])
408	if SGData.get('SGGray',False): SGText[-1] += " 1'"
409	CentStr = 'centrosymmetric'
410	if not SGData['SGInv']:
411	CentStr = 'non'+CentStr
412	if SGData['SGLatt'] in 'ABCIFR':
413	SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
414	else:
415	SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
416	SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
417	if 'SGPtGrp' in SGData: #patch
418	SGText.append(' The lattice point group is '+SGData['SGPtGrp'])
419	SGText.append(' Multiplicity of a general site is '+str(Mult))
420	if SGData['SGUniq'] in ['a','b','c']:
421	SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
422	if SGData['SGInv']:
423	SGText.append(' The inversion center is located at 0,0,0')
424	if SGData['SGPolax']:
425	SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
426	SGText.append(' ')
427	if len(SGData['SGCen']) == 1:
428	SGText.append(' The equivalent positions are:\n')
429	else:
430	SGText.append(' The equivalent positions are:\n')
431	SGText.append(' ('+Latt2text(SGData['SGCen'])+')+\n')
432	SGTable = []
433	for i,Opr in enumerate(SGData['SGOps']):
434	SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
435	if AddInv and SGData['SGInv']:
436	for i,Opr in enumerate(SGData['SGOps']):
437	IOpr = [-Opr[0],-Opr[1]]
438	SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))
439	# if SGData.get('SGGray',False) and not SGData.get('SGFixed',False):
440	# SGTable.append(" for 1'")
441	# for i,Opr in enumerate(SGData['SGOps']):
442	# SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
443	# if AddInv and SGData['SGInv']:
444	# for i,Opr in enumerate(SGData['SGOps']):
445	# IOpr = [-Opr[0],-Opr[1]]
446	# SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))
447	return SGText,SGTable
448
449	def AllOps(SGData):
450	'''
451	Returns a list of all operators for a space group, including those for
452	centering and a center of symmetry
453
454	:param SGData: from :func:`SpcGroup`
455	:returns: (SGTextList,offsetList,symOpList,G2oprList) where
456
457	* SGTextList: a list of strings with formatted and normalized
458	symmetry operators.
459	* offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
460	symmetry operation to the operator in SGTextList and symOpList.
461	these dx (etc.) values are added to the GSAS-II generated
462	positions to provide the positions that are generated
463	by the normalized symmetry operators.
464	* symOpList: a list of tuples with the normalized symmetry
465	operations as (M,T) values
466	(see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
467	* G2oprList: The GSAS-II operations for each symmetry operation as
468	a tuple with (center,mult,opnum,opcode), where center is (0,0,0), (0.5,0,0),
469	(0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
470	where opnum is the number for the symmetry operation, in ``SGOps``
471	(starting with 0) and opcode is mult(100icen+j+1).
472	'''
473	SGTextList = []
474	offsetList = []
475	symOpList = []
476	G2oprList = []
477	G2opcodes = []
478	onebar = (1,)
479	if SGData['SGInv']:
480	onebar += (-1,)
481	for icen,cen in enumerate(SGData['SGCen']):
482	for mult in onebar:
483	for j,(M,T) in enumerate(SGData['SGOps']):
484	offset = [0,0,0]
485	Tprime = (mult*T)+cen
486	for i in range(3):
487	while Tprime[i] < 0:
488	Tprime[i] += 1
489	offset[i] += 1
490	while Tprime[i] >= 1:
491	Tprime[i] += -1
492	offset[i] += -1
493	Opr = [mult*M,Tprime]
494	OPtxt = MT2text(Opr)
495	SGTextList.append(OPtxt.replace(' ',''))
496	offsetList.append(tuple(offset))
497	symOpList.append((mult*M,Tprime))
498	G2oprList.append((cen,mult,j))
499	G2opcodes.append(mult(100icen+j+1))
500	return SGTextList,offsetList,symOpList,G2oprList,G2opcodes
501
502	def MT2text(Opr,reverse=False):
503	"From space group matrix/translation operator returns text version"
504	XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
505	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
506	Fld = ''
507	M,T = Opr
508	for j in range(3):
509	IJ = int(round(2M[j][0]+3M[j][1]+4*M[j][2]+4))%12
510	IK = int(round(T[j]*12))%12
511	if IK:
512	if IJ < 3:
513	if reverse:
514	Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
515	else:
516	Fld += (TRA[IK]+XYZ[IJ]).rjust(5)
517	else:
518	if reverse:
519	Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
520	else:
521	Fld += (TRA[IK]+'+'+XYZ[IJ]).rjust(5)
522	else:
523	Fld += XYZ[IJ].rjust(5)
524	if j != 2: Fld += ', '
525	return Fld
526
527	def Latt2text(Cen):
528	"From lattice centering vectors returns ';' delimited cell centering vectors"
529	lattTxt = ''
530	fracList = ['1/2','1/3','2/3','1/4','3/4','1/5','2/5','3/5','4/5','1/6','5/6',
531	'1/7','2/7','3/7','4/7','5/7','6/7','1/8','3/8','5/8','7/8','1/9','2/9','4/9','5/9','7/9','8/9']
532	mulList = [2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,9,9,9,9,9,9]
533	prodList = [1.,1.,2.,1.,3.,1.,2.,3.,4.,1.,5.,1.,2.,3.,4.,5.,6.,1.,3.,5.,7.,1.,2.,4.,5.,7.,8.]
534	nCen = len(Cen)
535	for i,cen in enumerate(Cen):
536	txt = ''
537	for icen in cen:
538	if not icen:
539	txt += '0,'
540	continue
541	for mul,prod,frac in zip(mulList,prodList,fracList):
542	if abs(icen*mul-prod) < 1.e-5:
543	txt += frac+','
544	break
545	lattTxt += txt[:-1]+'; '
546	if i and not i%8 and i < nCen-1: #not for the last cen!
547	lattTxt += '\n '
548	return lattTxt[:-2]
549
550	def SpaceGroup(SGSymbol):
551	'''
552	Print the output of SpcGroup in a nicely formatted way.
553
554	:param SGSymbol: space group symbol (string) with spaces between axial fields
555	:returns: nothing
556	'''
557	E,A = SpcGroup(SGSymbol)
558	if E > 0:
559	print (SGErrors(E))
560	return
561	for l in SGPrint(A):
562	print (l)
563	################################################################################
564	#### Magnetic space group stuff
565	################################################################################
566
567	def GetGenSym(SGData):
568	'''
569	Get the space group generator symbols
570	:param SGData: from :func:`SpcGroup`
571	LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
572	LattSym = ('P','A','B','C','I','F','R')
573
574	'''
575	OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
576	if SGData['SGInv']:
577	OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
578	Nsyms = len(SGData['SGOps'])
579	if SGData['SGInv'] and not SGData['SGFixed']: Nsyms *= 2
580	UsymOp = ['1',]
581	OprFlg = [0,]
582	if Nsyms == 2: #Centric triclinic or acentric monoclinic
583	UsymOp.append(OprNames[1])
584	OprFlg.append(SGData['SGGen'][1])
585	elif Nsyms == 4: #Point symmetry 2/m, 222, 22m, or 4
586	if '4z' in OprNames[1]: #Point symmetry 4 or -4
587	UsymOp.append(OprNames[1])
588	OprFlg.append(SGData['SGGen'][1])
589	elif not SGData['SGInv']: #Acentric Orthorhombic
590	if 'm' in OprNames[1:4]: #22m, 2m2 or m22
591	if '2' in OprNames[1]: #Acentric orthorhombic, 2mm
592	UsymOp.append(OprNames[2])
593	OprFlg.append(SGData['SGGen'][2])
594	UsymOp.append(OprNames[3])
595	OprFlg.append(SGData['SGGen'][3])
596	elif '2' in OprNames[2]: #Acentric orthorhombic, m2m
597	UsymOp.append(OprNames[1])
598	OprFlg.append(SGData['SGGen'][1])
599	UsymOp.append(OprNames[3])
600	OprFlg.append(SGData['SGGen'][3])
601	else: #Acentric orthorhombic, mm2
602	UsymOp.append(OprNames[1])
603	OprFlg.append(SGData['SGGen'][1])
604	UsymOp.append(OprNames[2])
605	OprFlg.append(SGData['SGGen'][2])
606	else: #Acentric orthorhombic, 222
607	SGData['SGGen'][1:] = [4,2,1]
608	UsymOp.append(OprNames[1])
609	OprFlg.append(SGData['SGGen'][1])
610	UsymOp.append(OprNames[2])
611	OprFlg.append(SGData['SGGen'][2])
612	UsymOp.append(OprNames[3])
613	OprFlg.append(SGData['SGGen'][3])
614	else: #Centric Monoclinic
615	UsymOp.append(OprNames[1])
616	OprFlg.append(SGData['SGGen'][1])
617	UsymOp.append(OprNames[3])
618	OprFlg.append(SGData['SGGen'][3])
619	elif Nsyms == 6: #Point symmetry 32, 3m or 6
620	if '6' in OprNames[1]: #Hexagonal 6/m Laue symmetry
621	UsymOp.append(OprNames[1])
622	OprFlg.append(SGData['SGGen'][1])
623	else: #Trigonal
624	UsymOp.append(OprNames[4])
625	OprFlg.append(SGData['SGGen'][3])
626	if '2110' in OprNames[1]: UsymOp[-1] = ' 2100 '
627	elif Nsyms == 8: #Point symmetry mmm, 4/m, or 422, etc
628	if '4' in OprNames[1]: #Tetragonal
629	if SGData['SGInv']: #4/m
630	UsymOp.append(OprNames[1])
631	OprFlg.append(SGData['SGGen'][1])
632	UsymOp.append(OprNames[6])
633	OprFlg.append(SGData['SGGen'][6])
634	else:
635	if 'x' in OprNames[4]: #4mm type group
636	UsymOp.append(OprNames[4])
637	OprFlg.append(6)
638	UsymOp.append(OprNames[7])
639	OprFlg.append(8)
640	else: #-42m, -4m2, and 422 type groups
641	UsymOp.append(OprNames[5])
642	OprFlg.append(8)
643	UsymOp.append(OprNames[6])
644	OprFlg.append(19)
645	else: #Orthorhombic, mmm
646	UsymOp.append(OprNames[1])
647	OprFlg.append(SGData['SGGen'][1])
648	UsymOp.append(OprNames[2])
649	OprFlg.append(SGData['SGGen'][2])
650	UsymOp.append(OprNames[7])
651	OprFlg.append(SGData['SGGen'][7])
652	elif Nsyms == 12 and '3' in OprNames[1] and SGData['SGSys'] != 'cubic': #Trigonal
653	UsymOp.append(OprNames[3])
654	OprFlg.append(SGData['SGGen'][3])
655	UsymOp.append(OprNames[9])
656	OprFlg.append(SGData['SGGen'][9])
657	elif Nsyms == 12 and '6' in OprNames[1]: #Hexagonal
658	if 'mz' in OprNames[9]: #6/m
659	UsymOp.append(OprNames[1])
660	OprFlg.append(SGData['SGGen'][1])
661	UsymOp.append(OprNames[6])
662	OprFlg.append(SGData['SGGen'][6])
663	else: #6mm, -62m, -6m2 or 622
664	UsymOp.append(OprNames[6])
665	OprFlg.append(18)
666	if 'm' in UsymOp[-1]: OprFlg[-1] = 20
667	UsymOp.append(OprNames[7])
668	OprFlg.append(24)
669	elif Nsyms in [16,24]:
670	if '3' in OprNames[1]:
671	UsymOp.append('')
672	OprFlg.append(SGData['SGGen'][3])
673	for i in range(Nsyms):
674	if 'mx' in OprNames[i]:
675	UsymOp[-1] = OprNames[i]
676	elif 'm11' in OprNames[i]:
677	UsymOp[-1] = OprNames[i]
678	elif '211' in OprNames[i]:
679	UsymOp[-1] = OprNames[i]
680	OprFlg[-1] = 24
681	else: #4/mmm or 6/mmm
682	UsymOp.append(' mz ')
683	OprFlg.append(1)
684	if '4' in OprNames[1]: #4/mmm
685	UsymOp.append(' mx ')
686	OprFlg.append(20)
687	UsymOp.append(' m110 ')
688	OprFlg.append(24)
689	else: #6/mmm
690	UsymOp.append(' m110 ')
691	OprFlg.append(4)
692	UsymOp.append(' m+-0 ')
693	OprFlg.append(8)
694	else: #System is cubic
695	if Nsyms == 48:
696	UsymOp.append(' mx ')
697	OprFlg.append(4)
698	UsymOp.append(' m110 ')
699	OprFlg.append(24)
700
701	if 'P' in SGData['SGLatt']:
702	if SGData['SGSys'] == 'triclinic':
703	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}
704	elif SGData['SGSys'] == 'monoclinic':
705	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}
706	if SGData['SGUniq'] == 'a':
707	BNSsym.update({'P_B':[.5,0,.5],'P_C':[.5,.5,0]})
708	elif SGData['SGUniq'] == 'b':
709	BNSsym.update({'P_A':[.5,.5,0],'P_C':[0,.5,.5]})
710	elif SGData['SGUniq'] == 'c':
711	BNSsym.update({'P_A':[0,.5,.5],'P_B':[.5,0,.5]})
712	elif SGData['SGSys'] == 'orthorhombic':
713	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5],
714	'P_A':[0,.5,.5],'P_B':[.5,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
715	elif SGData['SGSys'] == 'tetragonal':
716	BNSsym = {'P_c':[0,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
717	elif SGData['SGSys'] in ['trigonal','hexagonal']:
718	BNSsym = {'P_c':[0,0,.5]}
719	elif SGData['SGSys'] == 'cubic':
720	BNSsym = {'P_I':[.5,.5,.5]}
721
722	elif 'A' in SGData['SGLatt']:
723	if SGData['SGSys'] == 'monoclinic':
724	BNSsym = {}
725	if SGData['SGUniq'] == 'b':
726	BNSsym.update({'A_a':[.5,0,0],'A_c':[0,0,.5]})
727	elif SGData['SGUniq'] == 'c':
728	BNSsym.update({'A_a':[.5,0,0],'A_b':[0,.5,0]})
729	elif SGData['SGSys'] == 'orthorhombic':
730	BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5],
731	'A_B':[.5,0,.5],'A_C':[.5,.5,0]}
732	elif SGData['SGSys'] == 'triclinic':
733	BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5]}
734
735	elif 'B' in SGData['SGLatt']:
736	if SGData['SGSys'] == 'monoclinic':
737	BNSsym = {}
738	if SGData['SGUniq'] == 'a':
739	BNSsym.update({'B_b':[0,.5,0],'B_c':[0,0,.5]})
740	elif SGData['SGUniq'] == 'c':
741	BNSsym.update({'B_a':[.5,0,0],'B_b':[0,.5,0]})
742	elif SGData['SGSys'] == 'orthorhombic':
743	BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5],
744	'B_A':[0,.5,.5],'B_C':[.5,.5,0]}
745	elif SGData['SGSys'] == 'triclinic':
746	BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5]}
747
748	elif 'C' in SGData['SGLatt']:
749	if SGData['SGSys'] == 'monoclinic':
750	BNSsym = {}
751	if SGData['SGUniq'] == 'a':
752	BNSsym.update({'C_b':[0,.5,.0],'C_c':[0,0,.5]})
753	elif SGData['SGUniq'] == 'b':
754	BNSsym.update({'C_a':[.5,0,0],'C_c':[0,0,.5]})
755	elif SGData['SGSys'] == 'orthorhombic':
756	BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5],
757	'C_A':[0,.5,.5],'C_B':[.5,0,.5]}
758	elif SGData['SGSys'] == 'triclinic':
759	BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5]}
760
761	elif 'I' in SGData['SGLatt']:
762	if SGData['SGSys'] in ['monoclinic','orthorhombic','triclinic']:
763	BNSsym = {'I_a':[.5,0,0],'I_b':[0,.5,0],'I_c':[0,0,.5]}
764	elif SGData['SGSys'] == 'tetragonal':
765	BNSsym = {'I_c':[0,0,.5]}
766	elif SGData['SGSys'] == 'cubic':
767	BNSsym = {}
768
769	elif 'F' in SGData['SGLatt']:
770	if SGData['SGSys'] in ['monoclinic','orthorhombic','cubic','triclinic']:
771	BNSsym = {'F_S':[.5,.5,.5]}
772
773	elif 'R' in SGData['SGLatt']:
774	BNSsym = {'R_I':[0,0,.5]}
775	return UsymOp,OprFlg,BNSsym
776
777	def ApplyBNSlatt(SGData,BNSlatt):
778	Tmat = np.eye(3)
779	BNS = BNSlatt[0]
780	A = np.array(BNSlatt[1])
781	SGCen = SGData['SGCen']
782	if '_a' in BNS:
783	Tmat[0,0] = 2.0
784	elif '_b' in BNS:
785	Tmat[1,1] = 2.0
786	elif '_c' in BNS:
787	Tmat[2,2] = 2.0
788	elif '_A' in BNS:
789	Tmat[0,0] = 2.0
790	elif '_B' in BNS:
791	Tmat[1,1] = 2.0
792	elif '_C' in BNS:
793	Tmat[2,2] = 2.0
794	elif '_I' in BNS:
795	Tmat *= 2.0
796	SGData['SGSpin'][-1] = -1
797	if 'R' in BNS:
798	SGData['SGSpin'][-1] = -1
799	elif '_S' in BNS:
800	SGData['SGSpin'][-1] = -1
801	SGData['SGSpin'] += [-1,-1,-1,]
802	Tmat *= 2.0
803	SGData['SGSpin'].append(-1)
804	if 'P' not in BNS:
805	SGData['SGSpin'].append(-1)
806	C = SGCen+A
807	SGData['SGCen'] = np.vstack((SGCen,C))%1.
808	return Tmat
809
810	def CheckSpin(isym,SGData):
811	''' Check for exceptions in spin rules
812	'''
813	if SGData['SGPtGrp'] in ['222','mm2','2mm','m2m']: #only 2/3 can be red; not 1/3 or 3/3
814	if SGData['SGSpin'][1]SGData['SGSpin'][2]SGData['SGSpin'][3] < 0:
815	SGData['SGSpin'][(isym+1)%3+1] *= -1
816	if SGData['SpGrp'][0] == 'F' and isym > 2:
817	SGData['SGSpin'][(isym+1)%3+3] == 1
818	elif SGData['SGPtGrp'] == 'mmm':
819	if SGData['SpGrp'][0] == 'F' and isym > 2:
820	SGData['SGSpin'][(isym+1)%3+3] == 1
821
822	def MagSGSym(SGData): #needs to use SGPtGrp not SGLaue!
823	SGLaue = SGData['SGLaue']
824	if '1' not in SGData['GenSym']: #patch for old gpx files
825	SGData['GenSym'] = ['1',]+SGData['GenSym']
826	SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])
827	if len(SGData['SGSpin'])<len(SGData['GenSym']):
828	SGData['SGSpin'] = [1,]+list(SGData['SGSpin']) #end patch
829	GenSym = SGData['GenSym'][1:] #skip identity
830	SpnFlp = SGData['SGSpin']
831	# print('SpnFlp',SpnFlp)
832	SGPtGrp = SGData['SGPtGrp']
833	if len(SpnFlp) == 1:
834	SGData['MagPtGp'] = SGPtGrp
835	return SGData['SpGrp']
836	magSym = SGData['SpGrp'].split()
837	if SGLaue in ['-1',]:
838	SGData['MagPtGp'] = SGPtGrp
839	if SpnFlp[1] == -1:
840	magSym[1] += "'"
841	SGData['MagPtGp'] += "'"
842	elif SGLaue in ['2/m','4/m','6/m']: #all ok
843	Uniq = {'a':1,'b':2,'c':3,'':1}
844	id = [0,1]
845	if len(magSym) > 2:
846	id = [0,Uniq[SGData['SGUniq']]]
847	sym = magSym[id[1]].split('/')
848	Ptsym = SGLaue.split('/')
849	if len(GenSym) == 3:
850	for i in [0,1,2]:
851	if SpnFlp[i+1] < 0:
852	sym[i] += "'"
853	Ptsym[i] += "'"
854	else:
855	for i in range(len(GenSym)):
856	if SpnFlp[i+1] < 0:
857	sym[i] += "'"
858	Ptsym[i] += "'"
859	SGData['MagPtGp'] = '/'.join(Ptsym)
860	magSym[id[1]] = '/'.join(sym)
861	elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
862	SGData['MagPtGp'] = ''
863	for i in [0,1,2]:
864	SGData['MagPtGp'] += SGPtGrp[i]
865	if SpnFlp[i+1] < 0:
866	magSym[i+1] += "'"
867	SGData['MagPtGp'] += "'"
868	elif SGLaue == '6/mmm': #ok
869	magPtGp = list(SGPtGrp)
870	if len(GenSym) == 2:
871	for i in [0,1]:
872	if SpnFlp[i+1] < 0:
873	magSym[i+2] += "'"
874	magPtGp[i+1] += "'"
875	if SpnFlp[1]*SpnFlp[2] < 0:
876	magSym[1] += "'"
877	magPtGp[0] += "'"
878	else:
879	sym = magSym[1].split('/')
880	Ptsym = ['6','m']
881	magPtGp = ['','m','m']
882	for i in [0,1,2]:
883	if SpnFlp[i+1] < 0:
884	if i:
885	magSym[i+1] += "'"
886	magPtGp[i] += "'"
887	else:
888	sym[1] += "'"
889	Ptsym[0] += "'"
890	if SpnFlp[2]*SpnFlp[3] < 0:
891	sym[0] += "'"
892	Ptsym[0] += "'"
893	magSym[1] = '/'.join(sym)
894	magPtGp[0] = '/'.join(Ptsym)
895	SGData['MagPtGp'] = ''.join(magPtGp)
896	elif SGLaue == '4/mmm':
897	magPtGp = list(SGPtGrp)
898	if len(GenSym) == 2:
899	for i in [0,1]:
900	if SpnFlp[i+1] < 0:
901	magSym[i+2] += "'"
902	magPtGp[i+1] += "'"
903	if SpnFlp[1]*SpnFlp[2] < 0:
904	magSym[1] += "'"
905	magPtGp[0] += "'"
906	else:
907	if '/' in magSym[1]: #P 4/m m m, etc.
908	sym = magSym[1].split('/')
909	Ptsym = ['4','m']
910	magPtGp = ['','m','m']
911	for i in [0,1,2]:
912	if SpnFlp[i+1] < 0:
913	if i:
914	magSym[i+1] += "'"
915	magPtGp[i] += "'"
916	else:
917	sym[1] += "'"
918	Ptsym[1] += "'"
919	if SpnFlp[2]*SpnFlp[3] < 0:
920	sym[0] += "'"
921	Ptsym[0] += "'"
922	magSym[1] = '/'.join(sym)
923	magPtGp[0] = '/'.join(Ptsym)
924	else:
925	for i in [0,1]:
926	if SpnFlp[i+1] < 0:
927	magSym[i+2] += "'"
928	if SpnFlp[1]*SpnFlp[2] < 0:
929	magSym[1] += "'"
930	SGData['MagPtGp'] = ''.join(magPtGp)
931	elif SGLaue in ['3','3m1','31m']: #ok
932	Ptsym = list(SGLaue)
933	if len(GenSym) == 1: #all ok
934	id = 2
935	if (len(magSym) == 4) and (magSym[2] == '1'):
936	id = 3
937	if '3' in GenSym[0]:
938	id = 1
939	magSym[id].strip("'")
940	if SpnFlp[1] < 0:
941	magSym[id] += "'"
942	Ptsym[id-1] += "'"
943	elif len(GenSym) == 2:
944	if 'R' in GenSym[1]:
945	magSym[-1].strip("'")
946	if SpnFlp[1] < 0:
947	magSym[-1] += "'"
948	Ptsym[-1] += "'"
949	else:
950	i,j = [1,2]
951	if magSym[2] == '1':
952	i,j = [1,3]
953	magSym[i].strip("'")
954	Ptsym[i-1].strip("'")
955	magSym[j].strip("'")
956	Ptsym[j-1].strip("'")
957	if SpnFlp[1:3] == [1,-1]:
958	magSym[i] += "'"
959	Ptsym[i-1] += "'"
960	elif SpnFlp[1:3] == [-1,-1]:
961	magSym[j] += "'"
962	Ptsym[j-1] += "'"
963	elif SpnFlp[1:3] == [-1,1]:
964	magSym[i] += "'"
965	Ptsym[i-1] += "'"
966	magSym[j] += "'"
967	Ptsym[j-1] += "'"
968	elif len(GenSym):
969	if 'c' not in magSym[2]:
970	i,j = [1,2]
971	magSym[i].strip("'")
972	Ptsym[i-1].strip("'")
973	magSym[j].strip("'")
974	Ptsym[j-1].strip("'")
975	if SpnFlp[1:3] == [1,-1]:
976	magSym[i] += "'"
977	Ptsym[i-1] += "'"
978	elif SpnFlp[1:3] == [-1,-1]:
979	magSym[j] += "'"
980	Ptsym[j-1] += "'"
981	elif SpnFlp[2] == [-1,1]:
982	magSym[i] += "'"
983	Ptsym[i-1] += "'"
984	magSym[j] += "'"
985	Ptsym[j-1] += "'"
986	SGData['MagPtGp'] = ''.join(Ptsym)
987	elif SGData['SGPtGrp'] == '23' and len(magSym):
988	SGData['MagPtGp'] = '23'
989	elif SGData['SGPtGrp'] == 'm3':
990	SGData['MagPtGp'] = "m3"
991	if SpnFlp[1] < 0:
992	magSym[1] += "'"
993	magSym[2] += "'"
994	SGData['MagPtGp'] = "m'3'"
995	if SpnFlp[2] < 0:
996	if not 'm' in magSym[1]: #only Ia3
997	magSym[1].strip("'")
998	SGData['MagPtGp'] = "m3'"
999	elif SGData['SGPtGrp'] in ['432','-43m']:
1000	Ptsym = SGData['SGPtGrp'].split('3')
1001	if SpnFlp[1] < 0:
1002	magSym[1] += "'"
1003	Ptsym[0] += "'"
1004	magSym[3] += "'"
1005	Ptsym[1] += "'"
1006	SGData['MagPtGp'] = '3'.join(Ptsym)
1007	elif SGData['SGPtGrp'] == 'm3m':
1008	Ptsym = ['m','3','m']
1009	if SpnFlp[1:3] == [-1,1]:
1010	magSym[1] += "'"
1011	Ptsym[0] += "'"
1012	magSym[2] += "'"
1013	Ptsym[1] += "'"
1014	elif SpnFlp[1:3] == [1,-1]:
1015	magSym[3] += "'"
1016	Ptsym[2] += "'"
1017	elif SpnFlp[1:3] == [-1,-1]:
1018	magSym[1] += "'"
1019	Ptsym[0] += "'"
1020	magSym[2] += "'"
1021	Ptsym[1] += "'"
1022	magSym[3] += "'"
1023	Ptsym[2] += "'"
1024	SGData['MagPtGp'] = ''.join(Ptsym)
1025	# print SpnFlp
1026	magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
1027	return ' '.join(magSym)
1028
1029	def MagText2MTS(mcifOpr,CIF=True):
1030	"From magnetic space group cif text returns matrix/translation + spin flip"
1031	XYZ = {'x':[1,0,0],'+x':[1,0,0],'-x':[-1,0,0],'y':[0,1,0],'+y':[0,1,0],'-y':[0,-1,0],
1032	'z':[0,0,1],'+z':[0,0,1],'-z':[0,0,-1],'x-y':[1,-1,0],'-x+y':[-1,1,0],'y-x':[-1,1,0],
1033	'+x-y':[1,-1,0],'+y-x':[-1,1,0]}
1034	ops = mcifOpr.split(",")
1035	M = []
1036	T = []
1037	for op in ops[:3]:
1038	ip = len(op)
1039	if '/' in op:
1040	if CIF:
1041	opMT = op.split('+')
1042	T.append(eval(opMT[1]))
1043	else:
1044	ip = op.index('/')
1045	T.append(eval(op[:ip+2]))
1046	opMT = [op[ip+2:],'']
1047	else:
1048	opMT = [op,'']
1049	T.append(0.)
1050	M.append(XYZ[opMT[0].lower()])
1051	spnflp = 1
1052	if '-1' in ops[3]:
1053	spnflp = -1
1054	return np.array(M),np.array(T),spnflp
1055
1056	def MagSSText2MTS(mcifOpr):
1057	"From magnetic super space group cif text returns matrix/translation + spin flip"
1058	XYZ = {'x1':[1,0,0,0],'-x1':[-1,0,0,0],
1059	'x2':[0,1,0,0],'-x2':[0,-1,0,0],
1060	'x3':[0,0,1,0],'-x3':[0,0,-1,0],
1061	'x4':[0,0,0,1],'-x4':[0,0,0,-1],
1062	'x1-x2':[1,-1,0,0],'-x1+x2':[-1,1,0,0],
1063	'x1-x4':[1,0,0,-1],'-x1+x4':[-1,0,0,1],
1064	'x2-x4':[0,1,0,-1],'-x2+x4':[0,-1,0,1],
1065	'-x1-x2+x4':[-1,-1,0,1],'x1+x2-x4':[1,1,0,-1]}
1066	ops = mcifOpr.split(",")
1067	M = []
1068	T = []
1069	for op in ops[:4]:
1070	ip = len(op)
1071	if '/' in op:
1072	ip = op.index('/')-2
1073	T.append(eval(op[ip:]))
1074	else:
1075	T.append(0.)
1076	M.append(XYZ[op[:ip]])
1077	spnflp = 1
1078	if '-1' in ops[4]:
1079	spnflp = -1
1080	return np.array(M),np.array(T),spnflp
1081
1082	def GenMagOps(SGData):
1083	FlpSpn = SGData['SGSpin']
1084	Nsym = len(SGData['SGOps'])
1085	Ncv = len(SGData['SGCen'])
1086	sgOp = [M for M,T in SGData['SGOps']]
1087	oprName = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
1088	if SGData['SGInv'] and not SGData['SGFixed'] and not SGData['SGGray']:
1089	Nsym *= 2
1090	sgOp += [-M for M,T in SGData['SGOps']]
1091	oprName += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
1092	Nsyms = 0
1093	sgOps = []
1094	OprNames = []
1095	for incv in range(Ncv):
1096	Nsyms += Nsym
1097	sgOps += sgOp
1098	OprNames += oprName
1099	if SGData['SGFixed']:
1100	SpnFlp = SGData['SpnFlp']
1101	else:
1102	SpnFlp = np.ones(Nsym,dtype=np.int)
1103	GenFlg = SGData.get('GenFlg',[0])
1104	Ngen = len(SGData['SGGen'])
1105	# print ('GenFlg:',SGData['GenFlg'])
1106	# print ('GenSym:',SGData['GenSym'])
1107	Nfl = len(GenFlg)
1108	for ieqv in range(Nsym):
1109	for iunq in range(Nfl):
1110	if SGData['SGGen'][ieqv%Ngen] & GenFlg[iunq]:
1111	SpnFlp[ieqv] *= FlpSpn[iunq]
1112	# print ('\nMagSpGrp:',SGData['MagSpGrp'],Ncv)
1113	# print ('FlpSpn:',Nfl,FlpSpn)
1114	for incv in range(Ncv):
1115	if incv:
1116	try:
1117	SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))
1118	except IndexError:
1119	FlpSpn = [1,]+FlpSpn
1120	SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))
1121	detM = [nl.det(M) for M in sgOp]
1122	MagMom = SpnFlpnp.array(NcvdetM) #duplicate for no. centerings
1123	SGData['MagMom'] = MagMom
1124	# print ('SgOps:',OprNames)
1125	# print ('SGGen:',SGData['SGGen'])
1126	# print ('SpnFlp:',SpnFlp)
1127	# print ('MagMom:',MagMom)
1128	return OprNames,SpnFlp
1129
1130	def GetOpNum(Opr,SGData):
1131	Nops = len(SGData['SGOps'])
1132	opNum = abs(Opr)%100
1133	cent = abs(Opr)//100
1134	if Opr < 0 and not SGData['SGFixed']:
1135	opNum += Nops
1136	if SGData['SGInv'] and not SGData['SGFixed']:
1137	Nops *= 2
1138	opNum += cent*Nops
1139	return opNum
1140
1141	################################################################################
1142	#### Superspace group codes
1143	################################################################################
1144
1145	def SSpcGroup(SGData,SSymbol):
1146	"""
1147	Determines supersymmetry information from superspace group name; currently only for (3+1) superlattices
1148
1149	:param SGData: space group data structure as defined in SpcGroup above (see :ref:`SGData<SGData_table>`).
1150	:param SSymbol: superspace group symbol extension (string) defining modulation direction & generator info.
1151	:returns: (SSGError,SSGData)
1152
1153	* SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
1154	* SSGData - is a dict (see :ref:`Superspace Group object<SSGData_table>`) with entries:
1155
1156	* 'SSpGrp': full superspace group symbol, accidental spaces removed; for display only
1157	* 'SSGCen': 4D cell centering vectors [0,0,0,0] at least
1158	* 'SSGOps': 4D symmetry operations as [M,T] so that M*x+T = x'
1159
1160	"""
1161
1162	def fixMonoOrtho():
1163	mod = ''.join(modsym).replace('1/2','0').replace('1','0')
1164	if SGData['SGPtGrp'] in ['2','m']: #OK
1165	if mod in ['a00','0b0','00g']:
1166	result = [i*-1 for i in SGData['SSGKl']]
1167	else:
1168	result = SGData['SSGKl'][:]
1169	if '/' in mod:
1170	return [i*-1 for i in result]
1171	else:
1172	return result
1173	elif SGData['SGPtGrp'] == '2/m': #OK
1174	if mod in ['a00','0b0','00g']:
1175	result = SGData['SSGKl'][:]
1176	else:
1177	result = [i*-1 for i in SGData['SSGKl']]
1178	if '/' in mod:
1179	return [i*-1 for i in result]
1180	else:
1181	return result
1182	else: #orthorhombic
1183	return [-SSGKl[i] if mod[i] in ['a','b','g'] else SSGKl[i] for i in range(3)]
1184
1185	def extendSSGOps(SSGOps):
1186	for OpA in SSGOps:
1187	OpAtxt = SSMT2text(OpA)
1188	if 't' not in OpAtxt:
1189	continue
1190	for OpB in SSGOps:
1191	OpBtxt = SSMT2text(OpB)
1192	if 't' not in OpBtxt:
1193	continue
1194	OpC = list(SGProd(OpB,OpA))
1195	OpC[1] %= 1.
1196	OpCtxt = SSMT2text(OpC)
1197	# print OpAtxt.replace(' ','')+' * '+OpBtxt.replace(' ','')+' = '+OpCtxt.replace(' ','')
1198	for k,OpD in enumerate(SSGOps):
1199	OpDtxt = SSMT2text(OpD)
1200	OpDtxt2 = ''
1201	if SGData['SGGray']:
1202	OpDtxt2 = SSMT2text([OpD[0],OpD[1]+np.array([0.,0.,0.,.5])])
1203	# print ' ('+OpCtxt.replace(' ','')+' = ? '+OpDtxt.replace(' ','')+')'
1204	if OpCtxt == OpDtxt:
1205	continue
1206	elif OpCtxt == OpDtxt2:
1207	continue
1208	elif OpCtxt.split(',')[:3] == OpDtxt.split(',')[:3]:
1209	if 't' not in OpDtxt:
1210	SSGOps[k] = OpC
1211	# print k,' new:',OpCtxt.replace(' ','')
1212	break
1213	else:
1214	OpCtxt = OpCtxt.replace(' ','')
1215	OpDtxt = OpDtxt.replace(' ','')
1216	Txt = OpCtxt+' conflicts with '+OpDtxt
1217	# print (Txt)
1218	return False,Txt
1219	return True,SSGOps
1220
1221	def findMod(modSym):
1222	for a in ['a','b','g']:
1223	if a in modSym:
1224	return a
1225
1226	def genSSGOps():
1227	SSGOps = SSGData['SSGOps'][:]
1228	iFrac = {}
1229	for i,frac in enumerate(SSGData['modSymb']):
1230	if frac in ['1/2','1/3','1/4','1/6','1']:
1231	iFrac[i] = frac+'.'
1232	# print SGData['SpGrp']+SSymbol
1233	# print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac,'modSymb',SSGData['modSymb']
1234	# set identity & 1,-1; triclinic
1235	SSGOps[0][0][3,3] = 1.
1236	## expand if centrosymmetric
1237	# if SGData['SGInv']:
1238	# SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
1239	# monoclinic - all done & all checked
1240	if SGData['SGPtGrp'] in ['2','m']: #OK
1241	SSGOps[1][0][3,3] = SSGKl[0]
1242	SSGOps[1][1][3] = genQ[0]
1243	for i in iFrac:
1244	SSGOps[1][0][3,i] = -SSGKl[0]
1245	elif SGData['SGPtGrp'] == '2/m': #OK
1246	SSGOps[1][0][3,3] = SSGKl[1]
1247	if gensym:
1248	SSGOps[1][1][3] = 0.5
1249	for i in iFrac:
1250	SSGOps[1][0][3,i] = SSGKl[0]
1251
1252	# orthorhombic - all OK not fully checked
1253	elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm']: #OK
1254	if SGData['SGPtGrp'] == '222':
1255	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[1,2],2:[1,3]},'b':{2:[3,2],0:[1,2]}} #OK
1256	elif SGData['SGPtGrp'] == 'mm2':
1257	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} #OK
1258	elif SGData['SGPtGrp'] == 'm2m':
1259	OrOps = {'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]}} #OK
1260	elif SGData['SGPtGrp'] == '2mm':
1261	OrOps = {'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]}} #OK
1262	a = findMod(SSGData['modSymb'])
1263	OrFrac = OrOps[a]
1264	for j in iFrac:
1265	for i in OrFrac[j]:
1266	SSGOps[i][0][3,j] = -2.eval(iFrac[j])SSGKl[i-1]
1267	for i in [0,1,2]:
1268	SSGOps[i+1][0][3,3] = SSGKl[i]
1269	SSGOps[i+1][1][3] = genQ[i]
1270	E,SSGOps = extendSSGOps(SSGOps)
1271	if not E:
1272	return E,SSGOps
1273	elif SGData['SGPtGrp'] == 'mmm': #OK
1274	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}}
1275	a = findMod(SSGData['modSymb'])
1276	if a == 'g':
1277	SSkl = [1,1,1]
1278	elif a == 'a':
1279	SSkl = [-1,1,-1]
1280	else:
1281	SSkl = [1,-1,-1]
1282	OrFrac = OrOps[a]
1283	for j in iFrac:
1284	for i in OrFrac[j]:
1285	SSGOps[i][0][3,j] = -2.eval(iFrac[j])SSkl[i-1]
1286	for i in [0,1,2]:
1287	SSGOps[i+1][0][3,3] = SSkl[i]
1288	SSGOps[i+1][1][3] = genQ[i]
1289	E,SSGOps = extendSSGOps(SSGOps)
1290	if not E:
1291	return E,SSGOps
1292	# tetragonal - all done & checked
1293	elif SGData['SGPtGrp'] == '4': #OK
1294	SSGOps[1][0][3,3] = SSGKl[0]
1295	SSGOps[1][1][3] = genQ[0]
1296	if '1/2' in SSGData['modSymb']:
1297	SSGOps[1][0][3,1] = -1
1298	elif SGData['SGPtGrp'] == '-4': #OK
1299	SSGOps[1][0][3,3] = SSGKl[0]
1300	if '1/2' in SSGData['modSymb']:
1301	SSGOps[1][0][3,1] = 1
1302	elif SGData['SGPtGrp'] in ['4/m',]: #OK
1303	if '1/2' in SSGData['modSymb']:
1304	SSGOps[1][0][3,1] = -SSGKl[0]
1305	for i,j in enumerate([1,3]):
1306	SSGOps[j][0][3,3] = 1
1307	if genQ[i]:
1308	SSGOps[j][1][3] = genQ[i]
1309	E,SSGOps = extendSSGOps(SSGOps)
1310	if not E:
1311	return E,SSGOps
1312	elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]: #OK
1313	iGens = [1,4,5]
1314	if SGData['SGPtGrp'] in ['4mm','-4m2',]:
1315	iGens = [1,6,7]
1316	for i,j in enumerate(iGens):
1317	if '1/2' in SSGData['modSymb'] and i < 2:
1318	SSGOps[j][0][3,1] = SSGKl[i]
1319	SSGOps[j][0][3,3] = SSGKl[i]
1320	if genQ[i]:
1321	if 's' in gensym and j == 6:
1322	SSGOps[j][1][3] = -genQ[i]
1323	else:
1324	SSGOps[j][1][3] = genQ[i]
1325	E,SSGOps = extendSSGOps(SSGOps)
1326	if not E:
1327	return E,SSGOps
1328	elif SGData['SGPtGrp'] in ['4/mmm',]:#OK
1329	if '1/2' in SSGData['modSymb']:
1330	SSGOps[1][0][3,1] = -SSGKl[0]
1331	SSGOps[6][0][3,1] = SSGKl[1]
1332	if modsym:
1333	SSGOps[1][1][3] = -genQ[3]
1334	for i,j in enumerate([1,2,6,7]):
1335	SSGOps[j][0][3,3] = 1
1336	SSGOps[j][1][3] = genQ[i]
1337	E,Result = extendSSGOps(SSGOps)
1338	if not E:
1339	return E,Result
1340	else:
1341	SSGOps = Result
1342
1343	# trigonal - all done & checked
1344	elif SGData['SGPtGrp'] == '3': #OK
1345	SSGOps[1][0][3,3] = SSGKl[0]
1346	if '1/3' in SSGData['modSymb']:
1347	SSGOps[1][0][3,1] = -1
1348	SSGOps[1][1][3] = genQ[0]
1349	elif SGData['SGPtGrp'] == '-3': #OK
1350	SSGOps[1][0][3,3] = -SSGKl[0]
1351	if '1/3' in SSGData['modSymb']:
1352	SSGOps[1][0][3,1] = -1
1353	SSGOps[1][1][3] = genQ[0]
1354	elif SGData['SGPtGrp'] in ['312','3m','-3m','-3m1','3m1']: #OK
1355	if '1/3' in SSGData['modSymb']:
1356	SSGOps[1][0][3,1] = -1
1357	for i,j in enumerate([1,5]):
1358	if SGData['SGPtGrp'] in ['3m','-3m']:
1359	SSGOps[j][0][3,3] = 1
1360	else:
1361	SSGOps[j][0][3,3] = SSGKl[i+1]
1362	if genQ[i]:
1363	SSGOps[j][1][3] = genQ[i]
1364	elif SGData['SGPtGrp'] in ['321','32']: #OK
1365	for i,j in enumerate([1,4]):
1366	SSGOps[j][0][3,3] = SSGKl[i]
1367	if genQ[i]:
1368	SSGOps[j][1][3] = genQ[i]
1369	elif SGData['SGPtGrp'] in ['31m','-31m']: #OK
1370	ids = [1,3]
1371	if SGData['SGPtGrp'] == '-31m':
1372	ids = [1,3]
1373	if '1/3' in SSGData['modSymb']:
1374	SSGOps[ids[0]][0][3,1] = -SSGKl[0]
1375	for i,j in enumerate(ids):
1376	SSGOps[j][0][3,3] = 1
1377	if genQ[i+1]:
1378	SSGOps[j][1][3] = genQ[i+1]
1379
1380	# hexagonal all done & checked
1381	elif SGData['SGPtGrp'] == '6': #OK
1382	SSGOps[1][0][3,3] = SSGKl[0]
1383	SSGOps[1][1][3] = genQ[0]
1384	elif SGData['SGPtGrp'] == '-6': #OK
1385	SSGOps[1][0][3,3] = SSGKl[0]
1386	elif SGData['SGPtGrp'] in ['6/m',]: #OK
1387	SSGOps[1][0][3,3] = -SSGKl[1]
1388	SSGOps[1][1][3] = genQ[0]
1389	SSGOps[2][1][3] = genQ[1]
1390	elif SGData['SGPtGrp'] in ['622',]: #OK
1391	for i,j in enumerate([1,8,9]):
1392	SSGOps[j][0][3,3] = SSGKl[i]
1393	if genQ[i]:
1394	SSGOps[j][1][3] = genQ[i]
1395	E,SSGOps = extendSSGOps(SSGOps)
1396
1397	elif SGData['SGPtGrp'] in ['6mm','-62m','-6m2',]: #OK
1398	for i,j in enumerate([1,6,7]):
1399	SSGOps[j][0][3,3] = SSGKl[i]
1400	if genQ[i]:
1401	SSGOps[j][1][3] = genQ[i]
1402	E,SSGOps = extendSSGOps(SSGOps)
1403	elif SGData['SGPtGrp'] in ['6/mmm',]: # OK
1404	for i,j in enumerate([1,2,10,11]):
1405	SSGOps[j][0][3,3] = 1
1406	if genQ[i]:
1407	SSGOps[j][1][3] = genQ[i]
1408	E,SSGOps = extendSSGOps(SSGOps)
1409	elif SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
1410	return True,SSGOps
1411	E,SSGOps = extendSSGOps(SSGOps)
1412	return E,SSGOps
1413
1414	def specialGen(gensym,modsym):
1415	sym = ''.join(gensym)
1416	if SGData['SGPtGrp'] in ['2/m',] and 'n' in SGData['SpGrp']:
1417	if 's' in sym:
1418	gensym = 'ss'
1419	if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
1420	gensym = '0ss'
1421	elif SGData['SGPtGrp'] in ['222',]:
1422	if sym == '00s':
1423	gensym = '0ss'
1424	elif sym == '0s0':
1425	gensym = 'ss0'
1426	elif sym == 's00':
1427	gensym = 's0s'
1428	elif SGData['SGPtGrp'] in ['mmm',]:
1429	if 'g' in modsym:
1430	if sym == 's00':
1431	gensym = 's0s'
1432	elif sym == '0s0':
1433	gensym = '0ss'
1434	elif 'a' in modsym:
1435	if sym == '0s0':
1436	gensym = 'ss0'
1437	elif sym == '00s':
1438	gensym = 's0s'
1439	elif 'b' in modsym:
1440	if sym == '00s':
1441	gensym = '0ss'
1442	elif sym == 's00':
1443	gensym = 'ss0'
1444	return gensym
1445
1446	Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'s':.5,'t':1./3,'q':.25,'h':1./6,'a':0.,'b':0.,'g':0.}
1447	if SGData['SGLaue'] in ['m3','m3m']:
1448	return '(3+1) superlattices not defined for cubic space groups',None
1449	elif SGData['SGLaue'] in ['3R','3mR']:
1450	return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
1451	if SGData['SGGray'] and SSymbol[-1] == 's':
1452	SSymbol = SSymbol[:-1]
1453	try:
1454	modsym,gensym = splitSSsym(SSymbol)
1455	except ValueError:
1456	return 'Error in superspace symbol '+SSymbol,None
1457	modQ = [Fracs[mod] for mod in modsym]
1458	SSGKl = SGData['SSGKl'][:]
1459	if SGData['SGLaue'] in ['2/m','mmm']:
1460	SSGKl = fixMonoOrtho()
1461	Ngen = len(gensym)
1462	if len(gensym) and Ngen != len(SSGKl):
1463	return 'Wrong number of items in generator symbol '+''.join(gensym),None
1464	gensym = specialGen(gensym[:Ngen],modsym)
1465	genQ = [Fracs[mod] for mod in gensym[:Ngen]]
1466	if not genQ:
1467	genQ = [0,0,0,0]
1468	SSgSpc = SGData['SpGrp']+SSymbol
1469	if SGData['SGGray']:
1470	SSgSpc = SSgSpc.replace('('," 1'(")
1471	SSgSpc += 's'
1472	SSGData = {'SSpGrp':SSgSpc,'modQ':modQ,'modSymb':modsym,'SSGKl':SSGKl}
1473	SSCen = np.zeros((len(SGData['SGCen']),4))
1474	for icen,cen in enumerate(SGData['SGCen']):
1475	SSCen[icen,0:3] = cen
1476	SSCen[0] = np.zeros(4)
1477	SSGData['SSGCen'] = SSCen
1478	SSGData['SSGOps'] = []
1479	for iop,op in enumerate(SGData['SGOps']):
1480	T = np.zeros(4)
1481	ssop = np.zeros((4,4))
1482	ssop[:3,:3] = op[0]
1483	T[:3] = op[1]
1484	SSGData['SSGOps'].append([ssop,T])
1485	E,Result = genSSGOps()
1486	if E:
1487	SSGData['SSGOps'] = Result
1488	if DEBUG:
1489	print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1490	for Op in Result:
1491	print (SSMT2text(Op).replace(' ',''))
1492	if SGData['SGInv']:
1493	for Op in Result:
1494	Op = [-Op[0],-Op[1]%1.]
1495	print (SSMT2text(Op).replace(' ',''))
1496	return None,SSGData
1497	else:
1498	return Result+'\nOperator conflict - incorrect superspace symbol',None
1499
1500	def SSChoice(SGData):
1501	'''
1502	Gets the unique set of possible super space groups for a given space group
1503	'''
1504	laueSS = {'-1':['(abg)',],
1505	'2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1506	'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1507	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1508	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1509	'4/m':['(00g)','(1/21/2g)'],
1510	'4/mmm':['(00g)','(1/21/2g)'],
1511	'3':['(00g)','(1/31/3g)'],
1512	'3m1':['(00g)'],
1513	'31m':['(00g)','(1/31/3g)'],
1514	'6/m':['(00g)',],
1515	'6/mmm':['(00g)',]}
1516
1517	laueTS = {'-1':['',],
1518	'2/m':['','s','0s','ss','s0'],
1519	'mmm':['','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1520	'4/m':['','q','s','s0',],
1521	'4/mmm':['','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1522	'3':['','t'],
1523	'3m1':['','t0','0s','t00','0s0'],
1524	'31m':['','t00','0ss'],
1525	'6/m':['','h','t','s','s0'],
1526	'6/mmm':['','h00','t00','s00','ss0','0ss','s0s','s0s0','00ss','s00s','ss00','0ss0','0s0s']}
1527	laue = SGData['SGLaue']
1528	SSChoice = []
1529	for ax in laueSS[laue]:
1530	for sx in laueTS[laue]:
1531	SSChoice.append(ax+sx)
1532	ssChoice = []
1533	ssHash = []
1534	for item in SSChoice:
1535	E,SSG = SSpcGroup(SGData,item)
1536	if SSG:
1537	sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1538	if sshash not in ssHash:
1539	ssHash.append(sshash)
1540	ssChoice.append(item)
1541	return ssChoice
1542
1543	def fixGray(SGData,SSymbol):
1544	modsym,gensym = SSymbol.replace(' ','').split(')')
1545	modsym += ')'
1546	sgPtGp = SGData['SGPtGrp']
1547	if gensym:
1548	if sgPtGp in ['1','2','m','3','4','6'] and len(gensym) == 1:
1549	gensym += 's'
1550	elif sgPtGp in ['2/m','4/m','6/m'] and len(gensym) == 2:
1551	gensym += 's'
1552	elif sgPtGp in ['4/mmm','6/mmm'] and len(gensym) == 4:
1553	gensym += 's'
1554	elif len(gensym) == 3:
1555	gensym += 's'
1556	else:
1557	if sgPtGp in ['1','2','m','3','4','6']:
1558	gensym += '0s'
1559	elif sgPtGp in ['2/m','4/m','6/m']:
1560	gensym += '00s'
1561	elif sgPtGp in ['4/mmm','6/mmm']:
1562	gensym += '0000s'
1563	else:
1564	gensym += '000s'
1565	return modsym+gensym
1566
1567	def splitSSsym(SSymbol):
1568	'''
1569	Splits supersymmetry symbol into two lists of strings
1570	'''
1571	modsym,gensym = SSymbol.replace(' ','').split(')')
1572	modsym = modsym.replace(',','')
1573	if "1'" in modsym:
1574	gensym = gensym[:-1]
1575	modsym = modsym.replace("1'",'')
1576	if gensym in ['0','00','000','0000']: #get rid of extraneous symbols
1577	gensym = ''
1578	nfrac = modsym.count('/')
1579	modsym = modsym.lstrip('(')
1580	if nfrac == 0:
1581	modsym = list(modsym)
1582	elif nfrac == 1:
1583	pos = modsym.find('/')
1584	if pos == 1:
1585	modsym = [modsym[:3],modsym[3],modsym[4]]
1586	elif pos == 2:
1587	modsym = [modsym[0],modsym[1:4],modsym[4]]
1588	else:
1589	modsym = [modsym[0],modsym[1],modsym[2:]]
1590	else:
1591	lpos = modsym.find('/')
1592	rpos = modsym.rfind('/')
1593	if lpos == 1 and rpos == 4:
1594	modsym = [modsym[:3],modsym[3:6],modsym[6]]
1595	elif lpos == 1 and rpos == 5:
1596	modsym = [modsym[:3],modsym[3],modsym[4:]]
1597	else:
1598	modsym = [modsym[0],modsym[1:4],modsym[4:]]
1599	gensym = list(gensym)
1600	return modsym,gensym
1601
1602	def SSGPrint(SGData,SSGData):
1603	'''
1604	Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1605
1606	:param SGData: space group data structure as defined in SpcGroup above.
1607	:param SSGData: from :func:`SSpcGroup`
1608	:returns:
1609	SSGText - list of strings with the superspace group details
1610	SGTable - list of strings for each of the operations
1611	'''
1612	Mult = len(SSGData['SSGCen'])len(SSGData['SSGOps'])(int(SGData['SGInv'])+1)
1613	SSGText = []
1614	SSGText.append(' Superspace Group: '+SSGData['SSpGrp'])
1615	CentStr = 'centrosymmetric'
1616	if not SGData['SGInv']:
1617	CentStr = 'non'+CentStr
1618	if SGData['SGLatt'] in 'ABCIFR':
1619	SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1620	else:
1621	SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1622	SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1623	SGptGp = SGData['SGPtGrp']
1624	if SGData['SGGray']:
1625	SGptGp += "1'"
1626	SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1627	SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1628	SSGText.append(' Multiplicity of a general site is '+str(Mult))
1629	if SGData['SGUniq'] in ['a','b','c']:
1630	SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1631	if SGData['SGInv']:
1632	SSGText.append(' The inversion center is located at 0,0,0')
1633	if SGData['SGPolax']:
1634	SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1635	SSGText.append(' ')
1636	if len(SSGData['SSGCen']) > 1:
1637	SSGText.append(' The equivalent positions are:')
1638	SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
1639	else:
1640	SSGText.append(' The equivalent positions are:\n')
1641	SSGTable = []
1642	for i,Opr in enumerate(SSGData['SSGOps']):
1643	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1644	if SGData['SGGray']:
1645	SSGTable.append(" for 1'")
1646	for i,Opr in enumerate(SSGData['SSGOps']):
1647	Opr2 = [Opr[0],Opr[1]+np.array([0,0,0,.5])]
1648	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr2)))
1649	return SSGText,SSGTable
1650
1651	def SSGModCheck(Vec,modSymb,newMod=True):
1652	''' Checks modulation vector compatibility with supersymmetry space group symbol.
1653	if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1654	'''
1655	Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1656	modQ = [Fracs[mod] for mod in modSymb]
1657	if newMod:
1658	newVec = Vec
1659	if not np.any(Vec):
1660	newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1661	return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)], \
1662	[True if mod in ['a','b','g'] else False for mod in modSymb]
1663	else:
1664	return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1665
1666	def SSMT2text(Opr):
1667	"From superspace group matrix/translation operator returns text version"
1668	XYZS = ('x','y','z','t') #Stokes, Campbell & van Smaalen notation
1669	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1670	Fld = ''
1671	M,T = Opr
1672	for j in range(4):
1673	IJ = ''
1674	for k in range(4):
1675	txt = str(int(round(M[j][k])))
1676	txt = txt.replace('1',XYZS[k]).replace('0','')
1677	if '2' in txt:
1678	txt += XYZS[k]
1679	if IJ and M[j][k] > 0:
1680	IJ += '+'+txt
1681	else:
1682	IJ += txt
1683	IK = int(round(T[j]*12))%12
1684	if IK:
1685	if not IJ:
1686	break
1687	if IJ[0] == '-':
1688	Fld += (TRA[IK]+IJ).rjust(8)
1689	else:
1690	Fld += (TRA[IK]+'+'+IJ).rjust(8)
1691	else:
1692	Fld += IJ.rjust(8)
1693	if j != 3: Fld += ', '
1694	return Fld
1695
1696	def SSLatt2text(SSGCen):
1697	"Lattice centering vectors to text"
1698	lattTxt = ''
1699	lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1700	for vec in SSGCen:
1701	lattTxt += ' '
1702	for item in vec:
1703	lattTxt += '%s,'%(lattDir[int(item*12)])
1704	lattTxt = lattTxt.rstrip(',')
1705	lattTxt += ';'
1706	lattTxt = lattTxt.rstrip(';').lstrip(' ')
1707	return lattTxt
1708
1709	def SSpaceGroup(SGSymbol,SSymbol):
1710	'''
1711	Print the output of SSpcGroup in a nicely formatted way.
1712
1713	:param SGSymbol: space group symbol with spaces between axial fields.
1714	:param SSymbol: superspace group symbol extension (string).
1715	:returns: nothing
1716	'''
1717
1718	E,A = SpcGroup(SGSymbol)
1719	if E > 0:
1720	print (SGErrors(E))
1721	return
1722	E,B = SSpcGroup(A,SSymbol)
1723	if E > 0:
1724	print (E)
1725	return
1726	for l in SSGPrint(B):
1727	print (l)
1728
1729	def SGProd(OpA,OpB):
1730	'''
1731	Form space group operator product. OpA & OpB are [M,V] pairs;
1732	both must be of same dimension (3 or 4). Returns [M,V] pair
1733	'''
1734	A,U = OpA
1735	B,V = OpB
1736	M = np.inner(B,A.T)
1737	W = np.inner(B,U)+V
1738	return M,W
1739
1740	def MoveToUnitCell(xyz):
1741	'''
1742	Translates a set of coordinates so that all values are >=0 and < 1
1743
1744	:param xyz: a list or numpy array of fractional coordinates
1745	:returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1746	'''
1747	XYZ = (np.array(xyz)+10.)%1.
1748	cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1749	return XYZ,cell
1750
1751	def Opposite(XYZ,toler=0.0002):
1752	'''
1753	Gives opposite corner, edge or face of unit cell for position within tolerance.
1754	Result may be just outside the cell within tolerance
1755
1756	:param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1757	:param toler: unit cell fraction tolerance making opposite
1758	:returns:
1759	XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1760	'''
1761	perm3 = [[1,1,1],[0,1,1],[1,0,1],[1,1,0],[1,0,0],[0,1,0],[0,0,1],[0,0,0]]
1762	TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1763	perm = TB*perm3
1764	cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1765	D = dict(zip(cperm,perm))
1766	new = {}
1767	for key in D:
1768	new[key] = np.array(D[key])+np.array(XYZ)
1769	return new
1770
1771	def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1772	'''
1773	Generates the equivalent positions for a specified coordinate and space group
1774
1775	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
1776	:param SGData: from :func:`SpcGroup`
1777	:param All: True return all equivalent positions including duplicates;
1778	False return only unique positions
1779	:param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1780	:param Move: True move generated atom positions to be inside cell
1781	False do not move atoms
1782	:return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
1783
1784	* [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1785	* Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1786	* is centering operator number (1-4) and - is for inversion
1787	Cell = unit cell translations needed to put new positions inside cell
1788	[UijEquiv] - equivalent Uij; absent if no Uij given
1789	* +1/-1 for spin inversion of operator - empty if not magnetic
1790
1791	'''
1792	XYZEquiv = []
1793	UijEquiv = []
1794	Idup = []
1795	Cell = []
1796	inv = int(SGData['SGInv']+1)
1797	icen = SGData['SGCen']
1798	if SGData['SGFixed']:
1799	inv = 1
1800	SpnFlp = SGData.get('SpnFlp',[])
1801	spnflp = []
1802	X = np.array(XYZ)
1803	mj = 0
1804	if Move:
1805	X = MoveToUnitCell(X)[0]
1806	for ic,cen in enumerate(icen):
1807	C = np.array(cen)
1808	for invers in range(inv):
1809	for io,[M,T] in enumerate(SGData['SGOps']):
1810	idup = ((io+1)+100ic)(1-2*invers)
1811	XT = np.inner(M,X)+T
1812	if len(Uij):
1813	U = Uij2U(Uij)
1814	U = np.inner(M,np.inner(U,M).T)
1815	newUij = U2Uij(U)
1816	if invers:
1817	XT = -XT
1818	XT += C
1819	cell = np.zeros(3,dtype=np.int32)
1820	cellj = np.zeros(3,dtype=np.int32)
1821	if Move:
1822	newX,cellj = MoveToUnitCell(XT)
1823	else:
1824	newX = XT
1825	cell += cellj
1826	if All:
1827	if np.allclose(newX,X,atol=0.0002):
1828	idup = False
1829	else:
1830	if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1831	idup = False
1832	if All or idup:
1833	XYZEquiv.append(newX)
1834	Idup.append(idup)
1835	Cell.append(cell)
1836	if len(Uij):
1837	UijEquiv.append(newUij)
1838	if len(SpnFlp):
1839	spnflp.append(SpnFlp[mj])
1840	else:
1841	spnflp.append(1)
1842	mj += 1
1843	if len(Uij):
1844	return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
1845	else:
1846	return zip(XYZEquiv,Idup,Cell,spnflp)
1847
1848	def GenHKL(HKL,SGData):
1849	''' Generates all equivlent reflections including Friedel pairs
1850	:param HKL: [h,k,l] must be integral values
1851	:param SGData: space group data obtained from SpcGroup
1852	:returns: array Uniq: equivalent reflections
1853	'''
1854
1855	Ops = SGData['SGOps']
1856	OpM = np.array([op[0] for op in Ops])
1857	Uniq = np.inner(OpM,HKL)
1858	Uniq = list(Uniq)+list(-1*Uniq)
1859	return np.array(Uniq)
1860
1861	def GenHKLf(HKL,SGData):
1862	'''
1863	Uses old GSAS Fortran routine genhkl.for
1864
1865	:param HKL: [h,k,l] must be integral values for genhkl.for to work
1866	:param SGData: space group data obtained from SpcGroup
1867	:returns: iabsnt,mulp,Uniq,phi
1868
1869	* iabsnt = True if reflection is forbidden by symmetry
1870	* mulp = reflection multiplicity including Friedel pairs
1871	* Uniq = numpy array of equivalent hkl in descending order of h,k,l
1872	* phi = phase offset for each equivalent h,k,l
1873
1874	'''
1875	hklf = list(HKL)+[0,] #could be numpy array!
1876	Ops = SGData['SGOps']
1877	OpM = np.array([op[0] for op in Ops],order='F')
1878	OpT = np.array([op[1] for op in Ops])
1879	Cen = np.array([cen for cen in SGData['SGCen']],order='F')
1880
1881	import pyspg
1882	Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1883	h,k,l,f = Uniq
1884	Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
1885	phi = f[:Nuniq]
1886	return iabsnt,mulp,Uniq,phi
1887
1888	def MagHKLchk(HKL,SGData):
1889	SpnFlp = SGData['SpnFlp']
1890	print(HKL)
1891	Uniq = GenHKL(HKL,SGData)
1892
1893	def checkSSLaue(HKL,SGData,SSGData):
1894	#Laue check here - Toss HKL if outside unique Laue part
1895	h,k,l,m = HKL
1896	if SGData['SGLaue'] == '2/m':
1897	if SGData['SGUniq'] == 'a':
1898	if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
1899	return False
1900	elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
1901	return False
1902	else:
1903	return True
1904	elif SGData['SGUniq'] == 'b':
1905	if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
1906	return False
1907	elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
1908	return False
1909	else:
1910	return True
1911	elif SGData['SGUniq'] == 'c':
1912	if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
1913	return False
1914	elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
1915	return False
1916	else:
1917	return True
1918	elif SGData['SGLaue'] == 'mmm':
1919	if 'a' in SSGData['modSymb']:
1920	if h == 0 and m < 0:
1921	return False
1922	else:
1923	return True
1924	elif 'b' in SSGData['modSymb']:
1925	if k == 0 and m < 0:
1926	return False
1927	else:
1928	return True
1929	elif 'g' in SSGData['modSymb']:
1930	if l == 0 and m < 0:
1931	return False
1932	else:
1933	return True
1934	else: #tetragonal, trigonal, hexagonal (& triclinic?)
1935	if l == 0 and m < 0:
1936	return False
1937	else:
1938	return True
1939
1940
1941	def checkSSextc(HKL,SSGData):
1942	Ops = SSGData['SSGOps']
1943	OpM = np.array([op[0] for op in Ops])
1944	OpT = np.array([op[1] for op in Ops])
1945	HKLS = np.array([HKL,-HKL]) #Freidel's Law
1946	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
1947	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
1948	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
1949	if dhkl.any():
1950	continue
1951	else:
1952	if phkl%1.:
1953	return False
1954	return True
1955
1956	################################################################################
1957	#### Site symmetry tables
1958	################################################################################
1959
1960	OprName = {
1961	'-6643': ['-1',1],'6479' : ['2(z)',2],'-6479': ['m(z)',3],
1962	'6481' : ['m(y)',4],'-6481': ['2(y)',5],'6641' : ['m(x)',6],
1963	'-6641': ['2(x)',7],'6591' : ['m(+-0)',8],'-6591': ['2(+-0)',9],
1964	'6531' : ['m(110)',10],'-6531': ['2(110)',11],'6537' : ['4(z)',12],
1965	'-6537': ['-4(z)',13],'975' : ['3(111)',14],'6456' : ['3',15],
1966	'-489' : ['3(+--)',16],'483' : ['3(-+-)',17],'-969' : ['3(--+)',18],
1967	'819' : ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' : ['m(0+-)',21],
1968	'-2431': ['2(0+-)',22],'-657' : ['m(xz)',23],'657' : ['2(xz)',24],
1969	'1943' : ['-4(x)',25],'-1943': ['4(x)',26],'-2429': ['m(yz)',27],
1970	'2429' : ['2(yz)',28],'639' : ['-4(y)',29],'-639' : ['4(y)',30],
1971	'-6484': ['2(010)',4],'6484' : ['m(010)',5],'-6668': ['2(100)',6],
1972	'6668' : ['m(100)',7],'-6454': ['2(120)',18],'6454' : ['m(120)',19],
1973	'-6638': ['2(210)',20],'6638' : ['m(210)',21], #search in SytSym ends at m(210)
1974	'2223' : ['3(+++)2',39],
1975	'6538' : ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
1976	'2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489' :['-3(+--)1',45],
1977	'801' : ['4(y)1',46],'1945' : ['4(x)3',47],'-6585': ['-4(z)3 ',48],
1978	'6585' : ['4(z)3',49],'6584' : ['3(z)2',50],'6666' : ['6(z)5 ',51],
1979	'6643' : ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
1980	'-6666': ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
1981	'-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
1982	'969' :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
1983	'-2151':['-3(+--)2',64], }
1984
1985	KNsym = {
1986	'0' :' 1 ','1' :' -1 ','64' :' 2(x)','32' :' m(x)',
1987	'97' :' 2/m(x)','16' :' 2(y)','8' :' m(y)','25' :' 2/m(y)',
1988	'2' :' 2(z)','4' :' m(z)','7' :' 2/m(z)','134217728' :' 2(yz)',
1989	'67108864' :' m(yz)','201326593' :' 2/m(yz)','2097152' :' 2(0+-)','1048576' :' m(0+-)',
1990	'3145729' :'2/m(0+-)','8388608' :' 2(xz)','4194304' :' m(xz)','12582913' :' 2/m(xz)',
1991	'524288' :' 2(+0-)','262144' :' m(+0-)','796433' :'2/m(+0-)','1024' :' 2(xy)',
1992	'512' :' m(xy)','1537' :' 2/m(xy)','256' :' 2(+-0)','128' :' m(+-0)',
1993	'385' :'2/m(+-0)','76' :' mm2(x)','52' :' mm2(y)','42' :' mm2(z)',
1994	'135266336' :' mm2(yz)','69206048' :'mm2(0+-)','8650760' :' mm2(xz)','4718600' :'mm2(+0-)',
1995	'1156' :' mm2(xy)','772' :'mm2(+-0)','82' :' 222 ','136314944' :' 222(x)',
1996	'8912912' :' 222(y)','1282' :' 222(z)','127' :' mmm ','204472417' :' mmm(x)',
1997	'13369369' :' mmm(y)','1927' :' mmm(z)','33554496' :' 4(x)','16777280' :' -4(x)',
1998	'50331745' :'4/m(x)' ,'169869394' :'422(x)','84934738' :'-42m(x)','101711948' :'4mm(x)',
1999	'254804095' :'4/mmm(x)','536870928 ':' 4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2000	'545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2001	'2050' :' 4(z)','4098' :' -4(z)','6151' :'4/m(z)','3410' :'422(z)',
2002	'4818' :'-42m(z)','2730' :'4mm(z)','8191' :'4/mmm(z)','8192' :' 3(111)',
2003	'8193' :' -3(111)','2629888' :' 32(111)','1319040' :' 3m(111)','3940737' :'-3m(111)',
2004	'32768' :' 3(+--)','32769' :' -3(+--)','10519552' :' 32(+--)','5276160' :' 3m(+--)',
2005	'15762945' :'-3m(+--)','65536' :' 3(-+-)','65537' :' -3(-+-)','134808576' :' 32(-+-)',
2006	'67437056' :' 3m(-+-)','202180097' :'-3m(-+-)','131072' :' 3(--+)','131073' :' -3(--+)',
2007	'142737664' :' 32(--+)','71434368' :' 3m(--+)','214040961' :'-3m(--+)','237650' :' 23 ',
2008	'237695' :' m3 ','715894098' :' 432 ','358068946' :' -43m ','1073725439':' m3m ',
2009	'68157504' :' mm2(d100)','4456464' :' mm2(d010)','642' :' mm2(d001)','153092172' :'-4m2(x)',
2010	'277348404' :'-4m2(y)','5418' :'-4m2(z)','1075726335':' 6/mmm ','1074414420':'-6m2(100)',
2011	'1075070124':'-6m2(120)','1075069650':' 6mm ','1074414890':' 622 ','1073758215':' 6/m ',
2012	'1073758212':' -6 ','1073758210':' 6 ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2013	'1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2014	'1073758209':' -3 ','1073758208':' 3 ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2015	'1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2016	'1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2017	'1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2018	'1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2019	'1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2020	'1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2021	'1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2022	'1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2023	'1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':' -1 ',
2024	}
2025
2026	NXUPQsym = {
2027	'1' :(28,29,28,28),'-1' :( 1,29,28, 0),'2(x)' :(12,18,12,25),'m(x)' :(25,18,12,25),
2028	'2/m(x)' :( 1,18, 0,-1),'2(y)' :(13,17,13,24),'m(y)' :(24,17,13,24),'2/m(y)' :( 1,17, 0,-1),
2029	'2(z)' :(14,16,14,23),'m(z)' :(23,16,14,23),'2/m(z)' :( 1,16, 0,-1),'2(yz)' :(10,23,10,22),
2030	'm(yz)' :(22,23,10,22),' 2/m(yz)' :( 1,23, 0,-1),'2(0+-)' :(11,24,11,21),'m(0+-)' :(21,24,11,21),
2031	'2/m(0+-)' :( 1,24, 0,-1),'2(xz)' :( 8,21, 8,20),'m(xz)' :(20,21, 8,20),'2/m(xz)' :( 1,21, 0,-1),
2032	'2(+0-)' :( 9,22, 9,19),'m(+0-)' :(19,22, 9,19),'2/m(+0-)' :( 1,22, 0,-1),'2(xy)' :( 6,19, 6,18),
2033	'm(xy)' :(18,19, 6,18),' 2/m(xy)' :( 1,19, 0,-1),'2(+-0)' :( 7,20, 7,17),'m(+-0)' :(17,20, 7,17),
2034	'2/m(+-0)' :( 1,20, 17,-1),'mm2(x)' :(12,10, 0,-1),'mm2(y)' :(13,10, 0,-1),'mm2(z)' :(14,10, 0,-1),
2035	'mm2(yz)' :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)' :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2036	'mm2(xy)' :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222' :( 1,10, 0,-1),'222(x)' :( 1,13, 0,-1),
2037	'222(y)' :( 1,12, 0,-1),'222(z)' :( 1,11, 0,-1),'mmm' :( 1,10, 0,-1),'mmm(x)' :( 1,13, 0,-1),
2038	'mmm(y)' :( 1,12, 0,-1),'mmm(z)' :( 1,11, 0,-1),'4(x)' :(12, 4,12, 0),'-4(x)' :( 1, 4,12, 0),
2039	'4/m(x)' :( 1, 4,12,-1),'422(x)' :( 1, 4, 0,-1),'-42m(x)' :( 1, 4, 0,-1),'4mm(x)' :(12, 4, 0,-1),
2040	'4/mmm(x)' :( 1, 4, 0,-1),'4(y)' :(13, 3,13, 0),'-4(y)' :( 1, 3,13, 0),'4/m(y)' :( 1, 3,13,-1),
2041	'422(y)' :( 1, 3, 0,-1),'-42m(y)' :( 1, 3, 0,-1),'4mm(y)' :(13, 3, 0,-1),'4/mmm(y)' :(1, 3, 0,-1,),
2042	'4(z)' :(14, 2,14, 0),'-4(z)' :( 1, 2,14, 0),'4/m(z)' :( 1, 2,14,-1),'422(z)' :( 1, 2, 0,-1),
2043	'-42m(z)' :( 1, 2, 0,-1),'4mm(z)' :(14, 2, 0,-1),'4/mmm(z)' :( 1, 2, 0,-1),'3(111)' :( 2, 5, 2, 0),
2044	'-3(111)' :( 1, 5, 2, 0),'32(111)' :( 1, 5, 0, 2),'3m(111)' :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2045	'3(+--)' :( 5, 8, 5, 0),'-3(+--)' :( 1, 8, 5, 0),'32(+--)' :( 1, 8, 0, 5),'3m(+--)' :( 5, 8, 0, 5),
2046	'-3m(+--)' :( 1, 8, 0,-1),'3(-+-)' :( 4, 7, 4, 0),'-3(-+-)' :( 1, 7, 4, 0),'32(-+-)' :( 1, 7, 0, 4),
2047	'3m(-+-)' :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)' :( 3, 6, 3, 0),'-3(--+)' :( 1, 6, 3, 0),
2048	'32(--+)' :( 1, 6, 0, 3),'3m(--+)' :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23' :( 1, 1, 0, 0),
2049	'm3' :( 1, 1, 0, 0),'432' :( 1, 1, 0, 0),'-43m' :( 1, 1, 0, 0),'m3m' :( 1, 1, 0, 0),
2050	'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)' :( 1, 4, 0,-1),
2051	'-4m2(y)' :( 1, 3, 0,-1),'-4m2(z)' :( 1, 2, 0,-1),'6/mmm' :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2052	'-6m2(120)':( 1, 9, 0,-1),'6mm' :(14, 9, 0,-1),'622' :( 1, 9, 0,-1),'6/m' :( 1, 9,14,-1),
2053	'-6' :( 1, 9,14, 0),'6' :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2054	'3m(100)' :(14, 9, 0,14),'3m(120)' :(14, 9, 0,14),'32(100)' :( 1, 9, 0,14),'32(120)' :( 1, 9, 0,14),
2055	'-3' :( 1, 9,14, 0),'3' :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2056	'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2057	'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2058	'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2059	'222(110)' :( 1,11, 0,-1),'2/m(001)' :( 1,16,14,-1),'2/m(100)' :( 1,25,12,-1),'2/m(010)' :( 1,28,13,-1),
2060	'2/m(110)' :( 1,19, 6,-1),'2/m(120)' :( 1,27,15,-1),'2/m(210)' :( 1,26,16,-1),'2/m(+-0)' :( 1,20,17,-1),
2061	'm(001)' :(23,16,14,23),'m(100)' :(26,25,12,26),'m(010)' :(27,28,13,27),'m(110)' :(18,19, 6,18),
2062	'm(120)' :(24,27,15,24),'m(210)' :(25,26,16,25),'m(+-0)' :(17,20, 7,17),'2(001)' :(14,16,14,23),
2063	'2(100)' :(12,25,12,26),'2(010)' :(13,28,13,27),'2(110)' :( 6,19, 6,18),'2(120)' :(15,27,15,24),
2064	'2(210)' :(16,26,16,25),'2(+-0)' :( 7,20, 7,17),'-1' :( 1,29,28, 0)
2065	}
2066
2067	CSxinel = [[], # 0th empty - indices are Fortran style
2068	[[0,0,0],[ 0.0, 0.0, 0.0]], #1 0 0 0
2069	[[1,1,1],[ 1.0, 1.0, 1.0]], #2 X X X
2070	[[1,1,1],[ 1.0, 1.0,-1.0]], #3 X X -X
2071	[[1,1,1],[ 1.0,-1.0, 1.0]], #4 X -X X
2072	[[1,1,1],[ 1.0,-1.0,-1.0]], #5 -X X X
2073	[[1,1,0],[ 1.0, 1.0, 0.0]], #6 X X 0
2074	[[1,1,0],[ 1.0,-1.0, 0.0]], #7 X -X 0
2075	[[1,0,1],[ 1.0, 0.0, 1.0]], #8 X 0 X
2076	[[1,0,1],[ 1.0, 0.0,-1.0]], #9 X 0 -X
2077	[[0,1,1],[ 0.0, 1.0, 1.0]], #10 0 Y Y
2078	[[0,1,1],[ 0.0, 1.0,-1.0]], #11 0 Y -Y
2079	[[1,0,0],[ 1.0, 0.0, 0.0]], #12 X 0 0
2080	[[0,1,0],[ 0.0, 1.0, 0.0]], #13 0 Y 0
2081	[[0,0,1],[ 0.0, 0.0, 1.0]], #14 0 0 Z
2082	[[1,1,0],[ 1.0, 2.0, 0.0]], #15 X 2X 0
2083	[[1,1,0],[ 2.0, 1.0, 0.0]], #16 2X X 0
2084	[[1,1,2],[ 1.0, 1.0, 1.0]], #17 X X Z
2085	[[1,1,2],[ 1.0,-1.0, 1.0]], #18 X -X Z
2086	[[1,2,1],[ 1.0, 1.0, 1.0]], #19 X Y X
2087	[[1,2,1],[ 1.0, 1.0,-1.0]], #20 X Y -X
2088	[[1,2,2],[ 1.0, 1.0, 1.0]], #21 X Y Y
2089	[[1,2,2],[ 1.0, 1.0,-1.0]], #22 X Y -Y
2090	[[1,2,0],[ 1.0, 1.0, 0.0]], #23 X Y 0
2091	[[1,0,2],[ 1.0, 0.0, 1.0]], #24 X 0 Z
2092	[[0,1,2],[ 0.0, 1.0, 1.0]], #25 0 Y Z
2093	[[1,1,2],[ 1.0, 2.0, 1.0]], #26 X 2X Z
2094	[[1,1,2],[ 2.0, 1.0, 1.0]], #27 2X X Z
2095	[[1,2,3],[ 1.0, 1.0, 1.0]], #28 X Y Z
2096	]
2097
2098	CSuinel = [[], # 0th empty - indices are Fortran style
2099	[[1,1,1,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,0,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #1 A A A 0 0 0
2100	[[1,1,2,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,0,1,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #2 A A C 0 0 0
2101	[[1,2,1,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #3 A B A 0 0 0
2102	[[1,2,2,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,0,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #4 A B B 0 0 0
2103	[[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #5 A A A D D D
2104	[[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0,-1.0,-1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #6 A A A D -D -D
2105	[[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0,-1.0, 1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #7 A A A D -D D
2106	[[1,1,1,2,2,2],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,0,0,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #8 A A A D D -D
2107	[[1,1,2,1,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,0,1,0,0,0],[1.0,1.0,1.0,0.5,0.0,0.0]], #9 A A C A/2 0 0
2108	[[1,2,3,0,0,0],[ 1.0, 1.0, 1.0, 0.0, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #10 A B C 0 0 0
2109	[[1,1,2,3,0,0],[ 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],[1,0,1,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #11 A A C D 0 0
2110	[[1,2,1,0,3,0],[ 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],[1,1,0,0,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #12 A B A 0 E 0
2111	[[1,2,2,0,0,3],[ 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],[1,1,0,0,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]], #13 A B B 0 0 F
2112	[[1,2,3,2,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.5,0.0]], #14 A B C B/2 0 0
2113	[[1,2,3,1,0,0],[ 1.0, 1.0, 1.0, 0.5, 0.0, 0.0],[1,1,1,0,0,0],[1.0,1.0,1.0,0.0,0.5,0.0]], #15 A B C A/2 0 0
2114	[[1,2,3,4,0,0],[ 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],[1,1,1,1,0,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #16 A B C D 0 0
2115	[[1,2,3,0,4,0],[ 1.0, 1.0, 1.0, 0.0, 1.0, 0.0],[1,1,1,0,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #17 A B C 0 E 0
2116	[[1,2,3,0,0,4],[ 1.0, 1.0, 1.0, 0.0, 0.0, 1.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]], #18 A B C 0 0 F
2117	[[1,1,2,3,4,4],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,0,1,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #19 A A C D E -E
2118	[[1,1,2,3,4,4],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,0,1,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #20 A A C D E E
2119	[[1,2,1,3,4,3],[ 1.0, 1.0, 1.0, 1.0, 1.0,-1.0],[1,1,0,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #21 A B A D E -D
2120	[[1,2,1,3,4,3],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,0,1,1,0],[1.0,1.0,1.0,0.0,0.0,0.0]], #22 A B A D E D
2121	[[1,2,2,3,3,4],[ 1.0, 1.0, 1.0, 1.0,-1.0, 1.0],[1,1,0,1,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]], #23 A B B D -D F
2122	[[1,2,2,3,3,4],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,0,1,0,1],[1.0,1.0,1.0,0.0,0.0,0.0]], #24 A B B D D F
2123	[[1,2,3,2,4,4],[ 1.0, 1.0, 1.0, 0.5, 1.0, 2.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.5,0.0,0.0]], #25 A B C B/2 F/2 F
2124	[[1,2,3,1,0,4],[ 1.0, 1.0, 1.0, 0.5, 0.0, 1.0],[1,1,1,0,0,1],[1.0,1.0,1.0,0.5,0.0,0.0]], #26 A B C A/2 0 F
2125	[[1,2,3,2,4,0],[ 1.0, 1.0, 1.0, 0.5, 1.0, 0.0],[1,1,1,0,1,0],[1.0,1.0,1.0,0.5,0.0,0.0]], #27 A B C B/2 E 0
2126	[[1,2,3,1,4,4],[ 1.0, 1.0, 1.0, 0.5, 1.0, 0.5],[1,1,1,0,1,0],[1.0,1.0,1.0,0.5,0.0,0.0]], #28 A B C A/2 E E/2
2127	[[1,2,3,4,5,6],[ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],[1,1,1,1,1,1],[1.0,1.0,1.0,0.0,0.0,0.0]], #29 A B C D E F
2128	]
2129
2130	################################################################################
2131	#### Site symmetry routines
2132	################################################################################
2133
2134	def GetOprPtrName(key):
2135	'Needs a doc string'
2136	oprName = OprName[key][0]
2137	return oprName.replace('(','').replace(')','')
2138
2139	def GetOprPtrNumber(key):
2140	'Needs a doc string'
2141	return OprName[key][1]
2142
2143	def GetOprName(key):
2144	'Needs a doc string'
2145	return OprName[key][0]
2146
2147	def GetKNsym(key):
2148	'Needs a doc string'
2149	try:
2150	return KNsym[key].strip()
2151	except KeyError:
2152	return 'sp'
2153
2154	def GetNXUPQsym(siteSym):
2155	'''
2156	The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2157	'''
2158	return NXUPQsym[siteSym]
2159
2160	def GetCSxinel(siteSym):
2161	"returns Xyz terms, multipliers, GUI flags"
2162	indx = GetNXUPQsym(siteSym.strip())
2163	return CSxinel[indx[0]]
2164
2165	def GetCSuinel(siteSym):
2166	"returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2167	indx = GetNXUPQsym(siteSym.strip())
2168	return CSuinel[indx[1]]
2169
2170	def GetCSpqinel(SpnFlp,dupDir):
2171	"returns Mxyz terms, multipliers, GUI flags"
2172	CSI = [[1,2,3],[1.0,1.0,1.0]]
2173	for sopr in dupDir:
2174	# print (sopr,dupDir[sopr])
2175	opr = sopr.replace("'",'')
2176	indx = GetNXUPQsym(opr)
2177	if SpnFlp[dupDir[sopr]] > 0:
2178	csi = CSxinel[indx[2]] #P
2179	else:
2180	csi = CSxinel[indx[3]] #Q
2181	# print(opr,indx,csi,CSI)
2182	if not len(csi):
2183	return [[0,0,0],[0.,0.,0.]]
2184	for kcs in [0,1,2]:
2185	if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2186	jcs = CSI[0][kcs]
2187	for ics in [0,1,2]:
2188	if CSI[0][ics] == jcs:
2189	CSI[0][ics] = 0
2190	CSI[1][ics] = 0.
2191	elif CSI[0][ics] > jcs:
2192	CSI[0][ics] = CSI[0][ics]-1
2193	elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2194	CSI[1][kcs] = csi[1][kcs]
2195	elif CSI[0][kcs] >= csi[0][kcs]:
2196	CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2197	if CSI[1][kcs] != csi[1][kcs]:
2198	if CSI[1][kcs] == 1.:
2199	CSI[1][kcs] = csi[1][kcs]
2200	# print(CSI)
2201	return CSI
2202
2203	def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2204	phase = np.sum(XYZ*wave)
2205	ssopinv = nl.inv(ssop[0])
2206	mst = ssopinv[3][:3]
2207	epsinv = ssopinv[3][3]
2208	sdet = nl.det(sop[0])
2209	ssdet = nl.det(ssop[0])
2210	dtau = mst(XYZ-sop[1])-epsinvssop[1][3]
2211	dT = 1.0
2212	if np.any(dtau%.5):
2213	sumdtau = np.sum(dtau%.5)
2214	dT = 0.
2215	if np.abs(sumdtau-.5) > 1.e-4:
2216	dT = np.tan(np.pi*sumdtau)
2217	tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2218	return sdet,ssdet,dtau,dT,tauT
2219
2220	def OpsfromStringOps(A,SGData,SSGData):
2221	SGOps = SGData['SGOps']
2222	SSGOps = SSGData['SSGOps']
2223	Ax = A.split('+')
2224	Ax[0] = int(Ax[0])
2225	iC = 1
2226	if Ax[0] < 0:
2227	iC = -1
2228	Ax[0] = abs(Ax[0])
2229	nA = Ax[0]%100-1
2230	nC = Ax[0]//100
2231	unit = [0,0,0]
2232	if len(Ax) > 1:
2233	unit = eval('['+Ax[1]+']')
2234	return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2235
2236	def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2237
2238	def orderParms(CSI):
2239	parms = [0,]
2240	for csi in CSI:
2241	for i in [0,1,2]:
2242	if csi[i] not in parms:
2243	parms.append(csi[i])
2244	for csi in CSI:
2245	for i in [0,1,2]:
2246	csi[i] = parms.index(csi[i])
2247	return CSI
2248
2249	def fracCrenel(tau,Toff,Twid):
2250	Tau = (tau-Toff[:,np.newaxis])%1.
2251	A = np.where(Tau<Twid[:,np.newaxis],1.,0.)
2252	return A
2253
2254	def fracFourier(tau,nH,fsin,fcos):
2255	SA = np.sin(2.nHnp.pi*tau)
2256	CB = np.cos(2.nHnp.pi*tau)
2257	A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
2258	B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
2259	return A+B
2260
2261	def posFourier(tau,nH,psin,pcos):
2262	SA = np.sin(2nHnp.pi*tau)
2263	CB = np.cos(2nHnp.pi*tau)
2264	A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
2265	B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
2266	return A+B
2267
2268	def posSawtooth(tau,Toff,slopes):
2269	Tau = (tau-Toff)%1.
2270	A = slopes[:,np.newaxis]*Tau
2271	return A
2272
2273	def posZigZag(tau,Tmm,XYZmax):
2274	DT = Tmm[1]-Tmm[0]
2275	slopeUp = 2.*XYZmax/DT
2276	slopeDn = 2.*XYZmax/(1.-DT)
2277	A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp((t-Tmm[0])%1.),XYZmax-slopeDn((t-Tmm[1])%1.)) for t in tau])
2278	return A
2279
2280	def posBlock(tau,Tmm,XYZmax):
2281	A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2282	return A
2283
2284	def DoFrac():
2285	dF = None
2286	dFTP = None
2287	if siteSym == '1':
2288	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2289	elif siteSym == '-1':
2290	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2291	else:
2292	delt2 = np.eye(2)*0.001
2293	FSC = np.ones(2,dtype='i')
2294	CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2295	if 'Crenel' in waveType:
2296	dF = np.zeros_like(tau)
2297	else:
2298	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2299	dFT = np.zeros_like(dF)
2300	dFTP = []
2301	for i in SdIndx:
2302	sop = Sop[i]
2303	ssop = SSop[i]
2304	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2305	fsc = np.ones(2,dtype='i')
2306	if 'Crenel' in waveType:
2307	dFT = np.zeros_like(tau)
2308	fsc = [1,1]
2309	else: #Fourier
2310	dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2311	dFT = nl.det(sop[0])*dFT
2312	dFT = dFT[:,np.argsort(tauT)]
2313	dFT[0] *= ssdet
2314	dFT[1] *= sdet
2315	dFTP.append(dFT)
2316
2317	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2318	fsc = [1,1]
2319	if dT:
2320	CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2321	else:
2322	CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2323	FSC = np.zeros(2,dtype='i')
2324	return CSI,dF,dFTP
2325	else:
2326	for i in range(2):
2327	if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2328	fsc[i] = 1
2329	else:
2330	fsc[i] = 0
2331	FSC &= fsc
2332	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2333	n = -1
2334	for i,F in enumerate(FSC):
2335	if F:
2336	n += 1
2337	CSI[0][i] = n+1
2338	CSI[1][i] = 1.0
2339
2340	return CSI,dF,dFTP
2341
2342	def DoXYZ():
2343	dX,dXTP = None,None
2344	if siteSym == '1':
2345	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2346	elif siteSym == '-1':
2347	CSI = [[1,0,0],[2,0,0],[3,0,0], [0,0,0],[0,0,0],[0,0,0]],3[[1.,0.,0.],]+3[[0.,0.,0.],]
2348	else:
2349	delt4 = np.ones(4)*0.001
2350	delt5 = np.ones(5)*0.001
2351	delt6 = np.eye(6)*0.001
2352	if 'Fourier' in waveType:
2353	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2354	#3x6x12 modulated position array (X,Spos,tau)& force positive
2355	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2356	elif waveType == 'Sawtooth':
2357	dX = posSawtooth(tau,delt4[0],delt4[1:])
2358	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2359	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2360	elif waveType in ['ZigZag','Block']:
2361	if waveType == 'ZigZag':
2362	dX = posZigZag(tau,delt5[:2],delt5[2:])
2363	else:
2364	dX = posBlock(tau,delt5[:2],delt5[2:])
2365	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2366	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2367	XSC = np.ones(6,dtype='i')
2368	dXTP = []
2369	for i in SdIndx:
2370	sop = Sop[i]
2371	ssop = SSop[i]
2372	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2373	xsc = np.ones(6,dtype='i')
2374	if 'Fourier' in waveType:
2375	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2376	elif waveType == 'Sawtooth':
2377	dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2378	elif waveType == 'ZigZag':
2379	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2380	elif waveType == 'Block':
2381	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2382	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2383	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2384	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2385	dXTP.append(dXT)
2386	if waveType == 'Fourier':
2387	for i in range(3):
2388	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2389	xsc[i] = 0
2390	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2391	xsc[i+3] = 0
2392	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2393	xsc[3:6] = 0
2394	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2395	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2396	if dT:
2397	if '(x)' in siteSym:
2398	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2399	if 'm' in siteSym and len(SdIndx) == 1:
2400	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2401	elif '(y)' in siteSym:
2402	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2403	if 'm' in siteSym and len(SdIndx) == 1:
2404	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2405	elif '(z)' in siteSym:
2406	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2407	if 'm' in siteSym and len(SdIndx) == 1:
2408	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2409	else:
2410	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2411	if '4/mmm' in laue:
2412	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2413	if '(xy)' in siteSym:
2414	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2415	if dT:
2416	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2417	else:
2418	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2419	if '(xy)' in siteSym or '(+-0)' in siteSym:
2420	mul = 1
2421	if '(+-0)' in siteSym:
2422	mul = -1
2423	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2424	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2425	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2426	xsc[3:5] = 0
2427	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2428	CSI[0][:2] = [[12,0,0],[12,0,0]]
2429	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2430	xsc[:2] = 0
2431	XSC &= xsc
2432	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2433	if waveType == 'Fourier':
2434	n = -1
2435	if debug: print (XSC)
2436	for i,X in enumerate(XSC):
2437	if X:
2438	n += 1
2439	CSI[0][i][0] = n+1
2440	CSI[1][i][0] = 1.0
2441
2442	return list(CSI),dX,dXTP
2443
2444	def DoUij():
2445	dU,dUTP = None,None
2446	if siteSym == '1':
2447	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2448	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2449	elif siteSym == '-1':
2450	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2451	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2452	else:
2453	tau = np.linspace(0,1,49,True)
2454	delt12 = np.eye(12)*0.0001
2455	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2456	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2457	USC = np.ones(12,dtype='i')
2458	dUTP = []
2459	dtau = 0.
2460	for i in SdIndx:
2461	sop = Sop[i]
2462	ssop = SSop[i]
2463	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2464	usc = np.ones(12,dtype='i')
2465	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2466	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2467	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2468	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2469	dUT = dUT[:,:,np.argsort(tauT)]
2470	dUT[:,:6,:] =(ssdetsdet)
2471	dUTP.append(dUT)
2472	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2473	if dT:
2474	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2475	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2476	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2477	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2478	else:
2479	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2480	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2481	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2482	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2483	if 'mm2(x)' in siteSym and dT:
2484	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2485	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2486	elif '(xy)' in siteSym and dT:
2487	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2488	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2489	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2490	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2491	elif '(x)' in siteSym and dT:
2492	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2493	elif '(y)' in siteSym and dT:
2494	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2495	elif '(z)' in siteSym and dT:
2496	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2497	for i in range(6):
2498	if not USC[i]:
2499	CSI[0][i] = [0,0,0]
2500	CSI[1][i] = [0.,0.,0.]
2501	CSI[0][i+6] = [0,0,0]
2502	CSI[1][i+6] = [0.,0.,0.]
2503	else:
2504	for i in range(6):
2505	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2506	usc[i] = 0
2507	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2508	usc[i+6] = 0
2509	if np.any(dUT[1,0,:]):
2510	if '4/m' in siteSym:
2511	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2512	if ssop[1][3]:
2513	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2514	usc[9] = 1
2515	else:
2516	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2517	usc[9] = 0
2518	elif '4' in siteSym:
2519	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2520	CSI[0][:2] = [[11,0,0],[11,0,0]]
2521	if ssop[1][3]:
2522	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2523	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2524	usc[2] = 0
2525	usc[8] = 0
2526	usc[3] = 1
2527	usc[9] = 1
2528	else:
2529	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2530	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2531	usc[2] = 1
2532	usc[8] = 1
2533	usc[3] = 0
2534	usc[9] = 0
2535	elif 'xy' in siteSym or '+-0' in siteSym:
2536	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2537	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2538	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2539	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2540	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2541	usc[4:6] = 0
2542	usc[6:8] = 0
2543
2544	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2545	USC &= usc
2546	if debug: print (USC)
2547	if not np.any(dtau%.5):
2548	n = -1
2549	for i,U in enumerate(USC):
2550	if U:
2551	n += 1
2552	CSI[0][i][0] = n+1
2553	CSI[1][i][0] = 1.0
2554
2555	return list(CSI),dU,dUTP
2556
2557	def DoMag():
2558	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2559	if siteSym == '1':
2560	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2561	elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2562	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2563	else:
2564	delt6 = np.eye(6)*0.001
2565	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,np.newaxis,np.newaxis]
2566	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2567	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2568	MSC = np.ones(6,dtype='i')
2569	dMTP = []
2570	for i in SdIndx:
2571	sop = Sop[i]
2572	ssop = SSop[i]
2573	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2574	msc = np.ones(6,dtype='i')
2575	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2576	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2577	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2578	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2579	dMTP.append(dMT)
2580	for i in range(3):
2581	if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2582	msc[i] = 0
2583	if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2584	msc[i+3] = 0
2585	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2586	msc[3:6] = 0
2587	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2588	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2589	if dT:
2590	if '(x)' in siteSym:
2591	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2592	if 'm' in siteSym and len(SdIndx) == 1:
2593	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2594	elif '(y)' in siteSym:
2595	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2596	if 'm' in siteSym and len(SdIndx) == 1:
2597	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2598	elif '(z)' in siteSym:
2599	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2600	if 'm' in siteSym and len(SdIndx) == 1:
2601	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2602	else:
2603	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2604	if '4/mmm' in laue:
2605	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2606	if '(xy)' in siteSym:
2607	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2608	if dT:
2609	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2610	else:
2611	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2612	if '(xy)' in siteSym or '(+-0)' in siteSym:
2613	mul = 1
2614	if '(+-0)' in siteSym:
2615	mul = -1
2616	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2617	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2618	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2619	msc[3:5] = 0
2620	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2621	CSI[0][:2] = [[12,0,0],[12,0,0]]
2622	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2623	msc[:2] = 0
2624	MSC &= msc
2625	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2626	n = -1
2627	if debug: print (MSC)
2628	for i,M in enumerate(MSC):
2629	if M:
2630	n += 1
2631	CSI[0][i][0] = n+1
2632	CSI[1][i][0] = 1.0
2633
2634	return CSI,None,None
2635
2636	if debug: print ('super space group: '+SSGData['SSpGrp'])
2637	xyz = np.array(XYZ)%1.
2638	SGOps = copy.deepcopy(SGData['SGOps'])
2639	laue = SGData['SGLaue']
2640	siteSym = SytSym(XYZ,SGData)[0].strip()
2641	if debug: print ('siteSym: '+siteSym)
2642	SSGOps = copy.deepcopy(SSGData['SSGOps'])
2643	#expand ops to include inversions if any
2644	if SGData['SGInv']:
2645	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2646	SGOps.append([-op[0],-op[1]%1.])
2647	SSGOps.append([-sop[0],-sop[1]%1.])
2648	#build set of sym ops around special position
2649	SSop = []
2650	Sop = []
2651	Sdtau = []
2652	for iop,Op in enumerate(SGOps):
2653	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2654	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2655	SSop.append(SSGOps[iop])
2656	Sop.append(SGOps[iop])
2657	ssopinv = nl.inv(SSGOps[iop][0])
2658	mst = ssopinv[3][:3]
2659	epsinv = ssopinv[3][3]
2660	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
2661	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
2662	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2663	#setup displacement arrays
2664	tau = np.linspace(-1,1,49,True)
2665	#make modulation arrays - one parameter at a time
2666	if Stype == 'Sfrac':
2667	CSI,dF,dFTP = DoFrac()
2668	elif Stype == 'Spos':
2669	CSI,dF,dFTP = DoXYZ()
2670	CSI[0] = orderParms(CSI[0])
2671	elif Stype == 'Sadp':
2672	CSI,dF,dFTP = DoUij()
2673	CSI[0] = orderParms(CSI[0])
2674	elif Stype == 'Smag':
2675	CSI,dF,dFTP = DoMag()
2676
2677	if debug:
2678	return CSI,dF,dFTP
2679	else:
2680	return CSI
2681
2682	def MustrainNames(SGData):
2683	'Needs a doc string'
2684	laue = SGData['SGLaue']
2685	uniq = SGData['SGUniq']
2686	if laue in ['m3','m3m']:
2687	return ['S400','S220']
2688	elif laue in ['6/m','6/mmm','3m1']:
2689	return ['S400','S004','S202']
2690	elif laue in ['31m','3']:
2691	return ['S400','S004','S202','S211']
2692	elif laue in ['3R','3mR']:
2693	return ['S400','S220','S310','S211']
2694	elif laue in ['4/m','4/mmm']:
2695	return ['S400','S004','S220','S022']
2696	elif laue in ['mmm']:
2697	return ['S400','S040','S004','S220','S202','S022']
2698	elif laue in ['2/m']:
2699	SHKL = ['S400','S040','S004','S220','S202','S022']
2700	if uniq == 'a':
2701	SHKL += ['S013','S031','S211']
2702	elif uniq == 'b':
2703	SHKL += ['S301','S103','S121']
2704	elif uniq == 'c':
2705	SHKL += ['S130','S310','S112']
2706	return SHKL
2707	else:
2708	SHKL = ['S400','S040','S004','S220','S202','S022']
2709	SHKL += ['S310','S103','S031','S130','S301','S013']
2710	SHKL += ['S211','S121','S112']
2711	return SHKL
2712
2713	def HStrainVals(HSvals,SGData):
2714	laue = SGData['SGLaue']
2715	uniq = SGData['SGUniq']
2716	DIJ = np.zeros(6)
2717	if laue in ['m3','m3m']:
2718	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2719	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2720	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2721	elif laue in ['3R','3mR']:
2722	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2723	elif laue in ['4/m','4/mmm']:
2724	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2725	elif laue in ['mmm']:
2726	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2727	elif laue in ['2/m']:
2728	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2729	if uniq == 'a':
2730	DIJ[5] = HSvals[3]
2731	elif uniq == 'b':
2732	DIJ[4] = HSvals[3]
2733	elif uniq == 'c':
2734	DIJ[3] = HSvals[3]
2735	else:
2736	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
2737	return DIJ
2738
2739	def HStrainNames(SGData):
2740	'Needs a doc string'
2741	laue = SGData['SGLaue']
2742	uniq = SGData['SGUniq']
2743	if laue in ['m3','m3m']:
2744	return ['D11','eA'] #add cubic strain term
2745	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2746	return ['D11','D33']
2747	elif laue in ['3R','3mR']:
2748	return ['D11','D12']
2749	elif laue in ['4/m','4/mmm']:
2750	return ['D11','D33']
2751	elif laue in ['mmm']:
2752	return ['D11','D22','D33']
2753	elif laue in ['2/m']:
2754	Dij = ['D11','D22','D33']
2755	if uniq == 'a':
2756	Dij += ['D23']
2757	elif uniq == 'b':
2758	Dij += ['D13']
2759	elif uniq == 'c':
2760	Dij += ['D12']
2761	return Dij
2762	else:
2763	Dij = ['D11','D22','D33','D12','D13','D23']
2764	return Dij
2765
2766	def MustrainCoeff(HKL,SGData):
2767	'Needs a doc string'
2768	#NB: order of terms is the same as returned by MustrainNames
2769	laue = SGData['SGLaue']
2770	uniq = SGData['SGUniq']
2771	h,k,l = HKL
2772	Strm = []
2773	if laue in ['m3','m3m']:
2774	Strm.append(h4+k4+l**4)
2775	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
2776	elif laue in ['6/m','6/mmm','3m1']:
2777	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
2778	Strm.append(l**4)
2779	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
2780	elif laue in ['31m','3']:
2781	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
2782	Strm.append(l**4)
2783	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
2784	Strm.append(4.0hkl(h+k))
2785	elif laue in ['3R','3mR']:
2786	Strm.append(h4+k4+l**4)
2787	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
2788	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
2789	Strm.append(4.0(klh2+hlk2+hkl*2))
2790	elif laue in ['4/m','4/mmm']:
2791	Strm.append(h4+k4)
2792	Strm.append(l**4)
2793	Strm.append(3.0(hk)**2)
2794	Strm.append(3.0((hl)*2+(kl)**2))
2795	elif laue in ['mmm']:
2796	Strm.append(h**4)
2797	Strm.append(k**4)
2798	Strm.append(l**4)
2799	Strm.append(3.0(hk)**2)
2800	Strm.append(3.0(hl)**2)
2801	Strm.append(3.0(kl)**2)
2802	elif laue in ['2/m']:
2803	Strm.append(h**4)
2804	Strm.append(k**4)
2805	Strm.append(l**4)
2806	Strm.append(3.0(hk)**2)
2807	Strm.append(3.0(hl)**2)
2808	Strm.append(3.0(kl)**2)
2809	if uniq == 'a':
2810	Strm.append(2.0kl**3)
2811	Strm.append(2.0lk**3)
2812	Strm.append(4.0klh*2)
2813	elif uniq == 'b':
2814	Strm.append(2.0lh**3)
2815	Strm.append(2.0hl**3)
2816	Strm.append(4.0hlk*2)
2817	elif uniq == 'c':
2818	Strm.append(2.0hk**3)
2819	Strm.append(2.0kh**3)
2820	Strm.append(4.0hkl*2)
2821	else:
2822	Strm.append(h**4)
2823	Strm.append(k**4)
2824	Strm.append(l**4)
2825	Strm.append(3.0(hk)**2)
2826	Strm.append(3.0(hl)**2)
2827	Strm.append(3.0(kl)**2)
2828	Strm.append(2.0kh**3)
2829	Strm.append(2.0hl**3)
2830	Strm.append(2.0lk**3)
2831	Strm.append(2.0hk**3)
2832	Strm.append(2.0lh**3)
2833	Strm.append(2.0kl**3)
2834	Strm.append(4.0klh*2)
2835	Strm.append(4.0hlk*2)
2836	Strm.append(4.0khl*2)
2837	return Strm
2838
2839	def MuShklMean(SGData,Amat,Shkl):
2840
2841	def genMustrain(xyz,Shkl):
2842	uvw = np.inner(Amat.T,xyz)
2843	Strm = np.array(MustrainCoeff(uvw,SGData))
2844	Sum = np.sum(np.multiply(Shkl,Strm))
2845	Sum = np.where(Sum > 0.01,Sum,0.01)
2846	Sum = np.sqrt(Sum)
2847	return Sum*xyz
2848
2849	PHI = np.linspace(0.,360.,30,True)
2850	PSI = np.linspace(0.,180.,30,True)
2851	X = np.outer(npcosd(PHI),npsind(PSI))
2852	Y = np.outer(npsind(PHI),npsind(PSI))
2853	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2854	XYZ = np.dstack((X,Y,Z))
2855	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
2856	return np.sqrt(np.sum(XYZ**2)/900.)
2857
2858	def Muiso2Shkl(muiso,SGData,cell):
2859	"this is to convert isotropic mustrain to generalized Shkls"
2860	import GSASIIlattice as G2lat
2861	A = G2lat.cell2AB(cell)[0]
2862
2863	def minMus(Shkl,muiso,H,SGData,A):
2864	U = np.inner(A.T,H)
2865	S = np.array(MustrainCoeff(U,SGData))
2866	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
2867	rad = np.sqrt(np.sum((Sum[:,np.newaxis]H)*2,axis=1))
2868	return (muiso-rad)**2
2869
2870	laue = SGData['SGLaue']
2871	PHI = np.linspace(0.,360.,60,True)
2872	PSI = np.linspace(0.,180.,60,True)
2873	X = np.outer(npsind(PHI),npsind(PSI))
2874	Y = np.outer(npcosd(PHI),npsind(PSI))
2875	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2876	HKL = np.dstack((X,Y,Z))
2877	if laue in ['m3','m3m']:
2878	S0 = [1000.,1000.]
2879	elif laue in ['6/m','6/mmm','3m1']:
2880	S0 = [1000.,1000.,1000.]
2881	elif laue in ['31m','3']:
2882	S0 = [1000.,1000.,1000.,1000.]
2883	elif laue in ['3R','3mR']:
2884	S0 = [1000.,1000.,1000.,1000.]
2885	elif laue in ['4/m','4/mmm']:
2886	S0 = [1000.,1000.,1000.,1000.]
2887	elif laue in ['mmm']:
2888	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
2889	elif laue in ['2/m']:
2890	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
2891	else:
2892	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
2893	1000.,1000.,0.,0.,0.]
2894	S0 = np.array(S0)
2895	HKL = np.reshape(HKL,(-1,3))
2896	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
2897	return result[0]
2898
2899	def PackRot(SGOps):
2900	IRT = []
2901	for ops in SGOps:
2902	M = ops[0]
2903	irt = 0
2904	for j in range(2,-1,-1):
2905	for k in range(2,-1,-1):
2906	irt *= 3
2907	irt += M[k][j]
2908	IRT.append(int(irt))
2909	return IRT
2910
2911	def SytSym(XYZ,SGData):
2912	'''
2913	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
2914
2915	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
2916	:param SGData: from SpcGroup
2917	:Returns: a four element tuple:
2918
2919	* The 1st element is a code for the site symmetry (see GetKNsym)
2920	* The 2nd element is the site multiplicity
2921	* Ndup number of overlapping operators
2922	* dupDir Dict - dictionary of overlapping operators
2923
2924	'''
2925	Mult = 1
2926	Isym = 0
2927	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2928	Isym = 1073741824
2929	Jdup = 0
2930	Ndup = 0
2931	dupDir = {}
2932	inv = SGData['SGInv']+1
2933	icen = SGData['SGCen']
2934	if SGData['SGFixed']: #already in list of operators
2935	inv = 1
2936	Xeqv = list(GenAtom(XYZ,SGData,True))
2937	# for xeqv in Xeqv: print(xeqv)
2938	IRT = PackRot(SGData['SGOps'])
2939	L = -1
2940	for ic,cen in enumerate(icen):
2941	for invers in range(int(inv)):
2942	for io,ops in enumerate(SGData['SGOps']):
2943	irtx = (1-2invers)IRT[io]
2944	L += 1
2945	if not Xeqv[L][1]:
2946	Ndup = io
2947	Jdup += 1
2948	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
2949	if jx < 39:
2950	px = GetOprName(str(irtx))
2951	if Xeqv[L][-1] < 0:
2952	if '(' in px:
2953	px = px.split('(')
2954	px[0] += "'"
2955	px = '('.join(px)
2956	else:
2957	px += "'"
2958	dupDir[px] = L
2959	Isym += 2**(jx-1)
2960	if Isym == 1073741824: Isym = 0
2961	Mult = len(SGData['SGOps'])len(SGData['SGCen'])inv//Jdup
2962
2963	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
2964
2965	def MagSytSym(SytSym,dupDir,SGData):
2966	'''
2967	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
2968	'''
2969	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
2970	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
2971	# print('dupDir',dupDir)
2972	SplitSytSym = SytSym.split('(')
2973	if SGData['SGGray']:
2974	return SytSym+"1'"
2975	if SytSym == '1': #genersl position
2976	return SytSym
2977	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
2978	try:
2979	MagSytSym = SGData['MagSpGrp'].split()[1]
2980	except IndexError:
2981	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
2982	if len(SplitSytSym) > 1:
2983	MagSytSym += '('+SplitSytSym[1]
2984	return MagSytSym
2985	if len(dupDir) == 1:
2986	return dupDir.keys()[0]
2987
2988
2989	if '2/m' in SytSym: #done I think; last 2wo might be not needed
2990	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
2991	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
2992	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
2993	'(110)':['2(110)','m(110)']}
2994
2995	elif '4/mmm' in SytSym:
2996	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
2997	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
2998	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
2999	elif '4mm' in SytSym:
3000	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3001	elif '422' in SytSym:
3002	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3003	elif '-4m2' in SytSym:
3004	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3005	elif '-42m' in SytSym:
3006	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3007	elif '-4' in SytSym:
3008	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3009	elif '4' in SytSym:
3010	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3011
3012	elif '222' in SytSym:
3013	ops = {'':['2(x)','2(y)','2(z)'],
3014	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3015	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3016	'(110)':['2(z)','2(110)','2(+-0)',],}
3017	elif 'mm2' in SytSym:
3018	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3019	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3020	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3021	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3022	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3023	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3024	'(210)':['m(z)','m(010)','2(210)'],
3025	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3026	'(110)':['m(z)','m(+-0)','2(110)',],}
3027	elif 'mmm' in SytSym:
3028	ops = {'':['m(x)','m(y)','m(z)'],
3029	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3030	'(110)':['m(z)','m(110)','m(+-0)',],
3031	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3032
3033	elif '32' in SytSym:
3034	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3035	elif '23' in SytSym:
3036	ops = {'':['2(x)','3(111)']}
3037	elif 'm3' in SytSym:
3038	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3039	elif '3m' in SytSym:
3040	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3041	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3042	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3043
3044	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]: #not simple cases
3045	MagSytSym = SytSym
3046	if "-1'" in dupDir:
3047	if '-6' in SytSym:
3048	MagSytSym = MagSytSym.replace('-6',"-6'")
3049	elif '-3m' in SytSym:
3050	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3051	elif '-3' in SytSym:
3052	MagSytSym = MagSytSym.replace('-3',"-3'")
3053	elif '-6m2' in SytSym:
3054	if "m'(110)" in dupDir:
3055	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3056	elif '6/m' in SytSym:
3057	if "m'(z)" in dupDir:
3058	MagSytSym = "6'/m'"
3059	elif '6mm' in SytSym:
3060	if "m'(110)" in dupDir:
3061	MagSytSym = "6'm'm"
3062	return MagSytSym
3063	try:
3064	axis = '('+SytSym.split('(')[1]
3065	except IndexError:
3066	axis = ''
3067	MagSytSym = ''
3068	for m in ops[axis]:
3069	if m in dupDir:
3070	MagSytSym += m.split('(')[0]
3071	else:
3072	MagSytSym += m.split('(')[0]+"'"
3073	if '2/m' in SytSym and '2' in m:
3074	MagSytSym += '/'
3075	if '-3/m' in SytSym:
3076	MagSytSym = '-'+MagSytSym
3077
3078	MagSytSym += axis
3079	# some exceptions & special rules
3080	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3081	return MagSytSym
3082
3083	# if len(GenSym) == 3:
3084	# if SGSpin[1] < 0:
3085	# if 'mm2' in SytSym:
3086	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3087	# else: #bad rule for I41/a
3088	# MagSytSym = SplitSytSym[0]+"'"
3089	# if len(SplitSytSym) > 1:
3090	# MagSytSym += '('+SplitSytSym[1]
3091	# else:
3092	# MagSytSym = SytSym
3093	# if len(SplitSytSym) >1:
3094	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3095	# MagSytSym = MagSytSym.replace('-4',"-4'")
3096	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3097	# MagSytSym = MagSytSym.replace('-6',"-6'")
3098	# return MagSytSym
3099	#
3100	return SytSym
3101
3102	def ElemPosition(SGData):
3103	''' Under development.
3104	Object here is to return a list of symmetry element types and locations suitable
3105	for say drawing them.
3106	So far I have the element type... getting all possible locations without lookup may be impossible!
3107	'''
3108	Inv = SGData['SGInv']
3109	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3110	# get operators & expand if centrosymmetric
3111	Ops = SGData['SGOps']
3112	opM = np.array([op[0].T for op in Ops])
3113	opT = np.array([op[1] for op in Ops])
3114	if Inv:
3115	opM = np.concatenate((opM,-opM))
3116	opT = np.concatenate((opT,-opT))
3117	opMT = list(zip(opM,opT))
3118	for M,T in opMT[1:]: #skip I
3119	Dt = int(nl.det(M))
3120	Tr = int(np.trace(M))
3121	Dt = -(Dt-1)//2
3122	Es = eleSym[Tr][Dt]
3123	if Dt: #rotation-inversion
3124	I = np.eye(3)
3125	if Tr == 1: #mirrors/glides
3126	if np.any(T): #glide
3127	M2 = np.inner(M,M)
3128	MT = np.inner(M,T)+T
3129	print ('glide',Es,MT)
3130	print (M2)
3131	else: #mirror
3132	print ('mirror',Es,T)
3133	print (I-M)
3134	X = [-1,-1,-1]
3135	elif Tr == -3: # pure inversion
3136	X = np.inner(nl.inv(I-M),T)
3137	print ('inversion',Es,X)
3138	else: #other rotation-inversion
3139	M2 = np.inner(M,M)
3140	MT = np.inner(M,T)+T
3141	print ('rot-inv',Es,MT)
3142	print (M2)
3143	X = [-1,-1,-1]
3144	else: #rotations
3145	print ('rotation',Es)
3146	X = [-1,-1,-1]
3147	#SymElements.append([Es,X])
3148
3149	return #SymElements
3150
3151	def ApplyStringOps(A,SGData,X,Uij=[]):
3152	'Needs a doc string'
3153	SGOps = SGData['SGOps']
3154	SGCen = SGData['SGCen']
3155	Ax = A.split('+')
3156	Ax[0] = int(Ax[0])
3157	iC = 0
3158	if Ax[0] < 0:
3159	iC = 1
3160	Ax[0] = abs(Ax[0])
3161	nA = Ax[0]%100-1
3162	cA = Ax[0]//100
3163	Cen = SGCen[cA]
3164	M,T = SGOps[nA]
3165	if len(Ax)>1:
3166	cellA = Ax[1].split(',')
3167	cellA = np.array([int(a) for a in cellA])
3168	else:
3169	cellA = np.zeros(3)
3170	newX = Cen+(1-2iC)(np.inner(M,X).T+T)+cellA
3171	if len(Uij):
3172	U = Uij2U(Uij)
3173	U = np.inner(M,np.inner(U,M).T)
3174	newUij = U2Uij(U)
3175	return [newX,newUij]
3176	else:
3177	return newX
3178
3179	def ApplyStringOpsMom(A,SGData,Mom):
3180	'Needs a doc string'
3181	SGOps = SGData['SGOps']
3182	Ax = A.split('+')
3183	Ax[0] = int(Ax[0])
3184	Ax[0] = abs(Ax[0])
3185	nA = Ax[0]%100-1
3186	M,T = SGOps[nA]
3187	if SGData['SGGray']:
3188	newMom = -(np.inner(Mom,M).T)*nl.det(M)
3189	else:
3190	newMom = -(np.inner(Mom,M).T)SGData['MagMom'][nA-1]nl.det(M)
3191	return newMom
3192
3193	def StringOpsProd(A,B,SGData):
3194	"""
3195	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3196	where '-' indicates inversion, c(>0) is the cell centering operator,
3197	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3198	Should return resultant string - C. SGData - dictionary using entries:
3199
3200	* 'SGCen': cell centering vectors [0,0,0] at least
3201	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3202
3203	"""
3204	SGOps = SGData['SGOps']
3205	SGCen = SGData['SGCen']
3206	#1st split out the cell translation part & work on the operator parts
3207	Ax = A.split('+'); Bx = B.split('+')
3208	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3209	iC = 0
3210	if Ax[0]*Bx[0] < 0:
3211	iC = 1
3212	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3213	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3214	cA = Ax[0]//100; cB = Bx[0]//100
3215	Cen = (SGCen[cA]+SGCen[cB])%1.0
3216	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3217	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3218	Mc = np.inner(Ma,Mb.T)
3219	# print Ma,Mb,Mc
3220	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3221	# print Ta,Tb,Tc
3222	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3223	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3224	#now the cell translation part
3225	if len(Ax)>1:
3226	cellA = Ax[1].split(',')
3227	cellA = [int(a) for a in cellA]
3228	else:
3229	cellA = [0,0,0]
3230	if len(Bx)>1:
3231	cellB = Bx[1].split(',')
3232	cellB = [int(b) for b in cellB]
3233	else:
3234	cellB = [0,0,0]
3235	cellC = np.add(cellA,cellB)
3236	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3237	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3238	return C
3239
3240	def U2Uij(U):
3241	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3242	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3243
3244	def Uij2U(Uij):
3245	#returns the thermal motion tensor U from Uij as numpy array
3246	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3247
3248	def StandardizeSpcName(spcgroup):
3249	'''Accept a spacegroup name where spaces may have not been used
3250	in the names according to the GSAS convention (spaces between symmetry
3251	for each axis) and return the space group name as used in GSAS
3252	'''
3253	rspc = spcgroup.replace(' ','').upper()
3254	# deal with rhombohedral and hexagonal setting designations
3255	rhomb = ''
3256	if rspc[-1:] == 'R':
3257	rspc = rspc[:-1]
3258	rhomb = ' R'
3259	gray = ''
3260	if "1'" in rspc:
3261	gray = " 1'"
3262	rspc = rspc.replace("1'",'')
3263	elif rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3264	rspc = rspc[:-1]
3265	else:
3266	rspc = rspc.replace("'",'')
3267	# look for a match in the spacegroup lists
3268	for i in spglist.values():
3269	for spc in i:
3270	if rspc == spc.replace(' ','').upper():
3271	return spc+gray+rhomb
3272	# how about the post-2002 orthorhombic names?
3273	if rspc in sgequiv_2002_orthorhombic:
3274	return sgequiv_2002_orthorhombic[rspc]+gray
3275	else:
3276	# not found
3277	return ''
3278
3279	spgbyNum = []
3280	'''Space groups indexed by number'''
3281	spgbyNum = [None,
3282	'P 1','P -1', #1-2
3283	'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3284	'C 2/m','P 2/c','P 21/c','C 2/c', #3-15
3285	'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3286	'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3287	'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3288	'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3289	'C m m 2','C m c 21','C c c 2',
3290	'A m m 2','A b m 2','A m a 2','A b a 2',
3291	'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3292	'P m m m','P n n n','P c c m','P b a n',
3293	'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3294	'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3295	'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3296	'F m m m', 'F d d d',
3297	'I m m m','I b a m','I b c a','I m m a', #16-74
3298	'P 4','P 41','P 42','P 43',
3299	'I 4','I 41',
3300	'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3301	'I 4/m','I 41/a',
3302	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3303	'P 42 21 2','P 43 2 2','P 43 21 2',
3304	'I 4 2 2','I 41 2 2',
3305	'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3306	'P 42 m c','P 42 b c',
3307	'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3308	'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3309	'P -4 c 2','P -4 b 2','P -4 n 2',
3310	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3311	'P 4/m m m','P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m',
3312	'P 4/m n c','P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
3313	'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m','P 42/n m c',
3314	'P 42/n c m',
3315	'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3316	'P 3','P 31','P 32','R 3','P -3','R -3',
3317	'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3318	'R 3 2',
3319	'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3320	'R 3 m','R 3 c',
3321	'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3322	'R -3 m','R -3 c', #75-167
3323	'P 6','P 61',
3324	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3325	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3326	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3327	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c', #168-194
3328	'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3329	'F m -3','F d -3','I m -3',
3330	'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3331	'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3332	'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3333	'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3334	'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3335	'I m -3 m','I a -3 d',] #195-230
3336	spg2origins = {}
3337	''' A dictionary of all spacegroups that have 2nd settings; the value is the
3338	1st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3339	'''
3340	spg2origins = {
3341	'P n n n':[-.25,-.25,-.25],
3342	'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3343	'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3344	'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3345	'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3346	'F d d d':[-.125,-.125,-.125],
3347	'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3348	'P 4/n b m':[-.25,-.25,0],'P 4/n n c':[-.25,-.25,-.25],'P 4/n m m':[-.25,-.25,0],'P 4/n c c':[-.25,-.25,0],
3349	'P 42/n b c':[-.25,-.25,-.25],'P 42/n n m':[-.25,.25,-.25],'P 42/n m c':[-.25,.25,-.25],'P 42/n c m':[-.25,.25,-.25],
3350	'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3351	'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3352	'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3353	'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3354	'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3355	spglist = {}
3356	'''A dictionary of space groups as ordered and named in the pre-2002 International
3357	Tables Volume A, except that spaces are used following the GSAS convention to
3358	separate the different crystallographic directions.
3359	Note that the symmetry codes here will recognize many non-standard space group
3360	symbols with different settings. They are ordered by Laue group
3361	'''
3362	spglist = {
3363	'P1' : ('P 1','P -1',), # 1-2
3364	'C1' : ('C 1','C -1',),
3365	'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3366	'P 2/m','P 21/m','P 2/c','P 2/a','P 2/n','P 21/c','P 21/a','P 21/n',), #3-15
3367	'C2/m':('C 2','C m','C c','C n',
3368	'C 2/m','C 2/c','C 2/n',),
3369	'A2/m':('A 2','A m','A a','A n',
3370	'A 2/m','A 2/a','A 2/n',),
3371	'Pmmm':('P 2 2 2',
3372	'P 2 2 21','P 21 2 2','P 2 21 2',
3373	'P 21 21 2','P 2 21 21','P 21 2 21',
3374	'P 21 21 21',
3375	'P m m 2','P 2 m m','P m 2 m',
3376	'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3377	'P c c 2','P 2 a a','P b 2 b',
3378	'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3379	'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3380	'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3381	'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3382	'P b a 2','P 2 c b','P c 2 a',
3383	'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3384	'P n n 2','P 2 n n','P n 2 n',
3385	'P m m m','P n n n',
3386	'P c c m','P m a a','P b m b',
3387	'P b a n','P n c b','P c n a',
3388	'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3389	'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3390	'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3391	'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3392	'P b a m','P m c b','P c m a',
3393	'P c c n','P n a a','P b n b',
3394	'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3395	'P n n m','P m n n','P n m n',
3396	'P m m n','P n m m','P m n m',
3397	'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3398	'P b c a','P c a b',
3399	'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3400	),
3401	'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3402	'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3403	'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3404	'C m c m','C c m m','C m c a','C c m b',
3405	'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3406	'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3407	'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3408	'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3409	'A m m a','A m a m','A b m a','A c a m',
3410	'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3411	'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3412	'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3413	'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3414	'B b m m','B m m b','B b c m','B m a b',
3415	'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3416	'Immm':('I 2 2 2','I 21 21 21',
3417	'I m m 2','I m 2 m','I 2 m m',
3418	'I b a 2','I 2 c b','I c 2 a',
3419	'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3420	'I m m m','I b a m','I m c b','I c m a',
3421	'I b c a','I c a b',
3422	'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3423	'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3424	'F m m 2','F m 2 m','F 2 m m',
3425	'F d d 2','F d 2 d','F 2 d d',),
3426	'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3427	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3428	'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3429	'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3430	'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3431	'P -4 n 2','P 4/m m m','P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m',
3432	'P 4/m n c','P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
3433	'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m','P 42/n m c',
3434	'P 42/n c m',),
3435	'I4/mmm':('I 4','I 41','I -4','I 4/m','I 41/a','I 4 2 2','I 41 2 2','I 4 m m',
3436	'I 4 c m','I 41 m d','I 41 c d',
3437	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d','I 4/m m m','I 4/m c m',
3438	'I 41/a m d','I 41/a c d'),
3439	'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3440	'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3441	'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3442	'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3443	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3444	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3445	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3446	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3447	'Pm3m': ('P 2 3','P 21 3','P m 3','P m -3','P n 3','P n -3','P a 3','P a -3',
3448	'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3449	'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3450	'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3451	'Im3m':('I 2 3','I 21 3','I m 3','I m -3','I a 3','I a -3', 'I 4 3 2','I 41 3 2',
3452	'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a 3 d','I a -3 d','I n 3 n','I n -3 n'),
3453	'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3454	'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3455	'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3456	'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3457	}
3458	sgequiv_2002_orthorhombic = {}
3459	''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3460	along with the pre-2002 name. The e designates a double glide-plane
3461	'''
3462	sgequiv_2002_orthorhombic = {
3463	'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3464	'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3465	'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3466	'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3467	'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3468	'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3469	'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3470	'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3471
3472	#'A few non-standard space groups for test use'
3473	nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h',
3474	'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3475	'R 3 c r','R -3 c r','R -3 m r',),
3476
3477	#Use the space groups types in this order to list the symbols in the
3478	#order they are listed in the International Tables, vol. A'''
3479	symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal',
3480	'trigonal', 'hexagonal', 'cubic')
3481
3482	# self-test materials follow. Requires files in directory testinp
3483	selftestlist = []
3484	'''Defines a list of self-tests'''
3485	selftestquiet = True
3486	def _ReportTest():
3487	'Report name and doc string of current routine when ``selftestquiet`` is False'
3488	if not selftestquiet:
3489	import inspect
3490	caller = inspect.stack()[1][3]
3491	doc = eval(caller).__doc__
3492	if doc is not None:
3493	print('testing '+__file__+' with '+caller+' ('+doc+')')
3494	else:
3495	print('testing '+__file__()+" with "+caller)
3496	def test0():
3497	'''self-test #0: exercise MoveToUnitCell'''
3498	_ReportTest()
3499	msg = "MoveToUnitCell failed"
3500	assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3501	assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3502	assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3503	assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3504	selftestlist.append(test0)
3505
3506	def test1():
3507	'''self-test #1: SpcGroup against previous results'''
3508	#'''self-test #1: SpcGroup and SGPrint against previous results'''
3509	_ReportTest()
3510	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3511	if ospath.exists(testdir):
3512	if testdir not in sys.path: sys.path.insert(0,testdir)
3513	import spctestinp
3514	def CompareSpcGroup(spc, referr, refdict, reflist):
3515	'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3516	# if an error is reported, the dictionary can be ignored
3517	msg0 = "CompareSpcGroup failed on space group %s" % spc
3518	result = SpcGroup(spc)
3519	if result[0] == referr and referr > 0: return True
3520	# #print result[1]['SpGrp']
3521	#msg = msg0 + " in list lengths"
3522	#assert len(keys) == len(refdict.keys()), msg
3523	for key in refdict.keys():
3524	if key == 'SGOps' or key == 'SGCen':
3525	msg = msg0 + (" in key %s length" % key)
3526	assert len(refdict[key]) == len(result[1][key]), msg
3527	for i in range(len(refdict[key])):
3528	msg = msg0 + (" in key %s level 0" % key)
3529	assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
3530	msg = msg0 + (" in key %s level 1" % key)
3531	assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
3532	else:
3533	msg = msg0 + (" in key %s" % key)
3534	assert result[1][key] == refdict[key], msg
3535	msg = msg0 + (" in key %s reflist" % key)
3536	#for (l1,l2) in zip(reflist, SGPrint(result[1])):
3537	# assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
3538	# for now disable SGPrint testing, output has changed
3539	#assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
3540	for spc in spctestinp.SGdat:
3541	CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
3542	selftestlist.append(test1)
3543
3544	def test2():
3545	'''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
3546	_ReportTest()
3547	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3548	if ospath.exists(testdir):
3549	if testdir not in sys.path: sys.path.insert(0,testdir)
3550	import sgtbxtestinp
3551	def CompareWcctbx(spcname, cctbx_in, debug=0):
3552	'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
3553	cctbx = cctbx_in[:] # make copy so we don't delete from the original
3554	spc = (SpcGroup(spcname))[1]
3555	if debug: print (spc['SpGrp'])
3556	if debug: print (spc['SGCen'])
3557	latticetype = spcname.strip().upper()[0]
3558	# lattice type of R implies Hexagonal centering", fix the rhombohedral settings
3559	if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
3560	assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
3561	onebar = [1]
3562	if spc['SGInv']: onebar.append(-1)
3563	for (op,off) in spc['SGOps']:
3564	for inv in onebar:
3565	for cen in spc['SGCen']:
3566	noff = off + cen
3567	noff = MoveToUnitCell(noff)[0]
3568	mult = tuple((op*inv).ravel().tolist())
3569	if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
3570	for refop in cctbx:
3571	if debug: print (refop)
3572	# check the transform
3573	if refop[:9] != mult: continue
3574	if debug: print ("mult match")
3575	# check the translation
3576	reftrans = list(refop[-3:])
3577	reftrans = MoveToUnitCell(reftrans)[0]
3578	if all(abs(noff - reftrans) < 1.e-5):
3579	cctbx.remove(refop)
3580	break
3581	else:
3582	assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
3583	for key in sgtbxtestinp.sgtbx:
3584	CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
3585	selftestlist.append(test2)
3586
3587	def test3():
3588	'''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
3589	for selected space groups against info in IT Volume A '''
3590	_ReportTest()
3591	def ExerciseSiteSym (spc, crdlist):
3592	'compare site symmetries and multiplicities for a specified space group'
3593	msg = "failed on site sym test for %s" % spc
3594	(E,S) = SpcGroup(spc)
3595	assert not E, msg
3596	for t in crdlist:
3597	symb, m, n, od = SytSym(t[0],S)
3598	if symb.strip() != t[2].strip() or m != t[1]:
3599	print (spc,t[0],m,n,symb,t[2],od)
3600	assert m == t[1]
3601	#assert symb.strip() == t[2].strip()
3602
3603	ExerciseSiteSym('p 1',[
3604	((0.13,0.22,0.31),1,'1'),
3605	((0,0,0),1,'1'),
3606	])
3607	ExerciseSiteSym('p -1',[
3608	((0.13,0.22,0.31),2,'1'),
3609	((0,0.5,0),1,'-1'),
3610	])
3611	ExerciseSiteSym('C 2/c',[
3612	((0.13,0.22,0.31),8,'1'),
3613	((0.0,.31,0.25),4,'2(y)'),
3614	((0.25,.25,0.5),4,'-1'),
3615	((0,0.5,0),4,'-1'),
3616	])
3617	ExerciseSiteSym('p 2 2 2',[
3618	((0.13,0.22,0.31),4,'1'),
3619	((0,0.5,.31),2,'2(z)'),
3620	((0.5,.31,0.5),2,'2(y)'),
3621	((.11,0,0),2,'2(x)'),
3622	((0,0.5,0),1,'222'),
3623	])
3624	ExerciseSiteSym('p 4/n',[
3625	((0.13,0.22,0.31),8,'1'),
3626	((0.25,0.75,.31),4,'2(z)'),
3627	((0.5,0.5,0.5),4,'-1'),
3628	((0,0.5,0),4,'-1'),
3629	((0.25,0.25,.31),2,'4(001)'),
3630	((0.25,.75,0.5),2,'-4(001)'),
3631	((0.25,.75,0.0),2,'-4(001)'),
3632	])
3633	ExerciseSiteSym('p 31 2 1',[
3634	((0.13,0.22,0.31),6,'1'),
3635	((0.13,0.0,0.833333333),3,'2(100)'),
3636	((0.13,0.13,0.),3,'2(110)'),
3637	])
3638	ExerciseSiteSym('R 3 c',[
3639	((0.13,0.22,0.31),18,'1'),
3640	((0.0,0.0,0.31),6,'3'),
3641	])
3642	ExerciseSiteSym('R 3 c R',[
3643	((0.13,0.22,0.31),6,'1'),
3644	((0.31,0.31,0.31),2,'3(111)'),
3645	])
3646	ExerciseSiteSym('P 63 m c',[
3647	((0.13,0.22,0.31),12,'1'),
3648	((0.11,0.22,0.31),6,'m(100)'),
3649	((0.333333,0.6666667,0.31),2,'3m(100)'),
3650	((0,0,0.31),2,'3m(100)'),
3651	])
3652	ExerciseSiteSym('I a -3',[
3653	((0.13,0.22,0.31),48,'1'),
3654	((0.11,0,0.25),24,'2(x)'),
3655	((0.11,0.11,0.11),16,'3(111)'),
3656	((0,0,0),8,'-3(111)'),
3657	])
3658	selftestlist.append(test3)
3659
3660	if __name__ == '__main__':
3661	# run self-tests
3662	selftestquiet = False
3663	for test in selftestlist:
3664	test()
3665	print ("OK")

Note: See TracBrowser for help on using the repository browser.

Download in other formats: