Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: trunk/GSASIIspc.py

Last change on this file was 5662, checked in by toby, 4 days ago
docs formatting in g2spc
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 180.3 KB

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