Context Navigation

source: trunk/GSASIIspc.py @ 5119

Last change on this file since 5119 was 5119, checked in by toby, 2 years ago
fix (partial?) problems with reading short HM names in fulHM routine
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 176.6 KB

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