Context Navigation

source: trunk/GSASIIspc.py @ 4019

Last change on this file since 4019 was 4019, checked in by vondreele, 4 years ago
put directions for delete restraints in status bar; fix bug in G2spc
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 173.5 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-06-07 12:48:58 +0000 (Fri, 07 Jun 2019) $
12	# $Author: vondreele $
13	# $Revision: 4019 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 4019 2019-06-07 12:48:58Z 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: 4019 $")
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	if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1665	Ncen = len(SGData['SGCen'])
1666	for icen in range(Ncen//2,Ncen):
1667	SSCen[icen,3] = 0.5
1668	SSGData['SSGCen'] = SSCen%1.
1669	SSGData['SSGOps'] = []
1670	for iop,op in enumerate(SGData['SGOps']):
1671	T = np.zeros(4)
1672	ssop = np.zeros((4,4))
1673	ssop[:3,:3] = op[0]
1674	T[:3] = op[1]
1675	SSGData['SSGOps'].append([ssop,T])
1676	E,Result = genSSGOps()
1677	if E:
1678	SSGData['SSGOps'] = Result
1679	if DEBUG:
1680	print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1681	for Op in Result:
1682	print (SSMT2text(Op).replace(' ',''))
1683	if SGData['SGInv']:
1684	for Op in Result:
1685	Op = [-Op[0],-Op[1]%1.]
1686	print (SSMT2text(Op).replace(' ',''))
1687	return None,SSGData
1688	else:
1689	return Result+'\nOperator conflict - incorrect superspace symbol',None
1690
1691	def SSChoice(SGData):
1692	'''
1693	Gets the unique set of possible super space groups for a given space group
1694	'''
1695	ptgpSS = {'1':['(abg)',],'-1':['(abg)',],
1696
1697	'2':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1698	'm':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1699	'2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1700
1701	'222':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1702	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1703	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1704	'mm2':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1705	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1706	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1707	'm2m':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1708	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1709	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1710	'2mm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1711	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1712	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1713	'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1714	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1715	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1716
1717	'4':['(00g)','(1/21/2g)'],'4mm':['(00g)','(1/21/2g)'],
1718	'4/m':['(00g)','(1/21/2g)'],
1719	'422':['(00g)','(1/21/2g)'],'-4m2':['(00g)','(1/21/2g)'],'-42m':['(00g)','(1/21/2g)'],
1720	'4/mmm':['(00g)','(1/21/2g)'],
1721
1722	'3':['(00g)','(1/31/3g)'],'-3':['(00g)','(1/31/3g)'],
1723	'32':['(00g)'],'3m':['(00g)'],'-3m':['(00g)'],
1724	'321':['(00g)'],'3m1':['(00g)'],'-3m1':['(00g)'],
1725	'312':['(00g)','(1/31/3g)'],'31m':['(00g)','(1/31/3g)'],'-31m':['(00g)','(1/31/3g)'],
1726
1727	'6':['(00g)',],'6/m':['(00g)',],'-62m':['(00g)',],'-6m2':['(00g)',],
1728	'622':['(00g)',],'6/mmm':['(00g)',],'6mm':['(00g)',],
1729
1730	'23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1731
1732	ptgpTS = {'1':['0',],'-1':['0',],
1733
1734	'2':['0','s'],'m':['0','s'],
1735	'2/m':['00','0s','ss','s0'],
1736
1737	'222':['000','s00','0s0','00s',],
1738	'mm2':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1739	'm2m':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1740	'2mm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1741	'mmm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1742
1743	'4':['0','q','s'],'4mm':['000','q00','s00','s0s','ss0','0ss','qq0','qqs'],
1744	'4/m':['00','s0'],'-4m2':['000','0s0','0q0'],'-42m':['000','00s'],
1745	'422':['000','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0'],
1746	'4/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1747
1748	'3':['0','t'],'-3':['0','t'],
1749	'32':['00','t0'],'3m':['00','0s'],'-3m':['00','0s'],
1750	'321':['000','t00'],'3m1':['000','0s0'],'-3m1':['000','0s0'],
1751	'312':['000','t00'],'31m':['000','00s'],'-31m':['000','00s'],
1752
1753	'6':['0','h','t','s'],
1754	'6/m':['00','s0'],'-62m':['000','00s'],'-6m2':['000','0s0'],
1755	'622':['000','h00','t00','s00',],'6mm':['000','ss0','s0s','0ss',],
1756	'6/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1757
1758	'23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1759
1760	ptgp = SGData['SGPtGrp']
1761	SSChoice = []
1762	for ax in ptgpSS[ptgp]:
1763	for sx in ptgpTS[ptgp]:
1764	SSChoice.append(ax+sx)
1765	if SGData['SGGray']: SSChoice[-1] += 's'
1766	ssChoice = []
1767	ssHash = []
1768	for item in SSChoice:
1769	E,SSG = SSpcGroup(SGData,item)
1770	if SSG:
1771	sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1772	if sshash not in ssHash:
1773	ssHash.append(sshash)
1774	ssChoice.append(item)
1775	return ssChoice
1776
1777	def splitSSsym(SSymbol):
1778	'''
1779	Splits supersymmetry symbol into two lists of strings
1780	'''
1781	mssym = SSymbol.replace(' ','').split(')')
1782	if len(mssym) > 1:
1783	modsym,gensym = mssym
1784	else:
1785	modsym = mssym[0]
1786	gensym = ''
1787	modsym = modsym.replace(',','')
1788	if "1'" in modsym:
1789	gensym = gensym[:-1]
1790	modsym = modsym.replace("1'",'')
1791	if gensym in ['0','00','000','0000']: #get rid of extraneous symbols
1792	gensym = ''
1793	nfrac = modsym.count('/')
1794	modsym = modsym.lstrip('(')
1795	if nfrac == 0:
1796	modsym = list(modsym)
1797	elif nfrac == 1:
1798	pos = modsym.find('/')
1799	if pos == 1:
1800	modsym = [modsym[:3],modsym[3],modsym[4]]
1801	elif pos == 2:
1802	modsym = [modsym[0],modsym[1:4],modsym[4]]
1803	else:
1804	modsym = [modsym[0],modsym[1],modsym[2:]]
1805	else:
1806	lpos = modsym.find('/')
1807	rpos = modsym.rfind('/')
1808	if lpos == 1 and rpos == 4:
1809	modsym = [modsym[:3],modsym[3:6],modsym[6]]
1810	elif lpos == 1 and rpos == 5:
1811	modsym = [modsym[:3],modsym[3],modsym[4:]]
1812	else:
1813	modsym = [modsym[0],modsym[1:4],modsym[4:]]
1814	gensym = list(gensym)
1815	return modsym,gensym
1816
1817	def SSGPrint(SGData,SSGData,AddInv=False):
1818	'''
1819	Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1820
1821	:param SGData: space group data structure as defined in SpcGroup above.
1822	:param SSGData: from :func:`SSpcGroup`
1823	:returns:
1824	SSGText - list of strings with the superspace group details
1825	SGTable - list of strings for each of the operations
1826	'''
1827	nCen = len(SSGData['SSGCen'])
1828	Mult = nCenlen(SSGData['SSGOps'])(int(SGData['SGInv'])+1)
1829	if SGData.get('SGFixed',False):
1830	Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])
1831	SSsymb = SSGData['SSpGrp']
1832	if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1833	SSsymb = SGData['BNSlattsym'][0]+SSsymb[1:]
1834	if SGData.get('SGGray',False):
1835	if SGData.get('SGFixed',False): Mult //= 2
1836	else:
1837	if "1'" in SSsymb: #leftover in nonmag phase in mcif file
1838	nCen //= 2
1839	Mult //= 2
1840	SSsymb = SSsymb.replace("1'",'')[:-1]
1841	SSGText = []
1842	SSGText.append(' Superspace Group: '+SSsymb)
1843	CentStr = 'centrosymmetric'
1844	if not SGData['SGInv']:
1845	CentStr = 'non'+CentStr
1846	if SGData['SGLatt'] in 'ABCIFR':
1847	SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1848	else:
1849	SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1850	SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1851	SGptGp = SGData['SGPtGrp']
1852	if SGData['SGGray']:
1853	SGptGp += "1'"
1854	SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1855	SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1856	SSGText.append(' Multiplicity of a general site is '+str(Mult))
1857	if SGData['SGUniq'] in ['a','b','c']:
1858	SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1859	if SGData['SGInv']:
1860	SSGText.append(' The inversion center is located at 0,0,0')
1861	if SGData['SGPolax']:
1862	SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1863	SSGText.append(' ')
1864	if len(SSGData['SSGCen']) > 1:
1865	SSGText.append(' The equivalent positions are:')
1866	SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'][:nCen])+')+\n')
1867	else:
1868	SSGText.append(' The equivalent positions are:\n')
1869	SSGTable = []
1870	for i,Opr in enumerate(SSGData['SSGOps']):
1871	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1872	if AddInv and SGData['SGInv']:
1873	for i,Opr in enumerate(SSGData['SSGOps']):
1874	IOpr = [-Opr[0],-Opr[1]]
1875	SSGTable.append('(%2d) %s'%(i+1+len(SSGData['SSGOps']),SSMT2text(IOpr)))
1876	return SSGText,SSGTable
1877
1878	def SSGModCheck(Vec,modSymb,newMod=True):
1879	''' Checks modulation vector compatibility with supersymmetry space group symbol.
1880	if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1881	'''
1882	Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1883	modQ = [Fracs[mod] for mod in modSymb]
1884	if newMod:
1885	newVec = Vec
1886	if not np.any(Vec):
1887	newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1888	return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)], \
1889	[True if mod in ['a','b','g'] else False for mod in modSymb]
1890	else:
1891	return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1892
1893	def SSMT2text(Opr):
1894	"From superspace group matrix/translation operator returns text version"
1895	XYZS = ('x','y','z','t') #Stokes, Campbell & van Smaalen notation
1896	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1897	Fld = ''
1898	M,T = Opr
1899	for j in range(4):
1900	IJ = ''
1901	for k in range(4):
1902	txt = str(int(round(M[j][k])))
1903	txt = txt.replace('1',XYZS[k]).replace('0','')
1904	if '2' in txt:
1905	txt += XYZS[k]
1906	if IJ and M[j][k] > 0:
1907	IJ += '+'+txt
1908	else:
1909	IJ += txt
1910	IK = int(round(T[j]*12))%12
1911	if IK:
1912	if not IJ:
1913	break
1914	if IJ[0] == '-':
1915	Fld += (TRA[IK]+IJ).rjust(8)
1916	else:
1917	Fld += (TRA[IK]+'+'+IJ).rjust(8)
1918	else:
1919	Fld += IJ.rjust(8)
1920	if j != 3: Fld += ', '
1921	return Fld
1922
1923	def SSLatt2text(SSGCen):
1924	"Lattice centering vectors to text"
1925	lattTxt = ''
1926	lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1927	for vec in SSGCen:
1928	lattTxt += ' '
1929	for item in vec:
1930	lattTxt += '%s,'%(lattDir[int(item*12)])
1931	lattTxt = lattTxt.rstrip(',')
1932	lattTxt += ';'
1933	lattTxt = lattTxt.rstrip(';').lstrip(' ')
1934	return lattTxt
1935
1936	def SSpaceGroup(SGSymbol,SSymbol):
1937	'''
1938	Print the output of SSpcGroup in a nicely formatted way.
1939
1940	:param SGSymbol: space group symbol with spaces between axial fields.
1941	:param SSymbol: superspace group symbol extension (string).
1942	:returns: nothing
1943	'''
1944
1945	E,A = SpcGroup(SGSymbol)
1946	if E > 0:
1947	print (SGErrors(E))
1948	return
1949	E,B = SSpcGroup(A,SSymbol)
1950	if E > 0:
1951	print (E)
1952	return
1953	for l in SSGPrint(B):
1954	print (l)
1955
1956	def SGProd(OpA,OpB):
1957	'''
1958	Form space group operator product. OpA & OpB are [M,V] pairs;
1959	both must be of same dimension (3 or 4). Returns [M,V] pair
1960	'''
1961	A,U = OpA
1962	B,V = OpB
1963	M = np.inner(B,A.T)
1964	W = np.inner(B,U)+V
1965	return M,W
1966
1967	def GetLittleGrpOps(SGData,vec):
1968	''' Find rotation part of operators that leave vec unchanged
1969
1970	:param SGData: space group data structure as defined in SpcGroup above.
1971	:param vec: a numpy array of fractional vector coordinates
1972	:returns: Little - list of operators [M,T] that form the little gropu
1973	'''
1974	Little = []
1975	Ops = SGData['SGOps'][:]
1976	if SGData['SGInv']:
1977	Ops += [[-M,-T] for [M,T] in Ops]
1978	for [M,T] in Ops:
1979	tvec = np.inner(M,vec)%1.
1980	if np.allclose(tvec,vec%1.):
1981	Little.append([M,T])
1982	return Little
1983
1984	def MoveToUnitCell(xyz):
1985	'''
1986	Translates a set of coordinates so that all values are >=0 and < 1
1987
1988	:param xyz: a list or numpy array of fractional coordinates
1989	:returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1990	'''
1991	XYZ = (np.array(xyz)+10.)%1.
1992	cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1993	return XYZ,cell
1994
1995	def Opposite(XYZ,toler=0.0002):
1996	'''
1997	Gives opposite corner, edge or face of unit cell for position within tolerance.
1998	Result may be just outside the cell within tolerance
1999
2000	:param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
2001	:param toler: unit cell fraction tolerance making opposite
2002	:returns:
2003	XYZ: dict of opposite positions; key=unit cell & always contains XYZ
2004	'''
2005	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]]
2006	TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
2007	perm = TB*perm3
2008	cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
2009	D = dict(zip(cperm,perm))
2010	new = {}
2011	for key in D:
2012	new[key] = np.array(D[key])+np.array(XYZ)
2013	return new
2014
2015	def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
2016	'''
2017	Generates the equivalent positions for a specified coordinate and space group
2018
2019	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
2020	:param SGData: from :func:`SpcGroup`
2021	:param All: True return all equivalent positions including duplicates;
2022	False return only unique positions
2023	:param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
2024	:param Move: True move generated atom positions to be inside cell
2025	False do not move atoms
2026	:return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
2027
2028	* [XYZEquiv] is list of equivalent positions (XYZ is first entry)
2029	* Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
2030	* is centering operator number (1-4) and - is for inversion
2031	Cell = unit cell translations needed to put new positions inside cell
2032	[UijEquiv] - equivalent Uij; absent if no Uij given
2033	* +1/-1 for spin inversion of operator - empty if not magnetic
2034
2035	'''
2036	XYZEquiv = []
2037	UijEquiv = []
2038	Idup = []
2039	Cell = []
2040	inv = int(SGData['SGInv']+1)
2041	icen = SGData['SGCen']
2042	if SGData.get('SGFixed',False):
2043	inv = 1
2044	SpnFlp = SGData.get('SpnFlp',[])
2045	spnflp = []
2046	X = np.array(XYZ)
2047	mj = 0
2048	cell0 = np.zeros(3,dtype=np.int32)
2049	if Move:
2050	X,cell0 = MoveToUnitCell(X)
2051	for ic,cen in enumerate(icen):
2052	C = np.array(cen)
2053	for invers in range(inv):
2054	for io,[M,T] in enumerate(SGData['SGOps']):
2055	idup = ((io+1)+100ic)(1-2*invers)
2056	XT = np.inner(M,X)+T
2057	if len(Uij):
2058	U = Uij2U(Uij)
2059	U = np.inner(M,np.inner(U,M).T)
2060	newUij = U2Uij(U)
2061	if invers:
2062	XT = -XT
2063	XT += C
2064	cell = np.zeros(3,dtype=np.int32)+cell0
2065	cellj = np.zeros(3,dtype=np.int32)
2066	if Move:
2067	newX,cellj = MoveToUnitCell(XT)
2068	else:
2069	newX = XT
2070	cell += cellj
2071	if All:
2072	if np.allclose(newX,X,atol=0.0002):
2073	idup = False
2074	else:
2075	if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
2076	idup = False
2077	if All or idup:
2078	XYZEquiv.append(newX)
2079	Idup.append(idup)
2080	Cell.append(cell)
2081	if len(Uij):
2082	UijEquiv.append(newUij)
2083	if len(SpnFlp):
2084	spnflp.append(SpnFlp[mj])
2085	else:
2086	spnflp.append(1)
2087	mj += 1
2088	if len(Uij):
2089	return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2090	else:
2091	return zip(XYZEquiv,Idup,Cell,spnflp)
2092
2093	def GenHKL(HKL,SGData):
2094	''' Generates all equivlent reflections including Friedel pairs
2095	:param HKL: [h,k,l] must be integral values
2096	:param SGData: space group data obtained from SpcGroup
2097	:returns: array Uniq: equivalent reflections
2098	'''
2099
2100	Ops = SGData['SGOps']
2101	OpM = np.array([op[0] for op in Ops])
2102	Uniq = np.inner(OpM,HKL)
2103	Uniq = list(Uniq)+list(-1*Uniq)
2104	return np.array(Uniq)
2105
2106	def GenHKLf(HKL,SGData):
2107	'''
2108	Uses old GSAS Fortran routine genhkl.for
2109
2110	:param HKL: [h,k,l] must be integral values for genhkl.for to work
2111	:param SGData: space group data obtained from SpcGroup
2112	:returns: iabsnt,mulp,Uniq,phi
2113
2114	* iabsnt = True if reflection is forbidden by symmetry
2115	* mulp = reflection multiplicity including Friedel pairs
2116	* Uniq = numpy array of equivalent hkl in descending order of h,k,l
2117	* phi = phase offset for each equivalent h,k,l
2118
2119	'''
2120	hklf = list(HKL)+[0,] #could be numpy array!
2121	Ops = SGData['SGOps']
2122	OpM = np.array([op[0] for op in Ops],order='F')
2123	OpT = np.array([op[1] for op in Ops])
2124	Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2125
2126	import pyspg
2127	Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2128	h,k,l,f = Uniq
2129	Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2130	phi = f[:Nuniq]
2131	return iabsnt,mulp,Uniq,phi
2132
2133	def checkSSLaue(HKL,SGData,SSGData):
2134	#Laue check here - Toss HKL if outside unique Laue part
2135	h,k,l,m = HKL
2136	if SGData['SGLaue'] == '2/m':
2137	if SGData['SGUniq'] == 'a':
2138	if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2139	return False
2140	elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2141	return False
2142	else:
2143	return True
2144	elif SGData['SGUniq'] == 'b':
2145	if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2146	return False
2147	elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2148	return False
2149	else:
2150	return True
2151	elif SGData['SGUniq'] == 'c':
2152	if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2153	return False
2154	elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2155	return False
2156	else:
2157	return True
2158	elif SGData['SGLaue'] == 'mmm':
2159	if 'a' in SSGData['modSymb']:
2160	if h == 0 and m < 0:
2161	return False
2162	else:
2163	return True
2164	elif 'b' in SSGData['modSymb']:
2165	if k == 0 and m < 0:
2166	return False
2167	else:
2168	return True
2169	elif 'g' in SSGData['modSymb']:
2170	if l == 0 and m < 0:
2171	return False
2172	else:
2173	return True
2174	else: #tetragonal, trigonal, hexagonal (& triclinic?)
2175	if l == 0 and m < 0:
2176	return False
2177	else:
2178	return True
2179
2180	def checkHKLextc(HKL,SGData):
2181	'''
2182	Checks if reflection extinct - does not check centering
2183
2184	:param HKL: [h,k,l]
2185	:param SGData: space group data obtained from SpcGroup
2186	:returns: True if extinct; False if allowed
2187
2188	'''
2189	Ops = SGData['SGOps']
2190	OpM = np.array([op[0] for op in Ops])
2191	OpT = np.array([op[1] for op in Ops])
2192	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2193	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2194	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2195	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
2196	if dhkl.any():
2197	continue
2198	else:
2199	if phkl%1.:
2200	return True
2201	return False
2202
2203	def checkMagextc(HKL,SGData):
2204	'''
2205	Checks if reflection magnetically extinct; does fullcheck (centering, too)
2206	uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2207
2208	:param HKL: [h,k,l]
2209	:param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2210	:returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2211
2212	'''
2213	Ops = SGData['SGOps']
2214	Ncen = len(SGData['SGCen'])
2215	OpM = np.array([op[0] for op in Ops])
2216	OpT = np.array([op[1] for op in Ops])
2217	if SGData['SGInv'] and not SGData['SGFixed']:
2218	OpM = np.vstack((OpM,-OpM))
2219	OpT = np.vstack((OpT,-OpT))%1.
2220	OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2221	OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2222	Spn = SGData['SpnFlp'][:len(OpM)]
2223	Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2224	DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2225	PHKL = np.reshape(np.cos(twopinp.inner(HKL,OpT))Mag,(-1,))[:,nxs,nxs]*OpM #compute T(R,theta) eq(7)
2226	Ftest = np.random.rand(3) #random magnetic moment
2227	Psum = np.zeros(3)
2228	nsum = 0.
2229	nA = 0
2230	for dhkl,phkl in zip(DHKL,PHKL):
2231	if not np.allclose(dhkl,HKL): #test for eq(5)
2232	continue
2233	else:
2234	nA += 1
2235	nsum += np.trace(phkl) #eq(8)
2236	pterm = np.inner(Ftest,phkl) #eq(9)
2237	Psum += pterm
2238	if nsum/nA > 1.: #only need to look at nA=1 frok eq(8)
2239	return False
2240	if np.allclose(Psum,np.zeros(3)):
2241	return True
2242	else:
2243	if np.inner(HKL,Psum):
2244	return True
2245	return False
2246
2247	def checkSSextc(HKL,SSGData):
2248	Ops = SSGData['SSGOps']
2249	OpM = np.array([op[0] for op in Ops])
2250	OpT = np.array([op[1] for op in Ops])
2251	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2252	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2253	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2254	for dhkl,phkl in list(zip(DHKL,PHKL))[1:]: #skip identity
2255	if dhkl.any():
2256	continue
2257	else:
2258	if phkl%1.:
2259	return False
2260	return True
2261
2262	################################################################################
2263	#### Site symmetry tables
2264	################################################################################
2265
2266	OprName = {
2267	'-6643': ['-1',1],'6479' : ['2(z)',2],'-6479': ['m(z)',3],
2268	'6481' : ['m(y)',4],'-6481': ['2(y)',5],'6641' : ['m(x)',6],
2269	'-6641': ['2(x)',7],'6591' : ['m(+-0)',8],'-6591': ['2(+-0)',9],
2270	'6531' : ['m(110)',10],'-6531': ['2(110)',11],'6537' : ['4(z)',12],
2271	'-6537': ['-4(z)',13],'975' : ['3(111)',14],'6456' : ['3',15],
2272	'-489' : ['3(+--)',16],'483' : ['3(-+-)',17],'-969' : ['3(--+)',18],
2273	'819' : ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' : ['m(0+-)',21],
2274	'-2431': ['2(0+-)',22],'-657' : ['m(xz)',23],'657' : ['2(xz)',24],
2275	'1943' : ['-4(x)',25],'-1943': ['4(x)',26],'-2429': ['m(yz)',27],
2276	'2429' : ['2(yz)',28],'639' : ['-4(y)',29],'-639' : ['4(y)',30],
2277	'-6484': ['2(010)',4],'6484' : ['m(010)',5],'-6668': ['2(100)',6],
2278	'6668' : ['m(100)',7],'-6454': ['2(120)',18],'6454' : ['m(120)',19],
2279	'-6638': ['2(210)',20],'6638' : ['m(210)',21], #search in SytSym ends at m(210)
2280	'2223' : ['3(+++)2',39],
2281	'6538' : ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2282	'2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489' :['-3(+--)1',45],
2283	'801' : ['4(y)1',46],'1945' : ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2284	'6585' : ['4(z)3',49],'6584' : ['3(z)2',50],'6666' : ['6(z)5 ',51],
2285	'6643' : ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2286	'-6666': ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2287	'-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2288	'969' :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2289	'-2151':['-3(+--)2',64], }
2290
2291	KNsym = {
2292	'0' :' 1 ','1' :' -1 ','64' :' 2(x)','32' :' m(x)',
2293	'97' :' 2/m(x)','16' :' 2(y)','8' :' m(y)','25' :' 2/m(y)',
2294	'2' :' 2(z)','4' :' m(z)','7' :' 2/m(z)','134217728' :' 2(yz)',
2295	'67108864' :' m(yz)','201326593' :' 2/m(yz)','2097152' :' 2(0+-)','1048576' :' m(0+-)',
2296	'3145729' :'2/m(0+-)','8388608' :' 2(xz)','4194304' :' m(xz)','12582913' :' 2/m(xz)',
2297	'524288' :' 2(+0-)','262144' :' m(+0-)','796433' :'2/m(+0-)','1024' :' 2(xy)',
2298	'512' :' m(xy)','1537' :' 2/m(xy)','256' :' 2(+-0)','128' :' m(+-0)',
2299	'385' :'2/m(+-0)','76' :' mm2(x)','52' :' mm2(y)','42' :' mm2(z)',
2300	'135266336' :' mm2(yz)','69206048' :'mm2(0+-)','8650760' :' mm2(xz)','4718600' :'mm2(+0-)',
2301	'1156' :' mm2(xy)','772' :'mm2(+-0)','82' :' 222 ','136314944' :' 222(x)',
2302	'8912912' :' 222(y)','1282' :' 222(z)','127' :' mmm ','204472417' :' mmm(x)',
2303	'13369369' :' mmm(y)','1927' :' mmm(z)','33554496' :' 4(x)','16777280' :' -4(x)',
2304	'50331745' :'4/m(x)' ,'169869394' :'422(x)','84934738' :'-42m(x)','101711948' :'4mm(x)',
2305	'254804095' :'4/mmm(x)','536870928 ':' 4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2306	'545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2307	'2050' :' 4(z)','4098' :' -4(z)','6151' :'4/m(z)','3410' :'422(z)',
2308	'4818' :'-42m(z)','2730' :'4mm(z)','8191' :'4/mmm(z)','8192' :' 3(111)',
2309	'8193' :' -3(111)','2629888' :' 32(111)','1319040' :' 3m(111)','3940737' :'-3m(111)',
2310	'32768' :' 3(+--)','32769' :' -3(+--)','10519552' :' 32(+--)','5276160' :' 3m(+--)',
2311	'15762945' :'-3m(+--)','65536' :' 3(-+-)','65537' :' -3(-+-)','134808576' :' 32(-+-)',
2312	'67437056' :' 3m(-+-)','202180097' :'-3m(-+-)','131072' :' 3(--+)','131073' :' -3(--+)',
2313	'142737664' :' 32(--+)','71434368' :' 3m(--+)','214040961' :'-3m(--+)','237650' :' 23 ',
2314	'237695' :' m3 ','715894098' :' 432 ','358068946' :' -43m ','1073725439':' m3m ',
2315	'68157504' :' mm2(d100)','4456464' :' mm2(d010)','642' :' mm2(d001)','153092172' :'-4m2(x)',
2316	'277348404' :'-4m2(y)','5418' :'-4m2(z)','1075726335':' 6/mmm ','1074414420':'-6m2(100)',
2317	'1075070124':'-6m2(120)','1075069650':' 6mm ','1074414890':' 622 ','1073758215':' 6/m ',
2318	'1073758212':' -6 ','1073758210':' 6 ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2319	'1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2320	'1073758209':' -3 ','1073758208':' 3 ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2321	'1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2322	'1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2323	'1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2324	'1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2325	'1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2326	'1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2327	'1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2328	'1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2329	'1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':' -1 ',
2330	}
2331
2332	NXUPQsym = {
2333	'1' :(28,29,28,28),'-1' :( 1,29,28, 0),'2(x)' :(12,18,12,25),'m(x)' :(25,18,12,25),
2334	'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),
2335	'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),
2336	'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),
2337	'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),
2338	'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),
2339	'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),
2340	'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),
2341	'mm2(yz)' :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)' :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2342	'mm2(xy)' :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222' :( 1,10, 0,-1),'222(x)' :( 1,13, 0,-1),
2343	'222(y)' :( 1,12, 0,-1),'222(z)' :( 1,11, 0,-1),'mmm' :( 1,10, 0,-1),'mmm(x)' :( 1,13, 0,-1),
2344	'mmm(y)' :( 1,12, 0,-1),'mmm(z)' :( 1,11, 0,-1),'4(x)' :(12, 4,12, 0),'-4(x)' :( 1, 4,12, 0),
2345	'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),
2346	'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),
2347	'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,),
2348	'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),
2349	'-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),
2350	'-3(111)' :( 1, 5, 2, 0),'32(111)' :( 1, 5, 0, 2),'3m(111)' :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2351	'3(+--)' :( 5, 8, 5, 0),'-3(+--)' :( 1, 8, 5, 0),'32(+--)' :( 1, 8, 0, 5),'3m(+--)' :( 5, 8, 0, 5),
2352	'-3m(+--)' :( 1, 8, 0,-1),'3(-+-)' :( 4, 7, 4, 0),'-3(-+-)' :( 1, 7, 4, 0),'32(-+-)' :( 1, 7, 0, 4),
2353	'3m(-+-)' :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)' :( 3, 6, 3, 0),'-3(--+)' :( 1, 6, 3, 0),
2354	'32(--+)' :( 1, 6, 0, 3),'3m(--+)' :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23' :( 1, 1, 0, 0),
2355	'm3' :( 1, 1, 0, 0),'432' :( 1, 1, 0, 0),'-43m' :( 1, 1, 0, 0),'m3m' :( 1, 1, 0, 0),
2356	'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)' :( 1, 4, 0,-1),
2357	'-4m2(y)' :( 1, 3, 0,-1),'-4m2(z)' :( 1, 2, 0,-1),'6/mmm' :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2358	'-6m2(120)':( 1, 9, 0,-1),'6mm' :(14, 9, 0,-1),'622' :( 1, 9, 0,-1),'6/m' :( 1, 9,14,-1),
2359	'-6' :( 1, 9,14, 0),'6' :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2360	'3m(100)' :(14, 9, 0,14),'3m(120)' :(14, 9, 0,14),'32(100)' :( 1, 9, 0,14),'32(120)' :( 1, 9, 0,14),
2361	'-3' :( 1, 9,14, 0),'3' :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2362	'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2363	'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2364	'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2365	'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),
2366	'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),
2367	'm(001)' :(23,16,14,23),'m(100)' :(26,25,12,26),'m(010)' :(27,28,13,27),'m(110)' :(18,19, 6,18),
2368	'm(120)' :(24,27,15,24),'m(210)' :(25,26,16,25),'m(+-0)' :(17,20, 7,17),'2(001)' :(14,16,14,23),
2369	'2(100)' :(12,25,12,26),'2(010)' :(13,28,13,27),'2(110)' :( 6,19, 6,18),'2(120)' :(15,27,15,24),
2370	'2(210)' :(16,26,16,25),'2(+-0)' :( 7,20, 7,17),'-1' :( 1,29,28, 0)
2371	}
2372
2373	CSxinel = [[], # 0th empty - indices are Fortran style
2374	[[0,0,0],[ 0.0, 0.0, 0.0]], #1 0 0 0
2375	[[1,1,1],[ 1.0, 1.0, 1.0]], #2 X X X
2376	[[1,1,1],[ 1.0, 1.0,-1.0]], #3 X X -X
2377	[[1,1,1],[ 1.0,-1.0, 1.0]], #4 X -X X
2378	[[1,1,1],[ 1.0,-1.0,-1.0]], #5 -X X X
2379	[[1,1,0],[ 1.0, 1.0, 0.0]], #6 X X 0
2380	[[1,1,0],[ 1.0,-1.0, 0.0]], #7 X -X 0
2381	[[1,0,1],[ 1.0, 0.0, 1.0]], #8 X 0 X
2382	[[1,0,1],[ 1.0, 0.0,-1.0]], #9 X 0 -X
2383	[[0,1,1],[ 0.0, 1.0, 1.0]], #10 0 Y Y
2384	[[0,1,1],[ 0.0, 1.0,-1.0]], #11 0 Y -Y
2385	[[1,0,0],[ 1.0, 0.0, 0.0]], #12 X 0 0
2386	[[0,1,0],[ 0.0, 1.0, 0.0]], #13 0 Y 0
2387	[[0,0,1],[ 0.0, 0.0, 1.0]], #14 0 0 Z
2388	[[1,1,0],[ 1.0, 2.0, 0.0]], #15 X 2X 0
2389	[[1,1,0],[ 2.0, 1.0, 0.0]], #16 2X X 0
2390	[[1,1,2],[ 1.0, 1.0, 1.0]], #17 X X Z
2391	[[1,1,2],[ 1.0,-1.0, 1.0]], #18 X -X Z
2392	[[1,2,1],[ 1.0, 1.0, 1.0]], #19 X Y X
2393	[[1,2,1],[ 1.0, 1.0,-1.0]], #20 X Y -X
2394	[[1,2,2],[ 1.0, 1.0, 1.0]], #21 X Y Y
2395	[[1,2,2],[ 1.0, 1.0,-1.0]], #22 X Y -Y
2396	[[1,2,0],[ 1.0, 1.0, 0.0]], #23 X Y 0
2397	[[1,0,2],[ 1.0, 0.0, 1.0]], #24 X 0 Z
2398	[[0,1,2],[ 0.0, 1.0, 1.0]], #25 0 Y Z
2399	[[1,1,2],[ 1.0, 2.0, 1.0]], #26 X 2X Z
2400	[[1,1,2],[ 2.0, 1.0, 1.0]], #27 2X X Z
2401	[[1,2,3],[ 1.0, 1.0, 1.0]], #28 X Y Z
2402	]
2403
2404	CSuinel = [[], # 0th empty - indices are Fortran style
2405	[[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
2406	[[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
2407	[[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
2408	[[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
2409	[[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
2410	[[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
2411	[[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
2412	[[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
2413	[[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
2414	[[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
2415	[[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
2416	[[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
2417	[[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
2418	[[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
2419	[[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
2420	[[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
2421	[[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
2422	[[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
2423	[[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
2424	[[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
2425	[[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
2426	[[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
2427	[[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
2428	[[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
2429	[[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
2430	[[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
2431	[[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
2432	[[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
2433	[[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
2434	]
2435
2436	################################################################################
2437	#### Site symmetry routines
2438	################################################################################
2439
2440	def GetOprPtrName(key):
2441	'Needs a doc string'
2442	try:
2443	oprName = OprName[key][0]
2444	except KeyError:
2445	return key
2446	return oprName.replace('(','').replace(')','')
2447
2448	def GetOprPtrNumber(key):
2449	'Needs a doc string'
2450	try:
2451	return OprName[key][1]
2452	except KeyError:
2453	return key
2454
2455	def GetOprName(key):
2456	'Needs a doc string'
2457	return OprName[key][0]
2458
2459	def GetKNsym(key):
2460	'Needs a doc string'
2461	try:
2462	return KNsym[key].strip()
2463	except KeyError:
2464	return 'sp'
2465
2466	def GetNXUPQsym(siteSym):
2467	'''
2468	The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2469	'''
2470	return NXUPQsym[siteSym]
2471
2472	def GetCSxinel(siteSym):
2473	"returns Xyz terms, multipliers, GUI flags"
2474	indx = GetNXUPQsym(siteSym.strip())
2475	return CSxinel[indx[0]]
2476
2477	def GetCSuinel(siteSym):
2478	"returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2479	indx = GetNXUPQsym(siteSym.strip())
2480	return CSuinel[indx[1]]
2481
2482	def GetCSpqinel(SpnFlp,dupDir):
2483	"returns Mxyz terms, multipliers, GUI flags"
2484	CSI = [[1,2,3],[1.0,1.0,1.0]]
2485	for sopr in dupDir:
2486	# print (sopr,dupDir[sopr])
2487	opr = sopr.replace("'",'')
2488	indx = GetNXUPQsym(opr)
2489	if SpnFlp[dupDir[sopr]] > 0:
2490	csi = CSxinel[indx[2]] #P
2491	else:
2492	csi = CSxinel[indx[3]] #Q
2493	# print(opr,indx,csi,CSI)
2494	if not len(csi):
2495	return [[0,0,0],[0.,0.,0.]]
2496	for kcs in [0,1,2]:
2497	if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2498	jcs = CSI[0][kcs]
2499	for ics in [0,1,2]:
2500	if CSI[0][ics] == jcs:
2501	CSI[0][ics] = 0
2502	CSI[1][ics] = 0.
2503	elif CSI[0][ics] > jcs:
2504	CSI[0][ics] = CSI[0][ics]-1
2505	elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2506	CSI[1][kcs] = csi[1][kcs]
2507	elif CSI[0][kcs] >= csi[0][kcs]:
2508	CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2509	if CSI[1][kcs] != csi[1][kcs]:
2510	if CSI[1][kcs] == 1.:
2511	CSI[1][kcs] = csi[1][kcs]
2512	# print(CSI)
2513	return CSI
2514
2515	def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2516	phase = np.sum(XYZ*wave)
2517	ssopinv = nl.inv(ssop[0])
2518	mst = ssopinv[3][:3]
2519	epsinv = ssopinv[3][3]
2520	sdet = nl.det(sop[0])
2521	ssdet = nl.det(ssop[0])
2522	dtau = mst(XYZ-sop[1])-epsinvssop[1][3]
2523	dT = 1.0
2524	if np.any(dtau%.5):
2525	sumdtau = np.sum(dtau%.5)
2526	dT = 0.
2527	if np.abs(sumdtau-.5) > 1.e-4:
2528	dT = np.tan(np.pi*sumdtau)
2529	tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2530	return sdet,ssdet,dtau,dT,tauT
2531
2532	def OpsfromStringOps(A,SGData,SSGData):
2533	SGOps = SGData['SGOps']
2534	SSGOps = SSGData['SSGOps']
2535	Ax = A.split('+')
2536	Ax[0] = int(Ax[0])
2537	iC = 1
2538	if Ax[0] < 0:
2539	iC = -1
2540	iAx = abs(Ax[0])
2541	nA = iAx%100-1
2542	nC = iAx//100
2543	unit = [0,0,0]
2544	if len(Ax) > 1:
2545	unit = eval('['+Ax[1]+']')
2546	return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2547
2548	def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2549
2550	def orderParms(CSI):
2551	parms = [0,]
2552	for csi in CSI:
2553	for i in [0,1,2]:
2554	if csi[i] not in parms:
2555	parms.append(csi[i])
2556	for csi in CSI:
2557	for i in [0,1,2]:
2558	csi[i] = parms.index(csi[i])
2559	return CSI
2560
2561	def fracCrenel(tau,Toff,Twid):
2562	Tau = (tau-Toff[:,nxs])%1.
2563	A = np.where(Tau<Twid[:,nxs],1.,0.)
2564	return A
2565
2566	def fracFourier(tau,nH,fsin,fcos):
2567	SA = np.sin(2.nHnp.pi*tau)
2568	CB = np.cos(2.nHnp.pi*tau)
2569	A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2570	B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2571	return A+B
2572
2573	def posFourier(tau,nH,psin,pcos):
2574	SA = np.sin(2nHnp.pi*tau)
2575	CB = np.cos(2nHnp.pi*tau)
2576	A = SA[nxs,nxs,:]*psin[:,:,nxs]
2577	B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2578	return A+B
2579
2580	def posZigZag(tau,Tmm,XYZmax):
2581	DT = Tmm[1]-Tmm[0]
2582	slopeUp = 2.*XYZmax/DT
2583	slopeDn = 2.*XYZmax/(1.-DT)
2584	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])
2585	return A
2586
2587	def posBlock(tau,Tmm,XYZmax):
2588	A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2589	return A
2590
2591	def DoFrac():
2592	delt2 = np.eye(2)*0.001
2593	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2594	dFTP = []
2595	if siteSym == '1':
2596	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2597	elif siteSym == '-1':
2598	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2599	else:
2600	FSC = np.ones(2,dtype='i')
2601	CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2602	if 'Crenel' in waveType:
2603	dF = np.zeros_like(tau)
2604	else:
2605	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2606	dFT = np.zeros_like(dF)
2607	dFTP = []
2608	for i in SdIndx:
2609	sop = Sop[i]
2610	ssop = SSop[i]
2611	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2612	fsc = np.ones(2,dtype='i')
2613	if 'Crenel' in waveType:
2614	dFT = np.zeros_like(tau)
2615	fsc = [1,1]
2616	else: #Fourier
2617	dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2618	dFT = nl.det(sop[0])*dFT
2619	dFT = dFT[:,np.argsort(tauT)]
2620	dFT[0] *= ssdet
2621	dFT[1] *= sdet
2622	dFTP.append(dFT)
2623
2624	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2625	fsc = [1,1]
2626	if dT:
2627	CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2628	else:
2629	CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2630	FSC = np.zeros(2,dtype='i')
2631	return CSI,dF,dFTP
2632	else:
2633	for i in range(2):
2634	if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2635	fsc[i] = 1
2636	else:
2637	fsc[i] = 0
2638	FSC &= fsc
2639	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2640	n = -1
2641	for i,F in enumerate(FSC):
2642	if F:
2643	n += 1
2644	CSI[0][i] = n+1
2645	CSI[1][i] = 1.0
2646
2647	return CSI,dF,dFTP
2648
2649	def DoXYZ():
2650	delt5 = np.ones(5)*0.001
2651	delt6 = np.eye(6)*0.001
2652	if 'Fourier' in waveType:
2653	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2654	#3x6x12 modulated position array (X,Spos,tau)& force positive
2655	elif waveType in ['ZigZag','Block']:
2656	if waveType == 'ZigZag':
2657	dX = posZigZag(tau,delt5[:2],delt5[2:])
2658	else:
2659	dX = posBlock(tau,delt5[:2],delt5[2:])
2660	dXTP = []
2661	if siteSym == '1':
2662	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2663	elif siteSym == '-1':
2664	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.],]
2665	else:
2666	if 'Fourier' in waveType:
2667	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2668	elif waveType in ['ZigZag','Block']:
2669	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2670	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2671	XSC = np.ones(6,dtype='i')
2672	dXTP = []
2673	for i in SdIndx:
2674	sop = Sop[i]
2675	ssop = SSop[i]
2676	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2677	xsc = np.ones(6,dtype='i')
2678	if 'Fourier' in waveType:
2679	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2680	elif waveType == 'ZigZag':
2681	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2682	elif waveType == 'Block':
2683	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2684	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2685	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2686	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2687	dXTP.append(dXT)
2688	if waveType == 'Fourier':
2689	for i in range(3):
2690	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2691	xsc[i] = 0
2692	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2693	xsc[i+3] = 0
2694	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2695	xsc[3:6] = 0
2696	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2697	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2698	if dT:
2699	if '(x)' in siteSym:
2700	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2701	if 'm' in siteSym and len(SdIndx) == 1:
2702	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2703	elif '(y)' in siteSym:
2704	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2705	if 'm' in siteSym and len(SdIndx) == 1:
2706	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2707	elif '(z)' in siteSym:
2708	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2709	if 'm' in siteSym and len(SdIndx) == 1:
2710	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2711	else:
2712	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2713	if '4/mmm' in laue:
2714	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2715	if '(xy)' in siteSym:
2716	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2717	if dT:
2718	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2719	else:
2720	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2721	if '(xy)' in siteSym or '(+-0)' in siteSym:
2722	mul = 1
2723	if '(+-0)' in siteSym:
2724	mul = -1
2725	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2726	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2727	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2728	xsc[3:5] = 0
2729	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2730	CSI[0][:2] = [[12,0,0],[12,0,0]]
2731	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2732	xsc[:2] = 0
2733	else:
2734	for i in range(3):
2735	if not np.allclose(dX[:,i],dXT[i,:,i]):
2736	xsc[i] = 0
2737	XSC &= xsc
2738	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2739	if waveType == 'Fourier':
2740	n = -1
2741	if debug: print (XSC)
2742	for i,X in enumerate(XSC):
2743	if X:
2744	n += 1
2745	CSI[0][i][0] = n+1
2746	CSI[1][i][0] = 1.0
2747
2748	return list(CSI),dX,dXTP
2749
2750	def DoUij():
2751	delt12 = np.eye(12)*0.0001
2752	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2753	dUTP = []
2754	if siteSym == '1':
2755	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2756	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2757	elif siteSym == '-1':
2758	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2759	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2760	else:
2761	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2762	USC = np.ones(12,dtype='i')
2763	dUTP = []
2764	dtau = 0.
2765	for i in SdIndx:
2766	sop = Sop[i]
2767	ssop = SSop[i]
2768	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2769	usc = np.ones(12,dtype='i')
2770	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2771	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2772	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2773	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2774	dUT = dUT[:,:,np.argsort(tauT)]
2775	dUT[:,:6,:] =(ssdetsdet)
2776	dUTP.append(dUT)
2777	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2778	if dT:
2779	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2780	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2781	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2782	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2783	else:
2784	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2785	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2786	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2787	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2788	if 'mm2(x)' in siteSym and dT:
2789	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2790	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2791	elif '(xy)' in siteSym and dT:
2792	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2793	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2794	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2795	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2796	elif '(x)' in siteSym and dT:
2797	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2798	elif '(y)' in siteSym and dT:
2799	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2800	elif '(z)' in siteSym and dT:
2801	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2802	for i in range(6):
2803	if not USC[i]:
2804	CSI[0][i] = [0,0,0]
2805	CSI[1][i] = [0.,0.,0.]
2806	CSI[0][i+6] = [0,0,0]
2807	CSI[1][i+6] = [0.,0.,0.]
2808	else:
2809	for i in range(6):
2810	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2811	usc[i] = 0
2812	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2813	usc[i+6] = 0
2814	if np.any(dUT[1,0,:]):
2815	if '4/m' in siteSym:
2816	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2817	if ssop[1][3]:
2818	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2819	usc[9] = 1
2820	else:
2821	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2822	usc[9] = 0
2823	elif '4' in siteSym:
2824	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2825	CSI[0][:2] = [[11,0,0],[11,0,0]]
2826	if ssop[1][3]:
2827	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2828	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2829	usc[2] = 0
2830	usc[8] = 0
2831	usc[3] = 1
2832	usc[9] = 1
2833	else:
2834	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2835	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2836	usc[2] = 1
2837	usc[8] = 1
2838	usc[3] = 0
2839	usc[9] = 0
2840	elif 'xy' in siteSym or '+-0' in siteSym:
2841	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2842	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2843	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2844	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2845	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2846	usc[4:6] = 0
2847	usc[6:8] = 0
2848
2849	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2850	USC &= usc
2851	if debug: print (USC)
2852	if not np.any(dtau%.5):
2853	n = -1
2854	for i,U in enumerate(USC):
2855	if U:
2856	n += 1
2857	CSI[0][i][0] = n+1
2858	CSI[1][i][0] = 1.0
2859
2860	return list(CSI),dU,dUTP
2861
2862	def DoMag():
2863	delt6 = np.eye(6)*0.001
2864	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2865	dMTP = []
2866	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2867	if siteSym == '1':
2868	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2869	elif siteSym in ['-1','mmm',]:
2870	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2871	elif siteSym in ['4(z)','422(z)']:
2872	CSI[0][0][0] = CSI[0][4][1] = 1
2873	CSI[1][0][0] = 1.0
2874	CSI[1][4][1] = -1.0
2875	elif siteSym in ['-4m2(z)','422(z)',]:
2876	CSI[0][5][0] = 1
2877	CSI[1][5][0] = 1.0
2878	elif siteSym in ['-32(100)','-3',]:
2879	CSI[0][2][0] = 1
2880	CSI[1][2][0] = 1.0
2881	elif siteSym in ['3',]:
2882	CSI[0][0][0] = CSI[0][3][0] = CSI[0][4][0] = 1
2883	CSI[1][0][0] = -np.sqrt(3.0)
2884	CSI[1][3][0] = 2.0
2885	CSI[1][4][0] = 1.0
2886	elif siteSym in ['622','2(100)','32(100)',]:
2887	CSI[0][0][0] = CSI[0][1][0] = CSI[0][3][0] = 1
2888	CSI[1][0][0] = 1.0
2889	CSI[1][1][0] = 2.0
2890	CSI[1][3][0] = np.sqrt(3.0)
2891	else:
2892	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2893	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2894	MSC = np.ones(6,dtype='i')
2895	dMTP = []
2896	for i in SdIndx:
2897	sop = Sop[i]
2898	ssop = SSop[i]
2899	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2900	msc = np.ones(6,dtype='i')
2901	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2902	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2903	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2904	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2905	dMTP.append(dMT)
2906	for i in range(3):
2907	if not np.allclose(dM[i,i,:],sdet*dMT[i,i,:]):
2908	msc[i] = 0
2909	if not np.allclose(dM[i,i+3,:],sdet*dMT[i,i+3,:]):
2910	msc[i+3] = 0
2911	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2912	msc[3:6] = 0
2913	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2914	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2915	if dT:
2916	if '(x)' in siteSym:
2917	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2918	if 'm' in siteSym and len(SdIndx) == 1:
2919	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2920	elif '(y)' in siteSym:
2921	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2922	if 'm' in siteSym and len(SdIndx) == 1:
2923	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2924	elif '(z)' in siteSym:
2925	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2926	if 'm' in siteSym and len(SdIndx) == 1:
2927	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2928	else:
2929	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2930	if '4/mmm' in laue:
2931	if siteSym in ['4/mmm(z)',]:
2932	CSI = 3[[0,0,0],]+[[0,0,0],[0,0,0],[1,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2933	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2934	if '(xy)' in siteSym:
2935	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2936	if dT:
2937	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2938	else:
2939	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2940	if '(xy)' in siteSym or '(+-0)' in siteSym:
2941	mul = 1
2942	if '(+-0)' in siteSym:
2943	mul = -1
2944	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2945	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2946	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2947	msc[3:5] = 0
2948	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2949	CSI[0][:2] = [[12,0,0],[12,0,0]]
2950	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2951	msc[:2] = 0
2952	MSC &= msc
2953	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2954	n = -1
2955	if debug: print (MSC)
2956	for i,M in enumerate(MSC):
2957	if M:
2958	n += 1
2959	CSI[0][i][0] = n+1
2960	CSI[1][i][0] = 1.0
2961
2962	return list(CSI),dM,dMTP
2963
2964	if debug: print ('super space group: '+SSGData['SSpGrp'])
2965	xyz = np.array(XYZ)%1.
2966	SGOps = copy.deepcopy(SGData['SGOps'])
2967	laue = SGData['SGLaue']
2968	siteSym = SytSym(XYZ,SGData)[0].strip()
2969	if debug: print ('siteSym: '+siteSym)
2970	SSGOps = copy.deepcopy(SSGData['SSGOps'])
2971	#expand ops to include inversions if any
2972	if SGData['SGInv'] and not SGData['SGFixed']:
2973	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2974	SGOps.append([-op[0],-op[1]%1.])
2975	SSGOps.append([-sop[0],-sop[1]%1.])
2976	#build set of sym ops around special position
2977	SSop = []
2978	Sop = []
2979	Sdtau = []
2980	for iop,Op in enumerate(SGOps):
2981	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2982	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2983	SSop.append(SSGOps[iop])
2984	Sop.append(SGOps[iop])
2985	ssopinv = nl.inv(SSGOps[iop][0])
2986	mst = ssopinv[3][:3]
2987	epsinv = ssopinv[3][3]
2988	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
2989	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
2990	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2991	#setup displacement arrays
2992	tau = np.linspace(-1,1,49,True)
2993	#make modulation arrays - one parameter at a time
2994	if Stype == 'Sfrac':
2995	CSI,dF,dFTP = DoFrac()
2996	elif Stype == 'Spos':
2997	CSI,dF,dFTP = DoXYZ()
2998	CSI[0] = orderParms(CSI[0])
2999	elif Stype == 'Sadp':
3000	CSI,dF,dFTP = DoUij()
3001	CSI[0] = orderParms(CSI[0])
3002	elif Stype == 'Smag':
3003	CSI,dF,dFTP = DoMag()
3004
3005	if debug:
3006	return CSI,dF,dFTP
3007	else:
3008	return CSI,[],[]
3009
3010	def MustrainNames(SGData):
3011	'Needs a doc string'
3012	laue = SGData['SGLaue']
3013	uniq = SGData['SGUniq']
3014	if laue in ['m3','m3m']:
3015	return ['S400','S220']
3016	elif laue in ['6/m','6/mmm','3m1']:
3017	return ['S400','S004','S202']
3018	elif laue in ['31m','3']:
3019	return ['S400','S004','S202','S301']
3020	elif laue in ['3R','3mR']:
3021	return ['S400','S220','S310','S211']
3022	elif laue in ['4/m','4/mmm']:
3023	return ['S400','S004','S220','S022']
3024	elif laue in ['mmm']:
3025	return ['S400','S040','S004','S220','S202','S022']
3026	elif laue in ['2/m']:
3027	SHKL = ['S400','S040','S004','S220','S202','S022']
3028	if uniq == 'a':
3029	SHKL += ['S013','S031','S211']
3030	elif uniq == 'b':
3031	SHKL += ['S301','S103','S121']
3032	elif uniq == 'c':
3033	SHKL += ['S130','S310','S112']
3034	return SHKL
3035	else:
3036	SHKL = ['S400','S040','S004','S220','S202','S022']
3037	SHKL += ['S310','S103','S031','S130','S301','S013']
3038	SHKL += ['S211','S121','S112']
3039	return SHKL
3040
3041	def HStrainVals(HSvals,SGData):
3042	laue = SGData['SGLaue']
3043	uniq = SGData['SGUniq']
3044	DIJ = np.zeros(6)
3045	if laue in ['m3','m3m']:
3046	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
3047	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3048	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
3049	elif laue in ['3R','3mR']:
3050	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3051	elif laue in ['4/m','4/mmm']:
3052	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3053	elif laue in ['mmm']:
3054	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3055	elif laue in ['2/m']:
3056	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3057	if uniq == 'a':
3058	DIJ[5] = HSvals[3]
3059	elif uniq == 'b':
3060	DIJ[4] = HSvals[3]
3061	elif uniq == 'c':
3062	DIJ[3] = HSvals[3]
3063	else:
3064	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3065	return DIJ
3066
3067	def HStrainNames(SGData):
3068	'Needs a doc string'
3069	laue = SGData['SGLaue']
3070	uniq = SGData['SGUniq']
3071	if laue in ['m3','m3m']:
3072	return ['D11','eA'] #add cubic strain term
3073	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3074	return ['D11','D33']
3075	elif laue in ['3R','3mR']:
3076	return ['D11','D12']
3077	elif laue in ['4/m','4/mmm']:
3078	return ['D11','D33']
3079	elif laue in ['mmm']:
3080	return ['D11','D22','D33']
3081	elif laue in ['2/m']:
3082	Dij = ['D11','D22','D33']
3083	if uniq == 'a':
3084	Dij += ['D23']
3085	elif uniq == 'b':
3086	Dij += ['D13']
3087	elif uniq == 'c':
3088	Dij += ['D12']
3089	return Dij
3090	else:
3091	Dij = ['D11','D22','D33','D12','D13','D23']
3092	return Dij
3093
3094	def MustrainCoeff(HKL,SGData):
3095	'Needs a doc string'
3096	#NB: order of terms is the same as returned by MustrainNames
3097	laue = SGData['SGLaue']
3098	uniq = SGData['SGUniq']
3099	h,k,l = HKL
3100	Strm = []
3101	if laue in ['m3','m3m']:
3102	Strm.append(h4+k4+l**4)
3103	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3104	elif laue in ['6/m','6/mmm','3m1']:
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	elif laue in ['31m','3']:
3109	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3110	Strm.append(l**4)
3111	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3112	Strm.append(4.0lh**3)
3113	elif laue in ['3R','3mR']:
3114	Strm.append(h4+k4+l**4)
3115	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3116	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
3117	Strm.append(4.0(klh2+hlk2+hkl*2))
3118	elif laue in ['4/m','4/mmm']:
3119	Strm.append(h4+k4)
3120	Strm.append(l**4)
3121	Strm.append(3.0(hk)**2)
3122	Strm.append(3.0((hl)*2+(kl)**2))
3123	elif laue in ['mmm']:
3124	Strm.append(h**4)
3125	Strm.append(k**4)
3126	Strm.append(l**4)
3127	Strm.append(3.0(hk)**2)
3128	Strm.append(3.0(hl)**2)
3129	Strm.append(3.0(kl)**2)
3130	elif laue in ['2/m']:
3131	Strm.append(h**4)
3132	Strm.append(k**4)
3133	Strm.append(l**4)
3134	Strm.append(3.0(hk)**2)
3135	Strm.append(3.0(hl)**2)
3136	Strm.append(3.0(kl)**2)
3137	if uniq == 'a':
3138	Strm.append(2.0kl**3)
3139	Strm.append(2.0lk**3)
3140	Strm.append(4.0klh*2)
3141	elif uniq == 'b':
3142	Strm.append(2.0lh**3)
3143	Strm.append(2.0hl**3)
3144	Strm.append(4.0hlk*2)
3145	elif uniq == 'c':
3146	Strm.append(2.0hk**3)
3147	Strm.append(2.0kh**3)
3148	Strm.append(4.0hkl*2)
3149	else:
3150	Strm.append(h**4)
3151	Strm.append(k**4)
3152	Strm.append(l**4)
3153	Strm.append(3.0(hk)**2)
3154	Strm.append(3.0(hl)**2)
3155	Strm.append(3.0(kl)**2)
3156	Strm.append(2.0kh**3)
3157	Strm.append(2.0hl**3)
3158	Strm.append(2.0lk**3)
3159	Strm.append(2.0hk**3)
3160	Strm.append(2.0lh**3)
3161	Strm.append(2.0kl**3)
3162	Strm.append(4.0klh*2)
3163	Strm.append(4.0hlk*2)
3164	Strm.append(4.0khl*2)
3165	return Strm
3166
3167	def MuShklMean(SGData,Amat,Shkl):
3168
3169	def genMustrain(xyz,Shkl):
3170	uvw = np.inner(Amat.T,xyz)
3171	Strm = np.array(MustrainCoeff(uvw,SGData))
3172	Sum = np.sum(np.multiply(Shkl,Strm))
3173	Sum = np.where(Sum > 0.01,Sum,0.01)
3174	Sum = np.sqrt(Sum)
3175	return Sum*xyz
3176
3177	PHI = np.linspace(0.,360.,30,True)
3178	PSI = np.linspace(0.,180.,30,True)
3179	X = np.outer(npcosd(PHI),npsind(PSI))
3180	Y = np.outer(npsind(PHI),npsind(PSI))
3181	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3182	XYZ = np.dstack((X,Y,Z))
3183	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3184	return np.sqrt(np.sum(XYZ**2)/900.)
3185
3186	def Muiso2Shkl(muiso,SGData,cell):
3187	"this is to convert isotropic mustrain to generalized Shkls"
3188	import GSASIIlattice as G2lat
3189	A = G2lat.cell2AB(cell)[0]
3190
3191	def minMus(Shkl,muiso,H,SGData,A):
3192	U = np.inner(A.T,H)
3193	S = np.array(MustrainCoeff(U,SGData))
3194	nS = S.shape[0]
3195	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:nS,nxs]),axis=0))
3196	rad = np.sqrt(np.sum((Sum[:,nxs]H)*2,axis=1))
3197	return (muiso-rad)**2
3198
3199	laue = SGData['SGLaue']
3200	PHI = np.linspace(0.,360.,60,True)
3201	PSI = np.linspace(0.,180.,60,True)
3202	X = np.outer(npsind(PHI),npsind(PSI))
3203	Y = np.outer(npcosd(PHI),npsind(PSI))
3204	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3205	HKL = np.dstack((X,Y,Z))
3206	if laue in ['m3','m3m']:
3207	S0 = [1000.,1000.]
3208	elif laue in ['6/m','6/mmm']:
3209	S0 = [1000.,1000.,1000.]
3210	elif laue in ['31m','3','3m1']:
3211	S0 = [1000.,1000.,1000.,1000.]
3212	elif laue in ['3R','3mR']:
3213	S0 = [1000.,1000.,1000.,1000.]
3214	elif laue in ['4/m','4/mmm']:
3215	S0 = [1000.,1000.,1000.,1000.]
3216	elif laue in ['mmm']:
3217	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3218	elif laue in ['2/m']:
3219	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3220	else:
3221	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
3222	1000.,1000.,0.,0.,0.]
3223	S0 = np.array(S0)
3224	HKL = np.reshape(HKL,(-1,3))
3225	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3226	return result[0]
3227
3228	def PackRot(SGOps):
3229	IRT = []
3230	for ops in SGOps:
3231	M = ops[0]
3232	irt = 0
3233	for j in range(2,-1,-1):
3234	for k in range(2,-1,-1):
3235	irt *= 3
3236	irt += M[k][j]
3237	IRT.append(int(irt))
3238	return IRT
3239
3240	def SytSym(XYZ,SGData):
3241	'''
3242	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3243
3244	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
3245	:param SGData: from SpcGroup
3246	:Returns: a four element tuple:
3247
3248	* The 1st element is a code for the site symmetry (see GetKNsym)
3249	* The 2nd element is the site multiplicity
3250	* Ndup number of overlapping operators
3251	* dupDir Dict - dictionary of overlapping operators
3252
3253	'''
3254	Mult = 1
3255	Isym = 0
3256	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3257	Isym = 1073741824
3258	Jdup = 0
3259	Ndup = 0
3260	dupDir = {}
3261	inv = SGData['SGInv']+1
3262	icen = SGData['SGCen']
3263	Ncen = len(icen)
3264	if SGData['SGFixed']: #already in list of operators
3265	inv = 1
3266	if SGData['SGGray'] and Ncen > 1: Ncen //= 2
3267	Xeqv = list(GenAtom(XYZ,SGData,True))
3268	# for xeqv in Xeqv: print(xeqv)
3269	IRT = PackRot(SGData['SGOps'])
3270	L = -1
3271	for ic,cen in enumerate(icen[:Ncen]):
3272	for invers in range(int(inv)):
3273	for io,ops in enumerate(SGData['SGOps']):
3274	irtx = (1-2invers)IRT[io]
3275	L += 1
3276	if not Xeqv[L][1]:
3277	Ndup = io
3278	Jdup += 1
3279	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
3280	if jx < 39:
3281	px = GetOprName(str(irtx))
3282	if Xeqv[L][-1] < 0:
3283	if '(' in px:
3284	px = px.split('(')
3285	px[0] += "'"
3286	px = '('.join(px)
3287	else:
3288	px += "'"
3289	dupDir[px] = L
3290	Isym += 2**(jx-1)
3291	if Isym == 1073741824: Isym = 0
3292	Mult = len(SGData['SGOps'])Nceninv//Jdup
3293
3294	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3295
3296	def MagSytSym(SytSym,dupDir,SGData):
3297	'''
3298	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3299	'''
3300	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3301	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3302	# print('dupDir',dupDir)
3303	SplitSytSym = SytSym.split('(')
3304	if SGData['SGGray']:
3305	return SytSym+"1'"
3306	if SytSym == '1': #genersl position
3307	return SytSym
3308	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
3309	try:
3310	MagSytSym = SGData['MagSpGrp'].split()[1]
3311	except IndexError:
3312	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3313	if len(SplitSytSym) > 1:
3314	MagSytSym += '('+SplitSytSym[1]
3315	return MagSytSym
3316	if len(dupDir) == 1:
3317	return list(dupDir.keys())[0]
3318
3319
3320	if '2/m' in SytSym: #done I think; last 2wo might be not needed
3321	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3322	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3323	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3324	'(110)':['2(110)','m(110)']}
3325
3326	elif '4/mmm' in SytSym:
3327	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
3328	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
3329	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
3330	elif '4mm' in SytSym:
3331	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3332	elif '422' in SytSym:
3333	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3334	elif '-4m2' in SytSym:
3335	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3336	elif '-42m' in SytSym:
3337	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3338	elif '-4' in SytSym:
3339	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3340	elif '4' in SytSym:
3341	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3342
3343	elif '222' in SytSym:
3344	ops = {'':['2(x)','2(y)','2(z)'],
3345	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3346	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3347	'(110)':['2(z)','2(110)','2(+-0)',],}
3348	elif 'mm2' in SytSym:
3349	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3350	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3351	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3352	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3353	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3354	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3355	'(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3356	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3357	'(110)':['m(z)','m(+-0)','2(110)',],}
3358	elif 'mmm' in SytSym:
3359	ops = {'':['m(x)','m(y)','m(z)'],
3360	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3361	'(110)':['m(z)','m(110)','m(+-0)',],
3362	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3363
3364	elif '32' in SytSym:
3365	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)'],'(111)':['3(111)','2(x)']}
3366	elif '23' in SytSym:
3367	ops = {'':['2(x)','3(111)']}
3368	elif 'm3' in SytSym:
3369	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3370	elif '3m' in SytSym:
3371	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3372	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3373	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3374
3375	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m','-43m',]: #not simple cases
3376	MagSytSym = SytSym
3377	if "-1'" in dupDir:
3378	if '-6' in SytSym:
3379	MagSytSym = MagSytSym.replace('-6',"-6'")
3380	elif '-3m' in SytSym:
3381	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3382	elif '-3' in SytSym:
3383	MagSytSym = MagSytSym.replace('-3',"-3'")
3384	elif '-6m2' in SytSym:
3385	if "m'(110)" in dupDir:
3386	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3387	elif '6/m' in SytSym:
3388	if "m'(z)" in dupDir:
3389	MagSytSym = "6'/m'"
3390	elif '6mm' in SytSym:
3391	if "m'(110)" in dupDir:
3392	MagSytSym = "6'm'm"
3393	elif '-43m' in SytSym:
3394	if "m'(110)" in dupDir:
3395	MagSytSym = "-43m'"
3396	return MagSytSym
3397	try:
3398	axis = '('+SytSym.split('(')[1]
3399	except IndexError:
3400	axis = ''
3401	MagSytSym = ''
3402	for m in ops[axis]:
3403	if m in dupDir:
3404	MagSytSym += m.split('(')[0]
3405	else:
3406	MagSytSym += m.split('(')[0]+"'"
3407	if '2/m' in SytSym and '2' in m:
3408	MagSytSym += '/'
3409	if '-3/m' in SytSym:
3410	MagSytSym = '-'+MagSytSym
3411
3412	MagSytSym += axis
3413	# some exceptions & special rules
3414	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3415	return MagSytSym
3416
3417	# if len(GenSym) == 3:
3418	# if SGSpin[1] < 0:
3419	# if 'mm2' in SytSym:
3420	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3421	# else: #bad rule for I41/a
3422	# MagSytSym = SplitSytSym[0]+"'"
3423	# if len(SplitSytSym) > 1:
3424	# MagSytSym += '('+SplitSytSym[1]
3425	# else:
3426	# MagSytSym = SytSym
3427	# if len(SplitSytSym) >1:
3428	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3429	# MagSytSym = MagSytSym.replace('-4',"-4'")
3430	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3431	# MagSytSym = MagSytSym.replace('-6',"-6'")
3432	# return MagSytSym
3433	#
3434	return SytSym
3435
3436	def UpdateSytSym(Phase):
3437	''' Update site symmetry/site multiplicity after space group/BNS lattice change
3438	'''
3439	generalData = Phase['General']
3440	SGData = generalData['SGData']
3441	Atoms = Phase['Atoms']
3442	cx,ct,cs,cia = generalData['AtomPtrs']
3443	for atom in Atoms:
3444	XYZ = atom[cx:cx+3]
3445	sytsym,Mult = SytSym(XYZ,SGData)[:2]
3446	sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3447	atom[cs] = sytsym
3448	if generalData['Type'] == 'magnetic':
3449	magSytSym = MagSytSym(sytSym,dupDir,SGData)
3450	atom[cs] = magSytSym
3451	atom[cs+1] = Mult
3452	return
3453
3454	def ElemPosition(SGData):
3455	''' Under development.
3456	Object here is to return a list of symmetry element types and locations suitable
3457	for say drawing them.
3458	So far I have the element type... getting all possible locations without lookup may be impossible!
3459	'''
3460	Inv = SGData['SGInv']
3461	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3462	# get operators & expand if centrosymmetric
3463	SymElements = []
3464	Ops = SGData['SGOps']
3465	opM = np.array([op[0].T for op in Ops])
3466	opT = np.array([op[1] for op in Ops])
3467	if Inv:
3468	opM = np.concatenate((opM,-opM))
3469	opT = np.concatenate((opT,-opT))
3470	opMT = list(zip(opM,opT))
3471	for M,T in opMT[1:]: #skip I
3472	Dt = int(nl.det(M))
3473	Tr = int(np.trace(M))
3474	Dt = -(Dt-1)//2
3475	Es = eleSym[Tr][Dt]
3476	if Dt: #rotation-inversion
3477	I = np.eye(3)
3478	if Tr == 1: #mirrors/glides
3479	if np.any(T): #glide
3480	M2 = np.inner(M,M)
3481	MT = np.inner(M,T)+T
3482	print ('glide',Es,MT)
3483	print (M2)
3484	else: #mirror
3485	print ('mirror',Es,T)
3486	print (I-M)
3487	X = [-1,-1,-1]
3488	elif Tr == -3: # pure inversion
3489	X = np.inner(nl.inv(I-M),T)
3490	print ('inversion',Es,X)
3491	else: #other rotation-inversion
3492	M2 = np.inner(M,M)
3493	MT = np.inner(M,T)+T
3494	print ('rot-inv',Es,MT)
3495	print (M2)
3496	X = [-1,-1,-1]
3497	else: #rotations
3498	print ('rotation',Es)
3499	X = [-1,-1,-1]
3500	SymElements.append([Es,X])
3501
3502	return SymElements
3503
3504	def ApplyStringOps(A,SGData,X,Uij=[]):
3505	'Needs a doc string'
3506	SGOps = SGData['SGOps']
3507	SGCen = SGData['SGCen']
3508	Ax = A.split('+')
3509	Ax[0] = int(Ax[0])
3510	iC = 1
3511	if Ax[0] < 0:
3512	iC = -1
3513	Ax[0] = abs(Ax[0])
3514	nA = Ax[0]%100-1
3515	cA = Ax[0]//100
3516	Cen = SGCen[cA]
3517	M,T = SGOps[nA]
3518	if len(Ax)>1:
3519	cellA = Ax[1].split(',')
3520	cellA = np.array([int(a) for a in cellA])
3521	else:
3522	cellA = np.zeros(3)
3523	newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3524	if len(Uij):
3525	U = Uij2U(Uij)
3526	U = np.inner(M,np.inner(U,M).T)
3527	newUij = U2Uij(U)
3528	return [newX,newUij]
3529	else:
3530	return newX
3531
3532	def ApplyStringOpsMom(A,SGData,Mom):
3533	'Needs a doc string'
3534	SGOps = SGData['SGOps']
3535	Ax = A.split('+')
3536	Ax[0] = int(Ax[0])
3537	iAx = abs(Ax[0])
3538	nA = iAx%100-1
3539	if SGData['SGInv'] and not SGData['SGFixed']:
3540	nC = 2len(SGOps)(iAx//100)
3541	else:
3542	nC = len(SGOps)*(iAx//100)
3543	NA = nA
3544	if Ax[0] < 0:
3545	NA += len(SGOps)
3546	M,T = SGOps[nA]
3547	if SGData['SGGray']:
3548	newMom = -np.inner(Mom,M).Tnl.det(M)SGData['SpnFlp'][NA+nC]
3549	else:
3550	newMom = np.inner(Mom,M).Tnl.det(M)SGData['SpnFlp'][NA+nC]
3551	# print(len(SGOps),Ax[0],iAx,nC,nA,NA,MT2text([M,T]).replace(' ',''),SGData['SpnFlp'][NA],Mom,newMom)
3552	# print(Mom,newMom,MT2text([M,T]),)
3553	return newMom
3554
3555	def StringOpsProd(A,B,SGData):
3556	"""
3557	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3558	where '-' indicates inversion, c(>0) is the cell centering operator,
3559	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3560	Should return resultant string - C. SGData - dictionary using entries:
3561
3562	* 'SGCen': cell centering vectors [0,0,0] at least
3563	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3564
3565	"""
3566	SGOps = SGData['SGOps']
3567	SGCen = SGData['SGCen']
3568	#1st split out the cell translation part & work on the operator parts
3569	Ax = A.split('+'); Bx = B.split('+')
3570	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3571	iC = 0
3572	if Ax[0]*Bx[0] < 0:
3573	iC = 1
3574	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3575	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3576	cA = Ax[0]//100; cB = Bx[0]//100
3577	Cen = (SGCen[cA]+SGCen[cB])%1.0
3578	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3579	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3580	Mc = np.inner(Ma,Mb.T)
3581	# print Ma,Mb,Mc
3582	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3583	# print Ta,Tb,Tc
3584	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3585	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3586	#now the cell translation part
3587	if len(Ax)>1:
3588	cellA = Ax[1].split(',')
3589	cellA = [int(a) for a in cellA]
3590	else:
3591	cellA = [0,0,0]
3592	if len(Bx)>1:
3593	cellB = Bx[1].split(',')
3594	cellB = [int(b) for b in cellB]
3595	else:
3596	cellB = [0,0,0]
3597	cellC = np.add(cellA,cellB)
3598	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3599	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3600	return C
3601
3602	def U2Uij(U):
3603	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3604	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3605
3606	def Uij2U(Uij):
3607	#returns the thermal motion tensor U from Uij as numpy array
3608	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3609
3610	def StandardizeSpcName(spcgroup):
3611	'''Accept a spacegroup name where spaces may have not been used
3612	in the names according to the GSAS convention (spaces between symmetry
3613	for each axis) and return the space group name as used in GSAS
3614	'''
3615	rspc = spcgroup.replace(' ','').upper()
3616	# deal with rhombohedral and hexagonal setting designations
3617	rhomb = ''