Context Navigation

source: trunk/GSASIIspc.py @ 3617

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