Context Navigation

source: trunk/GSASIIspc.py @ 3387

Last change on this file since 3387 was 3387, checked in by vondreele, 5 years ago
fix to GetLittleGrpOps?
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 153.8 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-15 13:33:56 +0000 (Tue, 15 May 2018) $
12	# $Author: vondreele $
13	# $Revision: 3387 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 3387 2018-05-15 13:33:56Z 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: 3387 $")
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 GetLittleGrpOps(SGData,vec):
1730	''' Find rotation part of operators that leave vec unchanged
1731
1732	:param SGData: space group data structure as defined in SpcGroup above.
1733	:param vec: a numpy array of fractional vector coordinates
1734	:returns: Little - list of operators [M,T] that form the little gropu
1735	'''
1736	Little = []
1737	Ops = SGData['SGOps'][:]
1738	if SGData['SGInv']:
1739	Ops += [[-M,-T] for [M,T] in Ops]
1740	for [M,T] in Ops:
1741	tvec = np.inner(M,vec)%1.
1742	if np.allclose(tvec,vec%1.):
1743	Little.append([M,T])
1744	return Little
1745
1746	def SGProd(OpA,OpB):
1747	'''
1748	Form space group operator product. OpA & OpB are [M,V] pairs;
1749	both must be of same dimension (3 or 4). Returns [M,V] pair
1750	'''
1751	A,U = OpA
1752	B,V = OpB
1753	M = np.inner(B,A.T)
1754	W = np.inner(B,U)+V
1755	return M,W
1756
1757	def MoveToUnitCell(xyz):
1758	'''
1759	Translates a set of coordinates so that all values are >=0 and < 1
1760
1761	:param xyz: a list or numpy array of fractional coordinates
1762	:returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1763	'''
1764	XYZ = (np.array(xyz)+10.)%1.
1765	cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1766	return XYZ,cell
1767
1768	def Opposite(XYZ,toler=0.0002):
1769	'''
1770	Gives opposite corner, edge or face of unit cell for position within tolerance.
1771	Result may be just outside the cell within tolerance
1772
1773	:param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1774	:param toler: unit cell fraction tolerance making opposite
1775	:returns:
1776	XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1777	'''
1778	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]]
1779	TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1780	perm = TB*perm3
1781	cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1782	D = dict(zip(cperm,perm))
1783	new = {}
1784	for key in D:
1785	new[key] = np.array(D[key])+np.array(XYZ)
1786	return new
1787
1788	def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1789	'''
1790	Generates the equivalent positions for a specified coordinate and space group
1791
1792	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
1793	:param SGData: from :func:`SpcGroup`
1794	:param All: True return all equivalent positions including duplicates;
1795	False return only unique positions
1796	:param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1797	:param Move: True move generated atom positions to be inside cell
1798	False do not move atoms
1799	:return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
1800
1801	* [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1802	* Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1803	* is centering operator number (1-4) and - is for inversion
1804	Cell = unit cell translations needed to put new positions inside cell
1805	[UijEquiv] - equivalent Uij; absent if no Uij given
1806	* +1/-1 for spin inversion of operator - empty if not magnetic
1807
1808	'''
1809	XYZEquiv = []
1810	UijEquiv = []
1811	Idup = []
1812	Cell = []
1813	inv = int(SGData['SGInv']+1)
1814	icen = SGData['SGCen']
1815	if SGData['SGFixed']:
1816	inv = 1
1817	SpnFlp = SGData.get('SpnFlp',[])
1818	spnflp = []
1819	X = np.array(XYZ)
1820	mj = 0
1821	if Move:
1822	X = MoveToUnitCell(X)[0]
1823	for ic,cen in enumerate(icen):
1824	C = np.array(cen)
1825	for invers in range(inv):
1826	for io,[M,T] in enumerate(SGData['SGOps']):
1827	idup = ((io+1)+100ic)(1-2*invers)
1828	XT = np.inner(M,X)+T
1829	if len(Uij):
1830	U = Uij2U(Uij)
1831	U = np.inner(M,np.inner(U,M).T)
1832	newUij = U2Uij(U)
1833	if invers:
1834	XT = -XT
1835	XT += C
1836	cell = np.zeros(3,dtype=np.int32)
1837	cellj = np.zeros(3,dtype=np.int32)
1838	if Move:
1839	newX,cellj = MoveToUnitCell(XT)
1840	else:
1841	newX = XT
1842	cell += cellj
1843	if All:
1844	if np.allclose(newX,X,atol=0.0002):
1845	idup = False
1846	else:
1847	if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1848	idup = False
1849	if All or idup:
1850	XYZEquiv.append(newX)
1851	Idup.append(idup)
1852	Cell.append(cell)
1853	if len(Uij):
1854	UijEquiv.append(newUij)
1855	if len(SpnFlp):
1856	spnflp.append(SpnFlp[mj])
1857	else:
1858	spnflp.append(1)
1859	mj += 1
1860	if len(Uij):
1861	return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
1862	else:
1863	return zip(XYZEquiv,Idup,Cell,spnflp)
1864
1865	def GenHKL(HKL,SGData):
1866	''' Generates all equivlent reflections including Friedel pairs
1867	:param HKL: [h,k,l] must be integral values
1868	:param SGData: space group data obtained from SpcGroup
1869	:returns: array Uniq: equivalent reflections
1870	'''
1871
1872	Ops = SGData['SGOps']
1873	OpM = np.array([op[0] for op in Ops])
1874	Uniq = np.inner(OpM,HKL)
1875	Uniq = list(Uniq)+list(-1*Uniq)
1876	return np.array(Uniq)
1877
1878	def GenHKLf(HKL,SGData):
1879	'''
1880	Uses old GSAS Fortran routine genhkl.for
1881
1882	:param HKL: [h,k,l] must be integral values for genhkl.for to work
1883	:param SGData: space group data obtained from SpcGroup
1884	:returns: iabsnt,mulp,Uniq,phi
1885
1886	* iabsnt = True if reflection is forbidden by symmetry
1887	* mulp = reflection multiplicity including Friedel pairs
1888	* Uniq = numpy array of equivalent hkl in descending order of h,k,l
1889	* phi = phase offset for each equivalent h,k,l
1890
1891	'''
1892	hklf = list(HKL)+[0,] #could be numpy array!
1893	Ops = SGData['SGOps']
1894	OpM = np.array([op[0] for op in Ops],order='F')
1895	OpT = np.array([op[1] for op in Ops])
1896	Cen = np.array([cen for cen in SGData['SGCen']],order='F')
1897
1898	import pyspg
1899	Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1900	h,k,l,f = Uniq
1901	Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
1902	phi = f[:Nuniq]
1903	return iabsnt,mulp,Uniq,phi
1904
1905	def MagHKLchk(HKL,SGData):
1906	SpnFlp = SGData['SpnFlp']
1907	print(HKL)
1908	Uniq = GenHKL(HKL,SGData)
1909
1910	def checkSSLaue(HKL,SGData,SSGData):
1911	#Laue check here - Toss HKL if outside unique Laue part
1912	h,k,l,m = HKL
1913	if SGData['SGLaue'] == '2/m':
1914	if SGData['SGUniq'] == 'a':
1915	if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
1916	return False
1917	elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
1918	return False
1919	else:
1920	return True
1921	elif SGData['SGUniq'] == 'b':
1922	if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
1923	return False
1924	elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
1925	return False
1926	else:
1927	return True
1928	elif SGData['SGUniq'] == 'c':
1929	if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
1930	return False
1931	elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
1932	return False
1933	else:
1934	return True
1935	elif SGData['SGLaue'] == 'mmm':
1936	if 'a' in SSGData['modSymb']:
1937	if h == 0 and m < 0:
1938	return False
1939	else:
1940	return True
1941	elif 'b' in SSGData['modSymb']:
1942	if k == 0 and m < 0:
1943	return False
1944	else:
1945	return True
1946	elif 'g' in SSGData['modSymb']:
1947	if l == 0 and m < 0:
1948	return False
1949	else:
1950	return True
1951	else: #tetragonal, trigonal, hexagonal (& triclinic?)
1952	if l == 0 and m < 0:
1953	return False
1954	else:
1955	return True
1956
1957
1958	def checkSSextc(HKL,SSGData):
1959	Ops = SSGData['SSGOps']
1960	OpM = np.array([op[0] for op in Ops])
1961	OpT = np.array([op[1] for op in Ops])
1962	HKLS = np.array([HKL,-HKL]) #Freidel's Law
1963	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
1964	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
1965	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
1966	if dhkl.any():
1967	continue
1968	else:
1969	if phkl%1.:
1970	return False
1971	return True
1972
1973	################################################################################
1974	#### Site symmetry tables
1975	################################################################################
1976
1977	OprName = {
1978	'-6643': ['-1',1],'6479' : ['2(z)',2],'-6479': ['m(z)',3],
1979	'6481' : ['m(y)',4],'-6481': ['2(y)',5],'6641' : ['m(x)',6],
1980	'-6641': ['2(x)',7],'6591' : ['m(+-0)',8],'-6591': ['2(+-0)',9],
1981	'6531' : ['m(110)',10],'-6531': ['2(110)',11],'6537' : ['4(z)',12],
1982	'-6537': ['-4(z)',13],'975' : ['3(111)',14],'6456' : ['3',15],
1983	'-489' : ['3(+--)',16],'483' : ['3(-+-)',17],'-969' : ['3(--+)',18],
1984	'819' : ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' : ['m(0+-)',21],
1985	'-2431': ['2(0+-)',22],'-657' : ['m(xz)',23],'657' : ['2(xz)',24],
1986	'1943' : ['-4(x)',25],'-1943': ['4(x)',26],'-2429': ['m(yz)',27],
1987	'2429' : ['2(yz)',28],'639' : ['-4(y)',29],'-639' : ['4(y)',30],
1988	'-6484': ['2(010)',4],'6484' : ['m(010)',5],'-6668': ['2(100)',6],
1989	'6668' : ['m(100)',7],'-6454': ['2(120)',18],'6454' : ['m(120)',19],
1990	'-6638': ['2(210)',20],'6638' : ['m(210)',21], #search in SytSym ends at m(210)
1991	'2223' : ['3(+++)2',39],
1992	'6538' : ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
1993	'2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489' :['-3(+--)1',45],
1994	'801' : ['4(y)1',46],'1945' : ['4(x)3',47],'-6585': ['-4(z)3 ',48],
1995	'6585' : ['4(z)3',49],'6584' : ['3(z)2',50],'6666' : ['6(z)5 ',51],
1996	'6643' : ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
1997	'-6666': ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
1998	'-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
1999	'969' :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2000	'-2151':['-3(+--)2',64], }
2001
2002	KNsym = {
2003	'0' :' 1 ','1' :' -1 ','64' :' 2(x)','32' :' m(x)',
2004	'97' :' 2/m(x)','16' :' 2(y)','8' :' m(y)','25' :' 2/m(y)',
2005	'2' :' 2(z)','4' :' m(z)','7' :' 2/m(z)','134217728' :' 2(yz)',
2006	'67108864' :' m(yz)','201326593' :' 2/m(yz)','2097152' :' 2(0+-)','1048576' :' m(0+-)',
2007	'3145729' :'2/m(0+-)','8388608' :' 2(xz)','4194304' :' m(xz)','12582913' :' 2/m(xz)',
2008	'524288' :' 2(+0-)','262144' :' m(+0-)','796433' :'2/m(+0-)','1024' :' 2(xy)',
2009	'512' :' m(xy)','1537' :' 2/m(xy)','256' :' 2(+-0)','128' :' m(+-0)',
2010	'385' :'2/m(+-0)','76' :' mm2(x)','52' :' mm2(y)','42' :' mm2(z)',
2011	'135266336' :' mm2(yz)','69206048' :'mm2(0+-)','8650760' :' mm2(xz)','4718600' :'mm2(+0-)',
2012	'1156' :' mm2(xy)','772' :'mm2(+-0)','82' :' 222 ','136314944' :' 222(x)',
2013	'8912912' :' 222(y)','1282' :' 222(z)','127' :' mmm ','204472417' :' mmm(x)',
2014	'13369369' :' mmm(y)','1927' :' mmm(z)','33554496' :' 4(x)','16777280' :' -4(x)',
2015	'50331745' :'4/m(x)' ,'169869394' :'422(x)','84934738' :'-42m(x)','101711948' :'4mm(x)',
2016	'254804095' :'4/mmm(x)','536870928 ':' 4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2017	'545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2018	'2050' :' 4(z)','4098' :' -4(z)','6151' :'4/m(z)','3410' :'422(z)',
2019	'4818' :'-42m(z)','2730' :'4mm(z)','8191' :'4/mmm(z)','8192' :' 3(111)',
2020	'8193' :' -3(111)','2629888' :' 32(111)','1319040' :' 3m(111)','3940737' :'-3m(111)',
2021	'32768' :' 3(+--)','32769' :' -3(+--)','10519552' :' 32(+--)','5276160' :' 3m(+--)',
2022	'15762945' :'-3m(+--)','65536' :' 3(-+-)','65537' :' -3(-+-)','134808576' :' 32(-+-)',
2023	'67437056' :' 3m(-+-)','202180097' :'-3m(-+-)','131072' :' 3(--+)','131073' :' -3(--+)',
2024	'142737664' :' 32(--+)','71434368' :' 3m(--+)','214040961' :'-3m(--+)','237650' :' 23 ',
2025	'237695' :' m3 ','715894098' :' 432 ','358068946' :' -43m ','1073725439':' m3m ',
2026	'68157504' :' mm2(d100)','4456464' :' mm2(d010)','642' :' mm2(d001)','153092172' :'-4m2(x)',
2027	'277348404' :'-4m2(y)','5418' :'-4m2(z)','1075726335':' 6/mmm ','1074414420':'-6m2(100)',
2028	'1075070124':'-6m2(120)','1075069650':' 6mm ','1074414890':' 622 ','1073758215':' 6/m ',
2029	'1073758212':' -6 ','1073758210':' 6 ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2030	'1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2031	'1073758209':' -3 ','1073758208':' 3 ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2032	'1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2033	'1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2034	'1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2035	'1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2036	'1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2037	'1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2038	'1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2039	'1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2040	'1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':' -1 ',
2041	}
2042
2043	NXUPQsym = {
2044	'1' :(28,29,28,28),'-1' :( 1,29,28, 0),'2(x)' :(12,18,12,25),'m(x)' :(25,18,12,25),
2045	'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),
2046	'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),
2047	'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),
2048	'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),
2049	'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),
2050	'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),
2051	'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),
2052	'mm2(yz)' :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)' :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2053	'mm2(xy)' :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222' :( 1,10, 0,-1),'222(x)' :( 1,13, 0,-1),
2054	'222(y)' :( 1,12, 0,-1),'222(z)' :( 1,11, 0,-1),'mmm' :( 1,10, 0,-1),'mmm(x)' :( 1,13, 0,-1),
2055	'mmm(y)' :( 1,12, 0,-1),'mmm(z)' :( 1,11, 0,-1),'4(x)' :(12, 4,12, 0),'-4(x)' :( 1, 4,12, 0),
2056	'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),
2057	'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),
2058	'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,),
2059	'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),
2060	'-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),
2061	'-3(111)' :( 1, 5, 2, 0),'32(111)' :( 1, 5, 0, 2),'3m(111)' :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2062	'3(+--)' :( 5, 8, 5, 0),'-3(+--)' :( 1, 8, 5, 0),'32(+--)' :( 1, 8, 0, 5),'3m(+--)' :( 5, 8, 0, 5),
2063	'-3m(+--)' :( 1, 8, 0,-1),'3(-+-)' :( 4, 7, 4, 0),'-3(-+-)' :( 1, 7, 4, 0),'32(-+-)' :( 1, 7, 0, 4),
2064	'3m(-+-)' :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)' :( 3, 6, 3, 0),'-3(--+)' :( 1, 6, 3, 0),
2065	'32(--+)' :( 1, 6, 0, 3),'3m(--+)' :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23' :( 1, 1, 0, 0),
2066	'm3' :( 1, 1, 0, 0),'432' :( 1, 1, 0, 0),'-43m' :( 1, 1, 0, 0),'m3m' :( 1, 1, 0, 0),
2067	'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)' :( 1, 4, 0,-1),
2068	'-4m2(y)' :( 1, 3, 0,-1),'-4m2(z)' :( 1, 2, 0,-1),'6/mmm' :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2069	'-6m2(120)':( 1, 9, 0,-1),'6mm' :(14, 9, 0,-1),'622' :( 1, 9, 0,-1),'6/m' :( 1, 9,14,-1),
2070	'-6' :( 1, 9,14, 0),'6' :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2071	'3m(100)' :(14, 9, 0,14),'3m(120)' :(14, 9, 0,14),'32(100)' :( 1, 9, 0,14),'32(120)' :( 1, 9, 0,14),
2072	'-3' :( 1, 9,14, 0),'3' :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2073	'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2074	'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2075	'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2076	'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),
2077	'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),
2078	'm(001)' :(23,16,14,23),'m(100)' :(26,25,12,26),'m(010)' :(27,28,13,27),'m(110)' :(18,19, 6,18),
2079	'm(120)' :(24,27,15,24),'m(210)' :(25,26,16,25),'m(+-0)' :(17,20, 7,17),'2(001)' :(14,16,14,23),
2080	'2(100)' :(12,25,12,26),'2(010)' :(13,28,13,27),'2(110)' :( 6,19, 6,18),'2(120)' :(15,27,15,24),
2081	'2(210)' :(16,26,16,25),'2(+-0)' :( 7,20, 7,17),'-1' :( 1,29,28, 0)
2082	}
2083
2084	CSxinel = [[], # 0th empty - indices are Fortran style
2085	[[0,0,0],[ 0.0, 0.0, 0.0]], #1 0 0 0
2086	[[1,1,1],[ 1.0, 1.0, 1.0]], #2 X X X
2087	[[1,1,1],[ 1.0, 1.0,-1.0]], #3 X X -X
2088	[[1,1,1],[ 1.0,-1.0, 1.0]], #4 X -X X
2089	[[1,1,1],[ 1.0,-1.0,-1.0]], #5 -X X X
2090	[[1,1,0],[ 1.0, 1.0, 0.0]], #6 X X 0
2091	[[1,1,0],[ 1.0,-1.0, 0.0]], #7 X -X 0
2092	[[1,0,1],[ 1.0, 0.0, 1.0]], #8 X 0 X
2093	[[1,0,1],[ 1.0, 0.0,-1.0]], #9 X 0 -X
2094	[[0,1,1],[ 0.0, 1.0, 1.0]], #10 0 Y Y
2095	[[0,1,1],[ 0.0, 1.0,-1.0]], #11 0 Y -Y
2096	[[1,0,0],[ 1.0, 0.0, 0.0]], #12 X 0 0
2097	[[0,1,0],[ 0.0, 1.0, 0.0]], #13 0 Y 0
2098	[[0,0,1],[ 0.0, 0.0, 1.0]], #14 0 0 Z
2099	[[1,1,0],[ 1.0, 2.0, 0.0]], #15 X 2X 0
2100	[[1,1,0],[ 2.0, 1.0, 0.0]], #16 2X X 0
2101	[[1,1,2],[ 1.0, 1.0, 1.0]], #17 X X Z
2102	[[1,1,2],[ 1.0,-1.0, 1.0]], #18 X -X Z
2103	[[1,2,1],[ 1.0, 1.0, 1.0]], #19 X Y X
2104	[[1,2,1],[ 1.0, 1.0,-1.0]], #20 X Y -X
2105	[[1,2,2],[ 1.0, 1.0, 1.0]], #21 X Y Y
2106	[[1,2,2],[ 1.0, 1.0,-1.0]], #22 X Y -Y
2107	[[1,2,0],[ 1.0, 1.0, 0.0]], #23 X Y 0
2108	[[1,0,2],[ 1.0, 0.0, 1.0]], #24 X 0 Z
2109	[[0,1,2],[ 0.0, 1.0, 1.0]], #25 0 Y Z
2110	[[1,1,2],[ 1.0, 2.0, 1.0]], #26 X 2X Z
2111	[[1,1,2],[ 2.0, 1.0, 1.0]], #27 2X X Z
2112	[[1,2,3],[ 1.0, 1.0, 1.0]], #28 X Y Z
2113	]
2114
2115	CSuinel = [[], # 0th empty - indices are Fortran style
2116	[[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
2117	[[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
2118	[[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
2119	[[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
2120	[[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
2121	[[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
2122	[[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
2123	[[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
2124	[[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
2125	[[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
2126	[[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
2127	[[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
2128	[[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
2129	[[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
2130	[[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
2131	[[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
2132	[[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
2133	[[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
2134	[[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
2135	[[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
2136	[[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
2137	[[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
2138	[[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
2139	[[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
2140	[[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
2141	[[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
2142	[[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
2143	[[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
2144	[[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
2145	]
2146
2147	################################################################################
2148	#### Site symmetry routines
2149	################################################################################
2150
2151	def GetOprPtrName(key):
2152	'Needs a doc string'
2153	oprName = OprName[key][0]
2154	return oprName.replace('(','').replace(')','')
2155
2156	def GetOprPtrNumber(key):
2157	'Needs a doc string'
2158	return OprName[key][1]
2159
2160	def GetOprName(key):
2161	'Needs a doc string'
2162	return OprName[key][0]
2163
2164	def GetKNsym(key):
2165	'Needs a doc string'
2166	try:
2167	return KNsym[key].strip()
2168	except KeyError:
2169	return 'sp'
2170
2171	def GetNXUPQsym(siteSym):
2172	'''
2173	The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2174	'''
2175	return NXUPQsym[siteSym]
2176
2177	def GetCSxinel(siteSym):
2178	"returns Xyz terms, multipliers, GUI flags"
2179	indx = GetNXUPQsym(siteSym.strip())
2180	return CSxinel[indx[0]]
2181
2182	def GetCSuinel(siteSym):
2183	"returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2184	indx = GetNXUPQsym(siteSym.strip())
2185	return CSuinel[indx[1]]
2186
2187	def GetCSpqinel(SpnFlp,dupDir):
2188	"returns Mxyz terms, multipliers, GUI flags"
2189	CSI = [[1,2,3],[1.0,1.0,1.0]]
2190	for sopr in dupDir:
2191	# print (sopr,dupDir[sopr])
2192	opr = sopr.replace("'",'')
2193	indx = GetNXUPQsym(opr)
2194	if SpnFlp[dupDir[sopr]] > 0:
2195	csi = CSxinel[indx[2]] #P
2196	else:
2197	csi = CSxinel[indx[3]] #Q
2198	# print(opr,indx,csi,CSI)
2199	if not len(csi):
2200	return [[0,0,0],[0.,0.,0.]]
2201	for kcs in [0,1,2]:
2202	if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2203	jcs = CSI[0][kcs]
2204	for ics in [0,1,2]:
2205	if CSI[0][ics] == jcs:
2206	CSI[0][ics] = 0
2207	CSI[1][ics] = 0.
2208	elif CSI[0][ics] > jcs:
2209	CSI[0][ics] = CSI[0][ics]-1
2210	elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2211	CSI[1][kcs] = csi[1][kcs]
2212	elif CSI[0][kcs] >= csi[0][kcs]:
2213	CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2214	if CSI[1][kcs] != csi[1][kcs]:
2215	if CSI[1][kcs] == 1.:
2216	CSI[1][kcs] = csi[1][kcs]
2217	# print(CSI)
2218	return CSI
2219
2220	def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2221	phase = np.sum(XYZ*wave)
2222	ssopinv = nl.inv(ssop[0])
2223	mst = ssopinv[3][:3]
2224	epsinv = ssopinv[3][3]
2225	sdet = nl.det(sop[0])
2226	ssdet = nl.det(ssop[0])
2227	dtau = mst(XYZ-sop[1])-epsinvssop[1][3]
2228	dT = 1.0
2229	if np.any(dtau%.5):
2230	sumdtau = np.sum(dtau%.5)
2231	dT = 0.
2232	if np.abs(sumdtau-.5) > 1.e-4:
2233	dT = np.tan(np.pi*sumdtau)
2234	tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2235	return sdet,ssdet,dtau,dT,tauT
2236
2237	def OpsfromStringOps(A,SGData,SSGData):
2238	SGOps = SGData['SGOps']
2239	SSGOps = SSGData['SSGOps']
2240	Ax = A.split('+')
2241	Ax[0] = int(Ax[0])
2242	iC = 1
2243	if Ax[0] < 0:
2244	iC = -1
2245	Ax[0] = abs(Ax[0])
2246	nA = Ax[0]%100-1
2247	nC = Ax[0]//100
2248	unit = [0,0,0]
2249	if len(Ax) > 1:
2250	unit = eval('['+Ax[1]+']')
2251	return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2252
2253	def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2254
2255	def orderParms(CSI):
2256	parms = [0,]
2257	for csi in CSI:
2258	for i in [0,1,2]:
2259	if csi[i] not in parms:
2260	parms.append(csi[i])
2261	for csi in CSI:
2262	for i in [0,1,2]:
2263	csi[i] = parms.index(csi[i])
2264	return CSI
2265
2266	def fracCrenel(tau,Toff,Twid):
2267	Tau = (tau-Toff[:,np.newaxis])%1.
2268	A = np.where(Tau<Twid[:,np.newaxis],1.,0.)
2269	return A
2270
2271	def fracFourier(tau,nH,fsin,fcos):
2272	SA = np.sin(2.nHnp.pi*tau)
2273	CB = np.cos(2.nHnp.pi*tau)
2274	A = SA[np.newaxis,np.newaxis,:]*fsin[:,:,np.newaxis]
2275	B = CB[np.newaxis,np.newaxis,:]*fcos[:,:,np.newaxis]
2276	return A+B
2277
2278	def posFourier(tau,nH,psin,pcos):
2279	SA = np.sin(2nHnp.pi*tau)
2280	CB = np.cos(2nHnp.pi*tau)
2281	A = SA[np.newaxis,np.newaxis,:]*psin[:,:,np.newaxis]
2282	B = CB[np.newaxis,np.newaxis,:]*pcos[:,:,np.newaxis]
2283	return A+B
2284
2285	def posSawtooth(tau,Toff,slopes):
2286	Tau = (tau-Toff)%1.
2287	A = slopes[:,np.newaxis]*Tau
2288	return A
2289
2290	def posZigZag(tau,Tmm,XYZmax):
2291	DT = Tmm[1]-Tmm[0]
2292	slopeUp = 2.*XYZmax/DT
2293	slopeDn = 2.*XYZmax/(1.-DT)
2294	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])
2295	return A
2296
2297	def posBlock(tau,Tmm,XYZmax):
2298	A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2299	return A
2300
2301	def DoFrac():
2302	dF = None
2303	dFTP = None
2304	if siteSym == '1':
2305	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2306	elif siteSym == '-1':
2307	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2308	else:
2309	delt2 = np.eye(2)*0.001
2310	FSC = np.ones(2,dtype='i')
2311	CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2312	if 'Crenel' in waveType:
2313	dF = np.zeros_like(tau)
2314	else:
2315	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2316	dFT = np.zeros_like(dF)
2317	dFTP = []
2318	for i in SdIndx:
2319	sop = Sop[i]
2320	ssop = SSop[i]
2321	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2322	fsc = np.ones(2,dtype='i')
2323	if 'Crenel' in waveType:
2324	dFT = np.zeros_like(tau)
2325	fsc = [1,1]
2326	else: #Fourier
2327	dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2328	dFT = nl.det(sop[0])*dFT
2329	dFT = dFT[:,np.argsort(tauT)]
2330	dFT[0] *= ssdet
2331	dFT[1] *= sdet
2332	dFTP.append(dFT)
2333
2334	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2335	fsc = [1,1]
2336	if dT:
2337	CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2338	else:
2339	CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2340	FSC = np.zeros(2,dtype='i')
2341	return CSI,dF,dFTP
2342	else:
2343	for i in range(2):
2344	if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2345	fsc[i] = 1
2346	else:
2347	fsc[i] = 0
2348	FSC &= fsc
2349	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2350	n = -1
2351	for i,F in enumerate(FSC):
2352	if F:
2353	n += 1
2354	CSI[0][i] = n+1
2355	CSI[1][i] = 1.0
2356
2357	return CSI,dF,dFTP
2358
2359	def DoXYZ():
2360	dX,dXTP = None,None
2361	if siteSym == '1':
2362	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2363	elif siteSym == '-1':
2364	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.],]
2365	else:
2366	delt4 = np.ones(4)*0.001
2367	delt5 = np.ones(5)*0.001
2368	delt6 = np.eye(6)*0.001
2369	if 'Fourier' in waveType:
2370	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2371	#3x6x12 modulated position array (X,Spos,tau)& force positive
2372	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2373	elif waveType == 'Sawtooth':
2374	dX = posSawtooth(tau,delt4[0],delt4[1:])
2375	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2376	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2377	elif waveType in ['ZigZag','Block']:
2378	if waveType == 'ZigZag':
2379	dX = posZigZag(tau,delt5[:2],delt5[2:])
2380	else:
2381	dX = posBlock(tau,delt5[:2],delt5[2:])
2382	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2383	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2384	XSC = np.ones(6,dtype='i')
2385	dXTP = []
2386	for i in SdIndx:
2387	sop = Sop[i]
2388	ssop = SSop[i]
2389	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2390	xsc = np.ones(6,dtype='i')
2391	if 'Fourier' in waveType:
2392	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2393	elif waveType == 'Sawtooth':
2394	dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2395	elif waveType == 'ZigZag':
2396	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2397	elif waveType == 'Block':
2398	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,np.newaxis,np.newaxis]
2399	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2400	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2401	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2402	dXTP.append(dXT)
2403	if waveType == 'Fourier':
2404	for i in range(3):
2405	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2406	xsc[i] = 0
2407	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2408	xsc[i+3] = 0
2409	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2410	xsc[3:6] = 0
2411	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2412	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2413	if dT:
2414	if '(x)' in siteSym:
2415	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2416	if 'm' in siteSym and len(SdIndx) == 1:
2417	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2418	elif '(y)' in siteSym:
2419	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2420	if 'm' in siteSym and len(SdIndx) == 1:
2421	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2422	elif '(z)' in siteSym:
2423	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2424	if 'm' in siteSym and len(SdIndx) == 1:
2425	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2426	else:
2427	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2428	if '4/mmm' in laue:
2429	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2430	if '(xy)' in siteSym:
2431	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2432	if dT:
2433	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2434	else:
2435	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2436	if '(xy)' in siteSym or '(+-0)' in siteSym:
2437	mul = 1
2438	if '(+-0)' in siteSym:
2439	mul = -1
2440	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2441	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2442	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2443	xsc[3:5] = 0
2444	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2445	CSI[0][:2] = [[12,0,0],[12,0,0]]
2446	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2447	xsc[:2] = 0
2448	XSC &= xsc
2449	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2450	if waveType == 'Fourier':
2451	n = -1
2452	if debug: print (XSC)
2453	for i,X in enumerate(XSC):
2454	if X:
2455	n += 1
2456	CSI[0][i][0] = n+1
2457	CSI[1][i][0] = 1.0
2458
2459	return list(CSI),dX,dXTP
2460
2461	def DoUij():
2462	dU,dUTP = None,None
2463	if siteSym == '1':
2464	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2465	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2466	elif siteSym == '-1':
2467	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2468	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2469	else:
2470	tau = np.linspace(0,1,49,True)
2471	delt12 = np.eye(12)*0.0001
2472	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2473	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2474	USC = np.ones(12,dtype='i')
2475	dUTP = []
2476	dtau = 0.
2477	for i in SdIndx:
2478	sop = Sop[i]
2479	ssop = SSop[i]
2480	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2481	usc = np.ones(12,dtype='i')
2482	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2483	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2484	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2485	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2486	dUT = dUT[:,:,np.argsort(tauT)]
2487	dUT[:,:6,:] =(ssdetsdet)
2488	dUTP.append(dUT)
2489	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2490	if dT:
2491	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2492	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2493	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2494	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2495	else:
2496	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2497	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2498	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2499	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2500	if 'mm2(x)' in siteSym and dT:
2501	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2502	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2503	elif '(xy)' in siteSym and dT:
2504	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2505	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2506	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2507	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2508	elif '(x)' in siteSym and dT:
2509	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2510	elif '(y)' in siteSym and dT:
2511	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2512	elif '(z)' in siteSym and dT:
2513	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2514	for i in range(6):
2515	if not USC[i]:
2516	CSI[0][i] = [0,0,0]
2517	CSI[1][i] = [0.,0.,0.]
2518	CSI[0][i+6] = [0,0,0]
2519	CSI[1][i+6] = [0.,0.,0.]
2520	else:
2521	for i in range(6):
2522	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2523	usc[i] = 0
2524	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2525	usc[i+6] = 0
2526	if np.any(dUT[1,0,:]):
2527	if '4/m' in siteSym:
2528	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2529	if ssop[1][3]:
2530	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2531	usc[9] = 1
2532	else:
2533	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2534	usc[9] = 0
2535	elif '4' in siteSym:
2536	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2537	CSI[0][:2] = [[11,0,0],[11,0,0]]
2538	if ssop[1][3]:
2539	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2540	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2541	usc[2] = 0
2542	usc[8] = 0
2543	usc[3] = 1
2544	usc[9] = 1
2545	else:
2546	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2547	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2548	usc[2] = 1
2549	usc[8] = 1
2550	usc[3] = 0
2551	usc[9] = 0
2552	elif 'xy' in siteSym or '+-0' in siteSym:
2553	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2554	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2555	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2556	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2557	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2558	usc[4:6] = 0
2559	usc[6:8] = 0
2560
2561	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2562	USC &= usc
2563	if debug: print (USC)
2564	if not np.any(dtau%.5):
2565	n = -1
2566	for i,U in enumerate(USC):
2567	if U:
2568	n += 1
2569	CSI[0][i][0] = n+1
2570	CSI[1][i][0] = 1.0
2571
2572	return list(CSI),dU,dUTP
2573
2574	def DoMag():
2575	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2576	if siteSym == '1':
2577	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2578	elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2579	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2580	else:
2581	delt6 = np.eye(6)*0.001
2582	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,np.newaxis,np.newaxis]
2583	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2584	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2585	MSC = np.ones(6,dtype='i')
2586	dMTP = []
2587	for i in SdIndx:
2588	sop = Sop[i]
2589	ssop = SSop[i]
2590	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2591	msc = np.ones(6,dtype='i')
2592	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,np.newaxis,np.newaxis]
2593	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2594	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2595	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2596	dMTP.append(dMT)
2597	for i in range(3):
2598	if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2599	msc[i] = 0
2600	if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2601	msc[i+3] = 0
2602	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2603	msc[3:6] = 0
2604	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2605	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2606	if dT:
2607	if '(x)' in siteSym:
2608	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2609	if 'm' in siteSym and len(SdIndx) == 1:
2610	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2611	elif '(y)' in siteSym:
2612	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2613	if 'm' in siteSym and len(SdIndx) == 1:
2614	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2615	elif '(z)' in siteSym:
2616	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2617	if 'm' in siteSym and len(SdIndx) == 1:
2618	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2619	else:
2620	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2621	if '4/mmm' in laue:
2622	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2623	if '(xy)' in siteSym:
2624	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2625	if dT:
2626	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2627	else:
2628	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2629	if '(xy)' in siteSym or '(+-0)' in siteSym:
2630	mul = 1
2631	if '(+-0)' in siteSym:
2632	mul = -1
2633	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2634	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2635	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2636	msc[3:5] = 0
2637	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2638	CSI[0][:2] = [[12,0,0],[12,0,0]]
2639	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2640	msc[:2] = 0
2641	MSC &= msc
2642	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2643	n = -1
2644	if debug: print (MSC)
2645	for i,M in enumerate(MSC):
2646	if M:
2647	n += 1
2648	CSI[0][i][0] = n+1
2649	CSI[1][i][0] = 1.0
2650
2651	return CSI,None,None
2652
2653	if debug: print ('super space group: '+SSGData['SSpGrp'])
2654	xyz = np.array(XYZ)%1.
2655	SGOps = copy.deepcopy(SGData['SGOps'])
2656	laue = SGData['SGLaue']
2657	siteSym = SytSym(XYZ,SGData)[0].strip()
2658	if debug: print ('siteSym: '+siteSym)
2659	SSGOps = copy.deepcopy(SSGData['SSGOps'])
2660	#expand ops to include inversions if any
2661	if SGData['SGInv']:
2662	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2663	SGOps.append([-op[0],-op[1]%1.])
2664	SSGOps.append([-sop[0],-sop[1]%1.])
2665	#build set of sym ops around special position
2666	SSop = []
2667	Sop = []
2668	Sdtau = []
2669	for iop,Op in enumerate(SGOps):
2670	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2671	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2672	SSop.append(SSGOps[iop])
2673	Sop.append(SGOps[iop])
2674	ssopinv = nl.inv(SSGOps[iop][0])
2675	mst = ssopinv[3][:3]
2676	epsinv = ssopinv[3][3]
2677	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
2678	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
2679	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2680	#setup displacement arrays
2681	tau = np.linspace(-1,1,49,True)
2682	#make modulation arrays - one parameter at a time
2683	if Stype == 'Sfrac':
2684	CSI,dF,dFTP = DoFrac()
2685	elif Stype == 'Spos':
2686	CSI,dF,dFTP = DoXYZ()
2687	CSI[0] = orderParms(CSI[0])
2688	elif Stype == 'Sadp':
2689	CSI,dF,dFTP = DoUij()
2690	CSI[0] = orderParms(CSI[0])
2691	elif Stype == 'Smag':
2692	CSI,dF,dFTP = DoMag()
2693
2694	if debug:
2695	return CSI,dF,dFTP
2696	else:
2697	return CSI
2698
2699	def MustrainNames(SGData):
2700	'Needs a doc string'
2701	laue = SGData['SGLaue']
2702	uniq = SGData['SGUniq']
2703	if laue in ['m3','m3m']:
2704	return ['S400','S220']
2705	elif laue in ['6/m','6/mmm','3m1']:
2706	return ['S400','S004','S202']
2707	elif laue in ['31m','3']:
2708	return ['S400','S004','S202','S211']
2709	elif laue in ['3R','3mR']:
2710	return ['S400','S220','S310','S211']
2711	elif laue in ['4/m','4/mmm']:
2712	return ['S400','S004','S220','S022']
2713	elif laue in ['mmm']:
2714	return ['S400','S040','S004','S220','S202','S022']
2715	elif laue in ['2/m']:
2716	SHKL = ['S400','S040','S004','S220','S202','S022']
2717	if uniq == 'a':
2718	SHKL += ['S013','S031','S211']
2719	elif uniq == 'b':
2720	SHKL += ['S301','S103','S121']
2721	elif uniq == 'c':
2722	SHKL += ['S130','S310','S112']
2723	return SHKL
2724	else:
2725	SHKL = ['S400','S040','S004','S220','S202','S022']
2726	SHKL += ['S310','S103','S031','S130','S301','S013']
2727	SHKL += ['S211','S121','S112']
2728	return SHKL
2729
2730	def HStrainVals(HSvals,SGData):
2731	laue = SGData['SGLaue']
2732	uniq = SGData['SGUniq']
2733	DIJ = np.zeros(6)
2734	if laue in ['m3','m3m']:
2735	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2736	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2737	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2738	elif laue in ['3R','3mR']:
2739	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2740	elif laue in ['4/m','4/mmm']:
2741	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2742	elif laue in ['mmm']:
2743	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2744	elif laue in ['2/m']:
2745	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
2746	if uniq == 'a':
2747	DIJ[5] = HSvals[3]
2748	elif uniq == 'b':
2749	DIJ[4] = HSvals[3]
2750	elif uniq == 'c':
2751	DIJ[3] = HSvals[3]
2752	else:
2753	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
2754	return DIJ
2755
2756	def HStrainNames(SGData):
2757	'Needs a doc string'
2758	laue = SGData['SGLaue']
2759	uniq = SGData['SGUniq']
2760	if laue in ['m3','m3m']:
2761	return ['D11','eA'] #add cubic strain term
2762	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2763	return ['D11','D33']
2764	elif laue in ['3R','3mR']:
2765	return ['D11','D12']
2766	elif laue in ['4/m','4/mmm']:
2767	return ['D11','D33']
2768	elif laue in ['mmm']:
2769	return ['D11','D22','D33']
2770	elif laue in ['2/m']:
2771	Dij = ['D11','D22','D33']
2772	if uniq == 'a':
2773	Dij += ['D23']
2774	elif uniq == 'b':
2775	Dij += ['D13']
2776	elif uniq == 'c':
2777	Dij += ['D12']
2778	return Dij
2779	else:
2780	Dij = ['D11','D22','D33','D12','D13','D23']
2781	return Dij
2782
2783	def MustrainCoeff(HKL,SGData):
2784	'Needs a doc string'
2785	#NB: order of terms is the same as returned by MustrainNames
2786	laue = SGData['SGLaue']
2787	uniq = SGData['SGUniq']
2788	h,k,l = HKL
2789	Strm = []
2790	if laue in ['m3','m3m']:
2791	Strm.append(h4+k4+l**4)
2792	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
2793	elif laue in ['6/m','6/mmm','3m1']:
2794	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
2795	Strm.append(l**4)
2796	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
2797	elif laue in ['31m','3']:
2798	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
2799	Strm.append(l**4)
2800	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
2801	Strm.append(4.0hkl(h+k))
2802	elif laue in ['3R','3mR']:
2803	Strm.append(h4+k4+l**4)
2804	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
2805	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
2806	Strm.append(4.0(klh2+hlk2+hkl*2))
2807	elif laue in ['4/m','4/mmm']:
2808	Strm.append(h4+k4)
2809	Strm.append(l**4)
2810	Strm.append(3.0(hk)**2)
2811	Strm.append(3.0((hl)*2+(kl)**2))
2812	elif laue in ['mmm']:
2813	Strm.append(h**4)
2814	Strm.append(k**4)
2815	Strm.append(l**4)
2816	Strm.append(3.0(hk)**2)
2817	Strm.append(3.0(hl)**2)
2818	Strm.append(3.0(kl)**2)
2819	elif laue in ['2/m']:
2820	Strm.append(h**4)
2821	Strm.append(k**4)
2822	Strm.append(l**4)
2823	Strm.append(3.0(hk)**2)
2824	Strm.append(3.0(hl)**2)
2825	Strm.append(3.0(kl)**2)
2826	if uniq == 'a':
2827	Strm.append(2.0kl**3)
2828	Strm.append(2.0lk**3)
2829	Strm.append(4.0klh*2)
2830	elif uniq == 'b':
2831	Strm.append(2.0lh**3)
2832	Strm.append(2.0hl**3)
2833	Strm.append(4.0hlk*2)
2834	elif uniq == 'c':
2835	Strm.append(2.0hk**3)
2836	Strm.append(2.0kh**3)
2837	Strm.append(4.0hkl*2)
2838	else:
2839	Strm.append(h**4)
2840	Strm.append(k**4)
2841	Strm.append(l**4)
2842	Strm.append(3.0(hk)**2)
2843	Strm.append(3.0(hl)**2)
2844	Strm.append(3.0(kl)**2)
2845	Strm.append(2.0kh**3)
2846	Strm.append(2.0hl**3)
2847	Strm.append(2.0lk**3)
2848	Strm.append(2.0hk**3)
2849	Strm.append(2.0lh**3)
2850	Strm.append(2.0kl**3)
2851	Strm.append(4.0klh*2)
2852	Strm.append(4.0hlk*2)
2853	Strm.append(4.0khl*2)
2854	return Strm
2855
2856	def MuShklMean(SGData,Amat,Shkl):
2857
2858	def genMustrain(xyz,Shkl):
2859	uvw = np.inner(Amat.T,xyz)
2860	Strm = np.array(MustrainCoeff(uvw,SGData))
2861	Sum = np.sum(np.multiply(Shkl,Strm))
2862	Sum = np.where(Sum > 0.01,Sum,0.01)
2863	Sum = np.sqrt(Sum)
2864	return Sum*xyz
2865
2866	PHI = np.linspace(0.,360.,30,True)
2867	PSI = np.linspace(0.,180.,30,True)
2868	X = np.outer(npcosd(PHI),npsind(PSI))
2869	Y = np.outer(npsind(PHI),npsind(PSI))
2870	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2871	XYZ = np.dstack((X,Y,Z))
2872	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
2873	return np.sqrt(np.sum(XYZ**2)/900.)
2874
2875	def Muiso2Shkl(muiso,SGData,cell):
2876	"this is to convert isotropic mustrain to generalized Shkls"
2877	import GSASIIlattice as G2lat
2878	A = G2lat.cell2AB(cell)[0]
2879
2880	def minMus(Shkl,muiso,H,SGData,A):
2881	U = np.inner(A.T,H)
2882	S = np.array(MustrainCoeff(U,SGData))
2883	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
2884	rad = np.sqrt(np.sum((Sum[:,np.newaxis]H)*2,axis=1))
2885	return (muiso-rad)**2
2886
2887	laue = SGData['SGLaue']
2888	PHI = np.linspace(0.,360.,60,True)
2889	PSI = np.linspace(0.,180.,60,True)
2890	X = np.outer(npsind(PHI),npsind(PSI))
2891	Y = np.outer(npcosd(PHI),npsind(PSI))
2892	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2893	HKL = np.dstack((X,Y,Z))
2894	if laue in ['m3','m3m']:
2895	S0 = [1000.,1000.]
2896	elif laue in ['6/m','6/mmm','3m1']:
2897	S0 = [1000.,1000.,1000.]
2898	elif laue in ['31m','3']:
2899	S0 = [1000.,1000.,1000.,1000.]
2900	elif laue in ['3R','3mR']:
2901	S0 = [1000.,1000.,1000.,1000.]
2902	elif laue in ['4/m','4/mmm']:
2903	S0 = [1000.,1000.,1000.,1000.]
2904	elif laue in ['mmm']:
2905	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
2906	elif laue in ['2/m']:
2907	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
2908	else:
2909	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
2910	1000.,1000.,0.,0.,0.]
2911	S0 = np.array(S0)
2912	HKL = np.reshape(HKL,(-1,3))
2913	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
2914	return result[0]
2915
2916	def PackRot(SGOps):
2917	IRT = []
2918	for ops in SGOps:
2919	M = ops[0]
2920	irt = 0
2921	for j in range(2,-1,-1):
2922	for k in range(2,-1,-1):
2923	irt *= 3
2924	irt += M[k][j]
2925	IRT.append(int(irt))
2926	return IRT
2927
2928	def SytSym(XYZ,SGData):
2929	'''
2930	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
2931
2932	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
2933	:param SGData: from SpcGroup
2934	:Returns: a four element tuple:
2935
2936	* The 1st element is a code for the site symmetry (see GetKNsym)
2937	* The 2nd element is the site multiplicity
2938	* Ndup number of overlapping operators
2939	* dupDir Dict - dictionary of overlapping operators
2940
2941	'''
2942	Mult = 1
2943	Isym = 0
2944	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2945	Isym = 1073741824
2946	Jdup = 0
2947	Ndup = 0
2948	dupDir = {}
2949	inv = SGData['SGInv']+1
2950	icen = SGData['SGCen']
2951	if SGData['SGFixed']: #already in list of operators
2952	inv = 1
2953	Xeqv = list(GenAtom(XYZ,SGData,True))
2954	# for xeqv in Xeqv: print(xeqv)
2955	IRT = PackRot(SGData['SGOps'])
2956	L = -1
2957	for ic,cen in enumerate(icen):
2958	for invers in range(int(inv)):
2959	for io,ops in enumerate(SGData['SGOps']):
2960	irtx = (1-2invers)IRT[io]
2961	L += 1
2962	if not Xeqv[L][1]:
2963	Ndup = io
2964	Jdup += 1
2965	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
2966	if jx < 39:
2967	px = GetOprName(str(irtx))
2968	if Xeqv[L][-1] < 0:
2969	if '(' in px:
2970	px = px.split('(')
2971	px[0] += "'"
2972	px = '('.join(px)
2973	else:
2974	px += "'"
2975	dupDir[px] = L
2976	Isym += 2**(jx-1)
2977	if Isym == 1073741824: Isym = 0
2978	Mult = len(SGData['SGOps'])len(SGData['SGCen'])inv//Jdup
2979
2980	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
2981
2982	def MagSytSym(SytSym,dupDir,SGData):
2983	'''
2984	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
2985	'''
2986	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
2987	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
2988	# print('dupDir',dupDir)
2989	SplitSytSym = SytSym.split('(')
2990	if SGData['SGGray']:
2991	return SytSym+"1'"
2992	if SytSym == '1': #genersl position
2993	return SytSym
2994	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
2995	try:
2996	MagSytSym = SGData['MagSpGrp'].split()[1]
2997	except IndexError:
2998	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
2999	if len(SplitSytSym) > 1:
3000	MagSytSym += '('+SplitSytSym[1]
3001	return MagSytSym
3002	if len(dupDir) == 1:
3003	return dupDir.keys()[0]
3004
3005
3006	if '2/m' in SytSym: #done I think; last 2wo might be not needed
3007	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3008	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3009	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3010	'(110)':['2(110)','m(110)']}
3011
3012	elif '4/mmm' in SytSym:
3013	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
3014	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
3015	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
3016	elif '4mm' in SytSym:
3017	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3018	elif '422' in SytSym:
3019	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3020	elif '-4m2' in SytSym:
3021	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3022	elif '-42m' in SytSym:
3023	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3024	elif '-4' in SytSym:
3025	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3026	elif '4' in SytSym:
3027	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3028
3029	elif '222' in SytSym:
3030	ops = {'':['2(x)','2(y)','2(z)'],
3031	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3032	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3033	'(110)':['2(z)','2(110)','2(+-0)',],}
3034	elif 'mm2' in SytSym:
3035	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3036	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3037	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3038	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3039	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3040	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3041	'(210)':['m(z)','m(010)','2(210)'],
3042	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3043	'(110)':['m(z)','m(+-0)','2(110)',],}
3044	elif 'mmm' in SytSym:
3045	ops = {'':['m(x)','m(y)','m(z)'],
3046	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3047	'(110)':['m(z)','m(110)','m(+-0)',],
3048	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3049
3050	elif '32' in SytSym:
3051	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3052	elif '23' in SytSym:
3053	ops = {'':['2(x)','3(111)']}
3054	elif 'm3' in SytSym:
3055	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3056	elif '3m' in SytSym:
3057	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3058	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3059	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3060
3061	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]: #not simple cases
3062	MagSytSym = SytSym
3063	if "-1'" in dupDir:
3064	if '-6' in SytSym:
3065	MagSytSym = MagSytSym.replace('-6',"-6'")
3066	elif '-3m' in SytSym:
3067	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3068	elif '-3' in SytSym:
3069	MagSytSym = MagSytSym.replace('-3',"-3'")
3070	elif '-6m2' in SytSym:
3071	if "m'(110)" in dupDir:
3072	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3073	elif '6/m' in SytSym:
3074	if "m'(z)" in dupDir:
3075	MagSytSym = "6'/m'"
3076	elif '6mm' in SytSym:
3077	if "m'(110)" in dupDir:
3078	MagSytSym = "6'm'm"
3079	return MagSytSym
3080	try:
3081	axis = '('+SytSym.split('(')[1]
3082	except IndexError:
3083	axis = ''
3084	MagSytSym = ''
3085	for m in ops[axis]:
3086	if m in dupDir:
3087	MagSytSym += m.split('(')[0]
3088	else:
3089	MagSytSym += m.split('(')[0]+"'"
3090	if '2/m' in SytSym and '2' in m:
3091	MagSytSym += '/'
3092	if '-3/m' in SytSym:
3093	MagSytSym = '-'+MagSytSym
3094
3095	MagSytSym += axis
3096	# some exceptions & special rules
3097	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3098	return MagSytSym
3099
3100	# if len(GenSym) == 3:
3101	# if SGSpin[1] < 0:
3102	# if 'mm2' in SytSym:
3103	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3104	# else: #bad rule for I41/a
3105	# MagSytSym = SplitSytSym[0]+"'"
3106	# if len(SplitSytSym) > 1:
3107	# MagSytSym += '('+SplitSytSym[1]
3108	# else:
3109	# MagSytSym = SytSym
3110	# if len(SplitSytSym) >1:
3111	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3112	# MagSytSym = MagSytSym.replace('-4',"-4'")
3113	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3114	# MagSytSym = MagSytSym.replace('-6',"-6'")
3115	# return MagSytSym
3116	#
3117	return SytSym
3118
3119	def ElemPosition(SGData):
3120	''' Under development.
3121	Object here is to return a list of symmetry element types and locations suitable
3122	for say drawing them.
3123	So far I have the element type... getting all possible locations without lookup may be impossible!
3124	'''
3125	Inv = SGData['SGInv']
3126	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3127	# get operators & expand if centrosymmetric
3128	Ops = SGData['SGOps']
3129	opM = np.array([op[0].T for op in Ops])
3130	opT = np.array([op[1] for op in Ops])
3131	if Inv:
3132	opM = np.concatenate((opM,-opM))
3133	opT = np.concatenate((opT,-opT))
3134	opMT = list(zip(opM,opT))
3135	for M,T in opMT[1:]: #skip I
3136	Dt = int(nl.det(M))
3137	Tr = int(np.trace(M))
3138	Dt = -(Dt-1)//2
3139	Es = eleSym[Tr][Dt]
3140	if Dt: #rotation-inversion
3141	I = np.eye(3)
3142	if Tr == 1: #mirrors/glides
3143	if np.any(T): #glide
3144	M2 = np.inner(M,M)
3145	MT = np.inner(M,T)+T
3146	print ('glide',Es,MT)
3147	print (M2)
3148	else: #mirror
3149	print ('mirror',Es,T)
3150	print (I-M)
3151	X = [-1,-1,-1]
3152	elif Tr == -3: # pure inversion
3153	X = np.inner(nl.inv(I-M),T)
3154	print ('inversion',Es,X)
3155	else: #other rotation-inversion
3156	M2 = np.inner(M,M)
3157	MT = np.inner(M,T)+T
3158	print ('rot-inv',Es,MT)
3159	print (M2)
3160	X = [-1,-1,-1]
3161	else: #rotations
3162	print ('rotation',Es)
3163	X = [-1,-1,-1]
3164	#SymElements.append([Es,X])
3165
3166	return #SymElements
3167
3168	def ApplyStringOps(A,SGData,X,Uij=[]):
3169	'Needs a doc string'
3170	SGOps = SGData['SGOps']
3171	SGCen = SGData['SGCen']
3172	Ax = A.split('+')
3173	Ax[0] = int(Ax[0])
3174	iC = 0
3175	if Ax[0] < 0:
3176	iC = 1
3177	Ax[0] = abs(Ax[0])
3178	nA = Ax[0]%100-1
3179	cA = Ax[0]//100
3180	Cen = SGCen[cA]
3181	M,T = SGOps[nA]
3182	if len(Ax)>1:
3183	cellA = Ax[1].split(',')
3184	cellA = np.array([int(a) for a in cellA])
3185	else:
3186	cellA = np.zeros(3)
3187	newX = Cen+(1-2iC)(np.inner(M,X).T+T)+cellA
3188	if len(Uij):
3189	U = Uij2U(Uij)
3190	U = np.inner(M,np.inner(U,M).T)
3191	newUij = U2Uij(U)
3192	return [newX,newUij]
3193	else:
3194	return newX
3195
3196	def ApplyStringOpsMom(A,SGData,Mom):
3197	'Needs a doc string'
3198	SGOps = SGData['SGOps']
3199	Ax = A.split('+')
3200	Ax[0] = int(Ax[0])
3201	Ax[0] = abs(Ax[0])
3202	nA = Ax[0]%100-1
3203	M,T = SGOps[nA]
3204	if SGData['SGGray']:
3205	newMom = -(np.inner(Mom,M).T)*nl.det(M)
3206	else:
3207	newMom = -(np.inner(Mom,M).T)SGData['MagMom'][nA-1]nl.det(M)
3208	return newMom
3209
3210	def StringOpsProd(A,B,SGData):
3211	"""
3212	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3213	where '-' indicates inversion, c(>0) is the cell centering operator,
3214	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3215	Should return resultant string - C. SGData - dictionary using entries:
3216
3217	* 'SGCen': cell centering vectors [0,0,0] at least
3218	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3219
3220	"""
3221	SGOps = SGData['SGOps']
3222	SGCen = SGData['SGCen']
3223	#1st split out the cell translation part & work on the operator parts
3224	Ax = A.split('+'); Bx = B.split('+')
3225	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3226	iC = 0
3227	if Ax[0]*Bx[0] < 0:
3228	iC = 1
3229	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3230	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3231	cA = Ax[0]//100; cB = Bx[0]//100
3232	Cen = (SGCen[cA]+SGCen[cB])%1.0
3233	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3234	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3235	Mc = np.inner(Ma,Mb.T)
3236	# print Ma,Mb,Mc
3237	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3238	# print Ta,Tb,Tc
3239	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3240	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3241	#now the cell translation part
3242	if len(Ax)>1:
3243	cellA = Ax[1].split(',')
3244	cellA = [int(a) for a in cellA]
3245	else:
3246	cellA = [0,0,0]
3247	if len(Bx)>1:
3248	cellB = Bx[1].split(',')
3249	cellB = [int(b) for b in cellB]
3250	else:
3251	cellB = [0,0,0]
3252	cellC = np.add(cellA,cellB)
3253	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3254	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3255	return C
3256
3257	def U2Uij(U):
3258	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3259	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3260
3261	def Uij2U(Uij):
3262	#returns the thermal motion tensor U from Uij as numpy array
3263	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3264
3265	def StandardizeSpcName(spcgroup):
3266	'''Accept a spacegroup name where spaces may have not been used
3267	in the names according to the GSAS convention (spaces between symmetry
3268	for each axis) and return the space group name as used in GSAS
3269	'''
3270	rspc = spcgroup.replace(' ','').upper()
3271	# deal with rhombohedral and hexagonal setting designations
3272	rhomb = ''
3273	if rspc[-1:] == 'R':
3274	rspc = rspc[:-1]
3275	rhomb = ' R'
3276	gray = ''
3277	if "1'" in rspc:
3278	gray = " 1'"
3279	rspc = rspc.replace("1'",'')
3280	elif rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3281	rspc = rspc[:-1]
3282	else:
3283	rspc = rspc.replace("'",'')
3284	# look for a match in the spacegroup lists
3285	for i in spglist.values():
3286	for spc in i:
3287	if rspc == spc.replace(' ','').upper():
3288	return spc+gray+rhomb
3289	# how about the post-2002 orthorhombic names?
3290	if rspc in sgequiv_2002_orthorhombic:
3291	return sgequiv_2002_orthorhombic[rspc]+gray
3292	else:
3293	# not found
3294	return ''
3295
3296	spgbyNum = []
3297	'''Space groups indexed by number'''
3298	spgbyNum = [None,
3299	'P 1','P -1', #1-2
3300	'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3301	'C 2/m','P 2/c','P 21/c','C 2/c', #3-15
3302	'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3303	'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3304	'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3305	'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3306	'C m m 2','C m c 21','C c c 2',
3307	'A m m 2','A b m 2','A m a 2','A b a 2',
3308	'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3309	'P m m m','P n n n','P c c m','P b a n',
3310	'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3311	'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3312	'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3313	'F m m m', 'F d d d',
3314	'I m m m','I b a m','I b c a','I m m a', #16-74
3315	'P 4','P 41','P 42','P 43',
3316	'I 4','I 41',
3317	'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3318	'I 4/m','I 41/a',
3319	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3320	'P 42 21 2','P 43 2 2','P 43 21 2',
3321	'I 4 2 2','I 41 2 2',
3322	'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3323	'P 42 m c','P 42 b c',
3324	'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3325	'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3326	'P -4 c 2','P -4 b 2','P -4 n 2',
3327	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3328	'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',
3329	'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',
3330	'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',
3331	'P 42/n c m',
3332	'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3333	'P 3','P 31','P 32','R 3','P -3','R -3',
3334	'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3335	'R 3 2',
3336	'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3337	'R 3 m','R 3 c',
3338	'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3339	'R -3 m','R -3 c', #75-167
3340	'P 6','P 61',
3341	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3342	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3343	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3344	'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
3345	'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3346	'F m -3','F d -3','I m -3',
3347	'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3348	'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3349	'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3350	'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3351	'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3352	'I m -3 m','I a -3 d',] #195-230
3353	spg2origins = {}
3354	''' A dictionary of all spacegroups that have 2nd settings; the value is the
3355	1st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3356	'''
3357	spg2origins = {
3358	'P n n n':[-.25,-.25,-.25],
3359	'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3360	'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3361	'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3362	'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3363	'F d d d':[-.125,-.125,-.125],
3364	'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3365	'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],
3366	'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],
3367	'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3368	'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3369	'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3370	'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3371	'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3372	spglist = {}
3373	'''A dictionary of space groups as ordered and named in the pre-2002 International
3374	Tables Volume A, except that spaces are used following the GSAS convention to
3375	separate the different crystallographic directions.
3376	Note that the symmetry codes here will recognize many non-standard space group
3377	symbols with different settings. They are ordered by Laue group
3378	'''
3379	spglist = {
3380	'P1' : ('P 1','P -1',), # 1-2
3381	'C1' : ('C 1','C -1',),
3382	'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3383	'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
3384	'C2/m':('C 2','C m','C c','C n',
3385	'C 2/m','C 2/c','C 2/n',),
3386	'A2/m':('A 2','A m','A a','A n',
3387	'A 2/m','A 2/a','A 2/n',),
3388	'Pmmm':('P 2 2 2',
3389	'P 2 2 21','P 21 2 2','P 2 21 2',
3390	'P 21 21 2','P 2 21 21','P 21 2 21',
3391	'P 21 21 21',
3392	'P m m 2','P 2 m m','P m 2 m',
3393	'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3394	'P c c 2','P 2 a a','P b 2 b',
3395	'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3396	'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3397	'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3398	'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3399	'P b a 2','P 2 c b','P c 2 a',
3400	'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3401	'P n n 2','P 2 n n','P n 2 n',
3402	'P m m m','P n n n',
3403	'P c c m','P m a a','P b m b',
3404	'P b a n','P n c b','P c n a',
3405	'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3406	'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3407	'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3408	'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3409	'P b a m','P m c b','P c m a',
3410	'P c c n','P n a a','P b n b',
3411	'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3412	'P n n m','P m n n','P n m n',
3413	'P m m n','P n m m','P m n m',
3414	'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3415	'P b c a','P c a b',
3416	'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3417	),
3418	'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3419	'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3420	'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3421	'C m c m','C c m m','C m c a','C c m b',
3422	'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3423	'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3424	'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3425	'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3426	'A m m a','A m a m','A b m a','A c a m',
3427	'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3428	'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3429	'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3430	'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3431	'B b m m','B m m b','B b c m','B m a b',
3432	'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3433	'Immm':('I 2 2 2','I 21 21 21',
3434	'I m m 2','I m 2 m','I 2 m m',
3435	'I b a 2','I 2 c b','I c 2 a',
3436	'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3437	'I m m m','I b a m','I m c b','I c m a',
3438	'I b c a','I c a b',
3439	'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3440	'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3441	'F m m 2','F m 2 m','F 2 m m',
3442	'F d d 2','F d 2 d','F 2 d d',),
3443	'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3444	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3445	'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3446	'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3447	'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3448	'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',
3449	'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',
3450	'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',
3451	'P 42/n c m',),
3452	'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',
3453	'I 4 c m','I 41 m d','I 41 c d',
3454	'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',
3455	'I 41/a m d','I 41/a c d'),
3456	'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3457	'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3458	'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3459	'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3460	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3461	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3462	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3463	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3464	'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',
3465	'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3466	'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3467	'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3468	'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',
3469	'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'),
3470	'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3471	'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3472	'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3473	'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3474	}
3475	sgequiv_2002_orthorhombic = {}
3476	''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3477	along with the pre-2002 name. The e designates a double glide-plane
3478	'''
3479	sgequiv_2002_orthorhombic = {
3480	'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3481	'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3482	'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3483	'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3484	'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3485	'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3486	'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3487	'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3488
3489	#'A few non-standard space groups for test use'
3490	nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h',
3491	'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3492	'R 3 c r','R -3 c r','R -3 m r',),
3493
3494	#Use the space groups types in this order to list the symbols in the
3495	#order they are listed in the International Tables, vol. A'''
3496	symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal',
3497	'trigonal', 'hexagonal', 'cubic')
3498
3499	# self-test materials follow. Requires files in directory testinp
3500	selftestlist = []
3501	'''Defines a list of self-tests'''
3502	selftestquiet = True
3503	def _ReportTest():
3504	'Report name and doc string of current routine when ``selftestquiet`` is False'
3505	if not selftestquiet:
3506	import inspect
3507	caller = inspect.stack()[1][3]
3508	doc = eval(caller).__doc__
3509	if doc is not None:
3510	print('testing '+__file__+' with '+caller+' ('+doc+')')
3511	else:
3512	print('testing '+__file__()+" with "+caller)
3513	def test0():
3514	'''self-test #0: exercise MoveToUnitCell'''
3515	_ReportTest()
3516	msg = "MoveToUnitCell failed"
3517	assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3518	assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3519	assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3520	assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3521	selftestlist.append(test0)
3522
3523	def test1():
3524	'''self-test #1: SpcGroup against previous results'''
3525	#'''self-test #1: SpcGroup and SGPrint against previous results'''
3526	_ReportTest()
3527	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3528	if ospath.exists(testdir):
3529	if testdir not in sys.path: sys.path.insert(0,testdir)
3530	import spctestinp
3531	def CompareSpcGroup(spc, referr, refdict, reflist):
3532	'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3533	# if an error is reported, the dictionary can be ignored
3534	msg0 = "CompareSpcGroup failed on space group %s" % spc
3535	result = SpcGroup(spc)
3536	if result[0] == referr and referr > 0: return True
3537	# #print result[1]['SpGrp']
3538	#msg = msg0 + " in list lengths"
3539	#assert len(keys) == len(refdict.keys()), msg
3540	for key in refdict.keys():
3541