Context Navigation

source: trunk/GSASIIspc.py @ 3795

Last change on this file since 3795 was 3795, checked in by vondreele, 4 years ago
fix space group & incommensurate space group displays for mag. incomensurate cases disable cell refine in Unit Cell List if no peaks another shot at incommensurate mag. structure factor calc.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 173.3 KB

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