Context Navigation

source: trunk/GSASIIspc.py @ 3561

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