Context Navigation

source: trunk/GSASIIspc.py @ 3768

Last change on this file since 3768 was 3768, checked in by vondreele, 4 years ago
name fix for single crystal 1D plots remove ZigZag? & Block from modulation derivatives - now done numerically fix ZigZag? modulation function & proper limits for Tmin & Tmax in GUI
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 173.6 KB

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