Context Navigation

source: trunk/GSASIIspc.py @ 3297

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