Context Navigation

source: trunk/GSASIIspc.py @ 4213

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