Context Navigation

source: trunk/GSASIIspc.py @ 3753

Last change on this file since 3753 was 3753, checked in by vondreele, 5 years ago
fix run k-SUBGROUPSMAG bug
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 174.0 KB

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