Context Navigation

source: trunk/GSASIIspc.py @ 3703

Last change on this file since 3703 was 3703, checked in by vondreele, 7 years ago
refactoring of GetSSfxuinel
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 170.3 KB

Line
1	# -- coding: utf-8 --
2	"""
3	GSASIIspc: Space group module
4	-------------------------------
5
6	Space group interpretation routines. Note that space group information is
7	stored in a :ref:`Space Group (SGData)<SGData_table>` object.
8
9	"""
10	########### SVN repository information ###################
11	# $Date: 2018-10-24 21:30:34 +0000 (Wed, 24 Oct 2018) $
12	# $Author: vondreele $
13	# $Revision: 3703 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 3703 2018-10-24 21:30:34Z 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: 3703 $")
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): SGText[-1] += " 1'"
411	CentStr = 'centrosymmetric'
412	if not SGData['SGInv']:
413	CentStr = 'non'+CentStr
414	if SGData['SGLatt'] in 'ABCIFR':
415	SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
416	else:
417	SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
418	SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
419	if 'SGPtGrp' in SGData: #patch
420	SGText.append(' The lattice point group is '+SGData['SGPtGrp'])
421	SGText.append(' Multiplicity of a general site is '+str(Mult))
422	if SGData['SGUniq'] in ['a','b','c']:
423	SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
424	if SGData['SGInv']:
425	SGText.append(' The inversion center is located at 0,0,0')
426	if SGData['SGPolax']:
427	SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
428	SGText.append(' ')
429	if len(SGData['SGCen']) == 1:
430	SGText.append(' The equivalent positions are:\n')
431	else:
432	SGText.append(' The equivalent positions are:\n')
433	SGText.append(' ('+Latt2text(SGData['SGCen'])+')+\n')
434	SGTable = []
435	for i,Opr in enumerate(SGData['SGOps']):
436	SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
437	if AddInv and SGData['SGInv']:
438	for i,Opr in enumerate(SGData['SGOps']):
439	IOpr = [-Opr[0],-Opr[1]]
440	SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))
441	# if SGData.get('SGGray',False) and not SGData.get('SGFixed',False):
442	# SGTable.append(" for 1'")
443	# for i,Opr in enumerate(SGData['SGOps']):
444	# SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
445	# if AddInv and SGData['SGInv']:
446	# for i,Opr in enumerate(SGData['SGOps']):
447	# IOpr = [-Opr[0],-Opr[1]]
448	# SGTable.append('(%2d) %s'%(i+1,MT2text(IOpr)))
449	return SGText,SGTable
450
451	def AllOps(SGData):
452	'''
453	Returns a list of all operators for a space group, including those for
454	centering and a center of symmetry
455
456	:param SGData: from :func:`SpcGroup`
457	:returns: (SGTextList,offsetList,symOpList,G2oprList) where
458
459	* SGTextList: a list of strings with formatted and normalized
460	symmetry operators.
461	* offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
462	symmetry operation to the operator in SGTextList and symOpList.
463	these dx (etc.) values are added to the GSAS-II generated
464	positions to provide the positions that are generated
465	by the normalized symmetry operators.
466	* symOpList: a list of tuples with the normalized symmetry
467	operations as (M,T) values
468	(see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
469	* G2oprList: The GSAS-II operations for each symmetry operation as
470	a tuple with (center,mult,opnum,opcode), where center is (0,0,0), (0.5,0,0),
471	(0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
472	where opnum is the number for the symmetry operation, in ``SGOps``
473	(starting with 0) and opcode is mult(100icen+j+1).
474	'''
475	SGTextList = []
476	offsetList = []
477	symOpList = []
478	G2oprList = []
479	G2opcodes = []
480	onebar = (1,)
481	if SGData['SGInv']:
482	onebar += (-1,)
483	for icen,cen in enumerate(SGData['SGCen']):
484	for mult in onebar:
485	for j,(M,T) in enumerate(SGData['SGOps']):
486	offset = [0,0,0]
487	Tprime = (mult*T)+cen
488	for i in range(3):
489	while Tprime[i] < 0:
490	Tprime[i] += 1
491	offset[i] += 1
492	while Tprime[i] >= 1:
493	Tprime[i] += -1
494	offset[i] += -1
495	Opr = [mult*M,Tprime]
496	OPtxt = MT2text(Opr)
497	SGTextList.append(OPtxt.replace(' ',''))
498	offsetList.append(tuple(offset))
499	symOpList.append((mult*M,Tprime))
500	G2oprList.append((cen,mult,j))
501	G2opcodes.append(mult(100icen+j+1))
502	return SGTextList,offsetList,symOpList,G2oprList,G2opcodes
503
504	def TextOps(text,table,reverse=False):
505	''' Makes formatted operator list
506	:param text,table: arrays of text made by SGPrint
507	:param reverse: True for x+1/2 form; False for 1/2+x form
508	:returns: OpText: full list of symmetry operators; one operation per line
509	generally printed to console for use via cut/paste in other programs, but
510	could be used for direct input
511	'''
512	OpText = []
513	Inv = True
514	if 'noncentro' in text[1]:
515	Inv = False
516	Cent = [[0,0,0],]
517	if '0,0,0' in text[-1]:
518	Cent = np.array(eval(text[-1].split('+')[0].replace(';','),(')))
519	OpsM = []
520	OpsT = []
521	for item in table:
522	M,T = Text2MT(item.split(')')[1].replace(' ',''),CIF=True)
523	OpsM.append(M)
524	OpsT.append(T)
525	OpsM = np.array(OpsM)
526	OpsT = np.array(OpsT)
527	if Inv:
528	OpsM = np.concatenate((OpsM,-OpsM))
529	OpsT = np.concatenate((OpsT,-OpsT%1.))
530	for cent in Cent:
531	for iop,opM in enumerate(list(OpsM)):
532	txt = MT2text([opM,(OpsT[iop]+cent)%1.],reverse)
533	OpText.append(txt.replace(' ','').lower())
534	return OpText
535
536	def TextGen(SGData,reverse=False):
537	GenSym,GenFlg,BNSsym = GetGenSym(SGData)
538	SGData['GenSym'] = GenSym
539	SGData['GenFlg'] = GenFlg
540	text,table = SGPrint(SGData)
541	GenText = []
542	OprNames = GetOprNames(SGData)
543	OpText = TextOps(text,table,reverse)
544	for name in SGData['GenSym']:
545	gid = OprNames.index(name.replace(' ',''))
546	GenText.append(OpText[gid])
547	if len(SGData['SGCen']) > 1:
548	GenText.append(OpText[-1])
549	return GenText
550
551	def GetOprNames(SGData):
552	OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
553	if SGData['SGInv']:
554	OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
555	return OprNames
556
557	def MT2text(Opr,reverse=False):
558	"From space group matrix/translation operator returns text version"
559	XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
560	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
561	Fld = ''
562	M,T = Opr
563	for j in range(3):
564	IJ = int(round(2M[j][0]+3M[j][1]+4*M[j][2]+4))%12
565	IK = int(round(T[j]*12))%12
566	if IK:
567	if IJ < 3:
568	if reverse:
569	Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
570	else:
571	Fld += (TRA[IK]+XYZ[IJ]).rjust(5)
572	else:
573	if reverse:
574	Fld += (XYZ[IJ]+'+'+TRA[IK]).rjust(5)
575	else:
576	Fld += (TRA[IK]+'+'+XYZ[IJ]).rjust(5)
577	else:
578	Fld += XYZ[IJ].rjust(5)
579	if j != 2: Fld += ', '
580	return Fld
581
582	def Latt2text(Cen):
583	"From lattice centering vectors returns ';' delimited cell centering vectors"
584	lattTxt = ''
585	fracList = ['1/2','1/3','2/3','1/4','3/4','1/5','2/5','3/5','4/5','1/6','5/6',
586	'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']
587	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]
588	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.]
589	nCen = len(Cen)
590	for i,cen in enumerate(Cen):
591	txt = ''
592	for icen in cen:
593	if icen == 1:
594	txt += '1,'
595	continue
596	if not icen:
597	txt += '0,'
598	continue
599	if icen < 0:
600	txt += '-'
601	icen *= -1
602	for mul,prod,frac in zip(mulList,prodList,fracList):
603	if abs(icen*mul-prod) < 1.e-5:
604	txt += frac+','
605	break
606	lattTxt += txt[:-1]+'; '
607	if i and not i%8 and i < nCen-1: #not for the last cen!
608	lattTxt += '\n '
609	return lattTxt[:-2]
610
611	def SpaceGroup(SGSymbol):
612	'''
613	Print the output of SpcGroup in a nicely formatted way.
614
615	:param SGSymbol: space group symbol (string) with spaces between axial fields
616	:returns: nothing
617	'''
618	E,A = SpcGroup(SGSymbol)
619	if E > 0:
620	print (SGErrors(E))
621	return
622	for l in SGPrint(A):
623	print (l)
624	################################################################################
625	#### Magnetic space group stuff
626	################################################################################
627
628	def SetMagnetic(SGData):
629	GenSym,GenFlg,BNSsym = GetGenSym(SGData)
630	SGData['GenSym'] = GenSym
631	SGData['GenFlg'] = GenFlg
632	OprNames,SpnFlp = GenMagOps(SGData)
633	SGData['SpnFlp'] = SpnFlp
634	SGData['MagSpGrp'] = MagSGSym(SGData)
635
636	def GetGenSym(SGData):
637	'''
638	Get the space group generator symbols
639	:param SGData: from :func:`SpcGroup`
640	LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
641	LattSym = ('P','A','B','C','I','F','R')
642
643	'''
644	OprNames = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
645	if SGData['SGInv']:
646	OprNames += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
647	Nsyms = len(SGData['SGOps'])
648	if SGData['SGInv'] and not SGData['SGFixed']: Nsyms *= 2
649	UsymOp = ['1',]
650	OprFlg = [0,]
651	if Nsyms == 2: #Centric triclinic or acentric monoclinic
652	UsymOp.append(OprNames[1])
653	OprFlg.append(SGData['SGGen'][1])
654	elif Nsyms == 4: #Point symmetry 2/m, 222, 22m, or 4
655	if '4z' in OprNames[1]: #Point symmetry 4 or -4
656	UsymOp.append(OprNames[1])
657	OprFlg.append(SGData['SGGen'][1])
658	elif not SGData['SGInv']: #Acentric Orthorhombic
659	if 'm' in OprNames[1:4]: #22m, 2m2 or m22
660	if '2' in OprNames[1]: #Acentric orthorhombic, 2mm
661	UsymOp.append(OprNames[2])
662	OprFlg.append(SGData['SGGen'][2])
663	UsymOp.append(OprNames[3])
664	OprFlg.append(SGData['SGGen'][3])
665	elif '2' in OprNames[2]: #Acentric orthorhombic, m2m
666	UsymOp.append(OprNames[1])
667	OprFlg.append(SGData['SGGen'][1])
668	UsymOp.append(OprNames[3])
669	OprFlg.append(SGData['SGGen'][3])
670	else: #Acentric orthorhombic, mm2
671	UsymOp.append(OprNames[1])
672	OprFlg.append(SGData['SGGen'][1])
673	UsymOp.append(OprNames[2])
674	OprFlg.append(SGData['SGGen'][2])
675	else: #Acentric orthorhombic, 222
676	SGData['SGGen'][1:] = [4,2,1]
677	UsymOp.append(OprNames[1])
678	OprFlg.append(SGData['SGGen'][1])
679	UsymOp.append(OprNames[2])
680	OprFlg.append(SGData['SGGen'][2])
681	UsymOp.append(OprNames[3])
682	OprFlg.append(SGData['SGGen'][3])
683	else: #Centric Monoclinic
684	UsymOp.append(OprNames[1])
685	OprFlg.append(SGData['SGGen'][1])
686	UsymOp.append(OprNames[3])
687	OprFlg.append(SGData['SGGen'][3])
688	elif Nsyms == 6: #Point symmetry 32, 3m or 6
689	if '6' in OprNames[1]: #Hexagonal 6/m Laue symmetry
690	UsymOp.append(OprNames[1])
691	OprFlg.append(SGData['SGGen'][1])
692	else: #Trigonal
693	UsymOp.append(OprNames[4])
694	OprFlg.append(SGData['SGGen'][3])
695	if '2110' in OprNames[1]: UsymOp[-1] = ' 2100 '
696	elif Nsyms == 8: #Point symmetry mmm, 4/m, or 422, etc
697	if '4' in OprNames[1]: #Tetragonal
698	if SGData['SGInv']: #4/m
699	UsymOp.append(OprNames[1])
700	OprFlg.append(SGData['SGGen'][1])
701	UsymOp.append(OprNames[6])
702	OprFlg.append(SGData['SGGen'][6])
703	else:
704	if 'x' in OprNames[4]: #4mm type group
705	UsymOp.append(OprNames[4])
706	OprFlg.append(6)
707	UsymOp.append(OprNames[7])
708	OprFlg.append(8)
709	else: #-42m, -4m2, and 422 type groups
710	UsymOp.append(OprNames[5])
711	OprFlg.append(8)
712	UsymOp.append(OprNames[6])
713	OprFlg.append(19)
714	else: #Orthorhombic, mmm
715	UsymOp.append(OprNames[1])
716	OprFlg.append(SGData['SGGen'][1])
717	UsymOp.append(OprNames[2])
718	OprFlg.append(SGData['SGGen'][2])
719	UsymOp.append(OprNames[7])
720	OprFlg.append(SGData['SGGen'][7])
721	elif Nsyms == 12 and '3' in OprNames[1] and SGData['SGSys'] != 'cubic': #Trigonal
722	UsymOp.append(OprNames[3])
723	OprFlg.append(SGData['SGGen'][3])
724	UsymOp.append(OprNames[9])
725	OprFlg.append(SGData['SGGen'][9])
726	elif Nsyms == 12 and '6' in OprNames[1]: #Hexagonal
727	if 'mz' in OprNames[9]: #6/m
728	UsymOp.append(OprNames[1])
729	OprFlg.append(SGData['SGGen'][1])
730	UsymOp.append(OprNames[6])
731	OprFlg.append(SGData['SGGen'][6])
732	else: #6mm, -62m, -6m2 or 622
733	UsymOp.append(OprNames[6])
734	OprFlg.append(18)
735	if 'm' in UsymOp[-1]: OprFlg[-1] = 20
736	UsymOp.append(OprNames[7])
737	OprFlg.append(24)
738	elif Nsyms in [16,24]:
739	if '3' in OprNames[1]:
740	UsymOp.append('')
741	OprFlg.append(SGData['SGGen'][3])
742	for i in range(Nsyms):
743	if 'mx' in OprNames[i]:
744	UsymOp[-1] = OprNames[i]
745	elif 'm11' in OprNames[i]:
746	UsymOp[-1] = OprNames[i]
747	elif '211' in OprNames[i]:
748	UsymOp[-1] = OprNames[i]
749	OprFlg[-1] = 24
750	else: #4/mmm or 6/mmm
751	UsymOp.append(' mz ')
752	OprFlg.append(1)
753	if '4' in OprNames[1]: #4/mmm
754	UsymOp.append(' mx ')
755	OprFlg.append(20)
756	UsymOp.append(' m110 ')
757	OprFlg.append(24)
758	else: #6/mmm
759	UsymOp.append(' m110 ')
760	OprFlg.append(4)
761	UsymOp.append(' m+-0 ')
762	OprFlg.append(8)
763	else: #System is cubic
764	if Nsyms == 48:
765	UsymOp.append(' mx ')
766	OprFlg.append(4)
767	UsymOp.append(' m110 ')
768	OprFlg.append(24)
769
770	if 'P' in SGData['SGLatt']:
771	if SGData['SGSys'] == 'triclinic':
772	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5]}
773	elif SGData['SGSys'] == 'monoclinic':
774	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5],'P_I':[.5,.5,.5]}
775	if SGData['SGUniq'] == 'a':
776	BNSsym.update({'P_B':[.5,0,.5],'P_C':[.5,.5,0]})
777	elif SGData['SGUniq'] == 'b':
778	BNSsym.update({'P_A':[.5,.5,0],'P_C':[0,.5,.5]})
779	elif SGData['SGUniq'] == 'c':
780	BNSsym.update({'P_A':[0,.5,.5],'P_B':[.5,0,.5]})
781	elif SGData['SGSys'] == 'orthorhombic':
782	BNSsym = {'P_a':[.5,0,0],'P_b':[0,.5,0],'P_c':[0,0,.5],
783	'P_A':[0,.5,.5],'P_B':[.5,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
784	elif SGData['SGSys'] == 'tetragonal':
785	BNSsym = {'P_c':[0,0,.5],'P_C':[.5,.5,0],'P_I':[.5,.5,.5]}
786	elif SGData['SGSys'] in ['trigonal','hexagonal']:
787	BNSsym = {'P_c':[0,0,.5]}
788	elif SGData['SGSys'] == 'cubic':
789	BNSsym = {'P_I':[.5,.5,.5]}
790
791	elif 'A' in SGData['SGLatt']:
792	if SGData['SGSys'] == 'monoclinic':
793	BNSsym = {}
794	if SGData['SGUniq'] == 'b':
795	BNSsym.update({'A_a':[.5,0,0],'A_c':[0,0,.5]})
796	elif SGData['SGUniq'] == 'c':
797	BNSsym.update({'A_a':[.5,0,0],'A_b':[0,.5,0]})
798	elif SGData['SGSys'] == 'orthorhombic':
799	BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5],
800	'A_B':[.5,0,.5],'A_C':[.5,.5,0]}
801	elif SGData['SGSys'] == 'triclinic':
802	BNSsym = {'A_a':[.5,0,0],'A_b':[0,.5,0],'A_c':[0,0,.5]}
803
804	elif 'B' in SGData['SGLatt']:
805	if SGData['SGSys'] == 'monoclinic':
806	BNSsym = {}
807	if SGData['SGUniq'] == 'a':
808	BNSsym.update({'B_b':[0,.5,0],'B_c':[0,0,.5]})
809	elif SGData['SGUniq'] == 'c':
810	BNSsym.update({'B_a':[.5,0,0],'B_b':[0,.5,0]})
811	elif SGData['SGSys'] == 'orthorhombic':
812	BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5],
813	'B_A':[0,.5,.5],'B_C':[.5,.5,0]}
814	elif SGData['SGSys'] == 'triclinic':
815	BNSsym = {'B_a':[.5,0,0],'B_b':[0,.5,0],'B_c':[0,0,.5]}
816
817	elif 'C' in SGData['SGLatt']:
818	if SGData['SGSys'] == 'monoclinic':
819	BNSsym = {}
820	if SGData['SGUniq'] == 'a':
821	BNSsym.update({'C_b':[0,.5,.0],'C_c':[0,0,.5]})
822	elif SGData['SGUniq'] == 'b':
823	BNSsym.update({'C_a':[.5,0,0],'C_c':[0,0,.5],'C_B':[.5,0.,.5]})
824	elif SGData['SGSys'] == 'orthorhombic':
825	BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5],
826	'C_A':[0,.5,.5],'C_B':[.5,0,.5]}
827	elif SGData['SGSys'] == 'triclinic':
828	BNSsym = {'C_a':[.5,0,0],'C_b':[0,.5,0],'C_c':[0,0,.5]}
829
830	elif 'I' in SGData['SGLatt']:
831	if SGData['SGSys'] in ['monoclinic','orthorhombic','triclinic']:
832	BNSsym = {'I_a':[.5,0,0],'I_b':[0,.5,0],'I_c':[0,0,.5]}
833	elif SGData['SGSys'] == 'tetragonal':
834	BNSsym = {'I_c':[0,0,.5]}
835	elif SGData['SGSys'] == 'cubic':
836	BNSsym = {}
837
838	elif 'F' in SGData['SGLatt']:
839	if SGData['SGSys'] in ['monoclinic','orthorhombic','cubic','triclinic']:
840	BNSsym = {'F_S':[.5,.5,.5]}
841
842	elif 'R' in SGData['SGLatt']:
843	BNSsym = {'R_I':[0,0,.5]}
844	return UsymOp,OprFlg,BNSsym
845
846	def ApplyBNSlatt(SGData,BNSlatt):
847	Tmat = np.eye(3)
848	BNS = BNSlatt[0]
849	A = np.array(BNSlatt[1])
850	Laue = SGData['SGLaue']
851	SGCen = SGData['SGCen']
852	if '_a' in BNS:
853	Tmat[0,0] = 2.0
854	elif '_b' in BNS:
855	Tmat[1,1] = 2.0
856	elif '_c' in BNS:
857	Tmat[2,2] = 2.0
858	elif '_A' in BNS:
859	Tmat[0,0] = 2.0
860	elif '_B' in BNS:
861	Tmat[1,1] = 2.0
862	elif '_C' in BNS:
863	Tmat[2,2] = 2.0
864	elif '_I' in BNS:
865	Tmat *= 2.0
866	if 'R' in Laue:
867	SGData['SGSpin'][-1] = -1
868	else:
869	SGData['SGSpin'].append(-1)
870	elif '_S' in BNS:
871	SGData['SGSpin'][-1] = -1
872	SGData['SGSpin'] += [-1,-1,-1,]
873	Tmat *= 2.0
874	else:
875	return Tmat
876	SGData['SGSpin'].append(-1)
877	C = SGCen+A
878	SGData['SGCen'] = np.vstack((SGCen,C))%1.
879	return Tmat
880
881	def CheckSpin(isym,SGData):
882	''' Check for exceptions in spin rules
883	'''
884	if SGData['SGPtGrp'] in ['222','mm2','2mm','m2m']: #only 2/3 can be red; not 1/3 or 3/3
885	if SGData['SGSpin'][1]SGData['SGSpin'][2]SGData['SGSpin'][3] < 0:
886	SGData['SGSpin'][(isym+1)%3+1] *= -1
887	if SGData['SpGrp'][0] == 'F' and isym > 2:
888	SGData['SGSpin'][(isym+1)%3+3] == 1
889	elif SGData['SGPtGrp'] == 'mmm':
890	if SGData['SpGrp'][0] == 'F' and isym > 2:
891	SGData['SGSpin'][(isym+1)%3+3] == 1
892
893	def MagSGSym(SGData): #needs to use SGPtGrp not SGLaue!
894	SGLaue = SGData['SGLaue']
895	if '1' not in SGData['GenSym']: #patch for old gpx files
896	SGData['GenSym'] = ['1',]+SGData['GenSym']
897	SGData['SGSpin'] = [1,]+list(SGData['SGSpin'])
898	if len(SGData['SGSpin'])<len(SGData['GenSym']):
899	SGData['SGSpin'] = [1,]+list(SGData['SGSpin']) #end patch
900	GenSym = SGData['GenSym'][1:] #skip identity
901	SpnFlp = SGData['SGSpin']
902	# print('SpnFlp',SpnFlp)
903	SGPtGrp = SGData['SGPtGrp']
904	if len(SpnFlp) == 1:
905	SGData['MagPtGp'] = SGPtGrp
906	return SGData['SpGrp']
907	magSym = SGData['SpGrp'].split()
908	if SGLaue in ['-1',]:
909	SGData['MagPtGp'] = SGPtGrp
910	if SpnFlp[1] == -1:
911	magSym[1] += "'"
912	SGData['MagPtGp'] += "'"
913	elif SGLaue in ['2/m','4/m','6/m']: #all ok
914	Uniq = {'a':1,'b':2,'c':3,'':1}
915	id = [0,1]
916	if len(magSym) > 2:
917	id = [0,Uniq[SGData['SGUniq']]]
918	sym = magSym[id[1]].split('/')
919	Ptsym = SGLaue.split('/')
920	if len(GenSym) == 3:
921	for i in [0,1,2]:
922	if SpnFlp[i+1] < 0:
923	sym[i] += "'"
924	Ptsym[i] += "'"
925	else:
926	for i in range(len(GenSym)):
927	if SpnFlp[i+1] < 0:
928	sym[i] += "'"
929	Ptsym[i] += "'"
930	SGData['MagPtGp'] = '/'.join(Ptsym)
931	magSym[id[1]] = '/'.join(sym)
932	elif SGPtGrp in ['mmm','mm2','m2m','2mm','222']:
933	SGData['MagPtGp'] = ''
934	for i in [0,1,2]:
935	SGData['MagPtGp'] += SGPtGrp[i]
936	if SpnFlp[i+1] < 0:
937	magSym[i+1] += "'"
938	SGData['MagPtGp'] += "'"
939	elif SGLaue == '6/mmm': #ok
940	magPtGp = list(SGPtGrp)
941	if len(GenSym) == 2:
942	for i in [0,1]:
943	if SpnFlp[i+1] < 0:
944	magSym[i+2] += "'"
945	magPtGp[i+1] += "'"
946	if SpnFlp[1]*SpnFlp[2] < 0:
947	magSym[1] += "'"
948	magPtGp[0] += "'"
949	else:
950	sym = magSym[1].split('/')
951	Ptsym = ['6','m']
952	magPtGp = ['','m','m']
953	for i in [0,1,2]:
954	if SpnFlp[i+1] < 0:
955	if i:
956	magSym[i+1] += "'"
957	magPtGp[i] += "'"
958	else:
959	sym[1] += "'"
960	Ptsym[0] += "'"
961	if SpnFlp[2]*SpnFlp[3] < 0:
962	sym[0] += "'"
963	Ptsym[0] += "'"
964	magSym[1] = '/'.join(sym)
965	magPtGp[0] = '/'.join(Ptsym)
966	SGData['MagPtGp'] = ''.join(magPtGp)
967	elif SGLaue == '4/mmm':
968	magPtGp = list(SGPtGrp)
969	if len(GenSym) == 2:
970	for i in [0,1]:
971	if SpnFlp[i+1] < 0:
972	magSym[i+2] += "'"
973	magPtGp[i+1] += "'"
974	if SpnFlp[1]*SpnFlp[2] < 0:
975	magSym[1] += "'"
976	magPtGp[0] += "'"
977	else:
978	if '/' in magSym[1]: #P 4/m m m, etc.
979	sym = magSym[1].split('/')
980	Ptsym = ['4','m']
981	magPtGp = ['','m','m']
982	for i in [0,1,2]:
983	if SpnFlp[i+1] < 0:
984	if i:
985	magSym[i+1] += "'"
986	magPtGp[i] += "'"
987	else:
988	sym[1] += "'"
989	Ptsym[1] += "'"
990	if SpnFlp[2]*SpnFlp[3] < 0:
991	sym[0] += "'"
992	Ptsym[0] += "'"
993	magSym[1] = '/'.join(sym)
994	magPtGp[0] = '/'.join(Ptsym)
995	else:
996	for i in [0,1]:
997	if SpnFlp[i+1] < 0:
998	magSym[i+2] += "'"
999	if SpnFlp[1]*SpnFlp[2] < 0:
1000	magSym[1] += "'"
1001	SGData['MagPtGp'] = ''.join(magPtGp)
1002	elif SGLaue in ['3','3m1','31m']: #ok
1003	Ptsym = list(SGLaue)
1004	if len(GenSym) == 1: #all ok
1005	id = 2
1006	if (len(magSym) == 4) and (magSym[2] == '1'):
1007	id = 3
1008	if '3' in GenSym[0]:
1009	id = 1
1010	magSym[id].strip("'")
1011	if SpnFlp[1] < 0:
1012	magSym[id] += "'"
1013	Ptsym[id-1] += "'"
1014	elif len(GenSym) == 2:
1015	if 'R' in GenSym[1]:
1016	magSym[-1].strip("'")
1017	if SpnFlp[1] < 0:
1018	magSym[-1] += "'"
1019	Ptsym[-1] += "'"
1020	else:
1021	i,j = [1,2]
1022	if magSym[2] == '1':
1023	i,j = [1,3]
1024	magSym[i].strip("'")
1025	Ptsym[i-1].strip("'")
1026	magSym[j].strip("'")
1027	Ptsym[j-1].strip("'")
1028	if SpnFlp[1:3] == [1,-1]:
1029	magSym[i] += "'"
1030	Ptsym[i-1] += "'"
1031	elif SpnFlp[1:3] == [-1,-1]:
1032	magSym[j] += "'"
1033	Ptsym[j-1] += "'"
1034	elif SpnFlp[1:3] == [-1,1]:
1035	magSym[i] += "'"
1036	Ptsym[i-1] += "'"
1037	magSym[j] += "'"
1038	Ptsym[j-1] += "'"
1039	elif len(GenSym):
1040	if 'c' not in magSym[2]:
1041	i,j = [1,2]
1042	magSym[i].strip("'")
1043	Ptsym[i-1].strip("'")
1044	magSym[j].strip("'")
1045	Ptsym[j-1].strip("'")
1046	if SpnFlp[1:3] == [1,-1]:
1047	magSym[i] += "'"
1048	Ptsym[i-1] += "'"
1049	elif SpnFlp[1:3] == [-1,-1]:
1050	magSym[j] += "'"
1051	Ptsym[j-1] += "'"
1052	elif SpnFlp[2] == [-1,1]:
1053	magSym[i] += "'"
1054	Ptsym[i-1] += "'"
1055	magSym[j] += "'"
1056	Ptsym[j-1] += "'"
1057	SGData['MagPtGp'] = ''.join(Ptsym)
1058	elif SGData['SGPtGrp'] == '23' and len(magSym):
1059	SGData['MagPtGp'] = '23'
1060	elif SGData['SGPtGrp'] == 'm3':
1061	SGData['MagPtGp'] = "m3"
1062	if SpnFlp[1] < 0:
1063	magSym[1] += "'"
1064	magSym[2] += "'"
1065	SGData['MagPtGp'] = "m'3'"
1066	if SpnFlp[2] < 0:
1067	if not 'm' in magSym[1]: #only Ia3
1068	magSym[1].strip("'")
1069	SGData['MagPtGp'] = "m3'"
1070	elif SGData['SGPtGrp'] in ['432','-43m']:
1071	Ptsym = SGData['SGPtGrp'].split('3')
1072	if SpnFlp[1] < 0:
1073	magSym[1] += "'"
1074	Ptsym[0] += "'"
1075	magSym[3] += "'"
1076	Ptsym[1] += "'"
1077	SGData['MagPtGp'] = '3'.join(Ptsym)
1078	elif SGData['SGPtGrp'] == 'm3m':
1079	Ptsym = ['m','3','m']
1080	if SpnFlp[1:3] == [-1,1]:
1081	magSym[1] += "'"
1082	Ptsym[0] += "'"
1083	magSym[2] += "'"
1084	Ptsym[1] += "'"
1085	elif SpnFlp[1:3] == [1,-1]:
1086	magSym[3] += "'"
1087	Ptsym[2] += "'"
1088	elif SpnFlp[1:3] == [-1,-1]:
1089	magSym[1] += "'"
1090	Ptsym[0] += "'"
1091	magSym[2] += "'"
1092	Ptsym[1] += "'"
1093	magSym[3] += "'"
1094	Ptsym[2] += "'"
1095	SGData['MagPtGp'] = ''.join(Ptsym)
1096	# print SpnFlp
1097	magSym[0] = SGData.get('BNSlattsym',[SGData['SGLatt'],[0,0,0]])[0]
1098	return ' '.join(magSym)
1099
1100	def fixMono(SpGrp):
1101	'fixes b-unique monoclinics in e.g. P 1 2/1c 1 --> P 21/c '
1102	Flds = SpGrp.split()
1103	if len(Flds) == 4:
1104	if Flds[2] != '1':
1105	return '%s %s'%(Flds[0],Flds[2])
1106	else:
1107	return None
1108	else:
1109	return SpGrp
1110
1111	def Trans2Text(Trans):
1112	"from transformation matrix to text"
1113	cells = ['a','b','c']
1114	Text = ''
1115	for row in Trans:
1116	Fld = ''
1117	for i in [0,1,2]:
1118	if row[i]:
1119	if Fld and row[i] > 0.:
1120	Fld += '+'
1121	Fld += '%3.1f'%(row[i])+cells[i]
1122	Text += Fld
1123	Text += ','
1124	Text = Text.replace('1.0','').replace('.0','').replace('0.5','1/2')
1125	return Text[:-1]
1126
1127	def getlattSym(Trans):
1128	Fives = {'ababc':'abc','bcbca':'cba','acacb':'acb','cabab':'cab','abcab':'acb'}
1129	transText = Trans2Text(Trans)
1130	lattSym = ''
1131	for fld in transText.split(','):
1132	if 'a' in fld: lattSym += 'a'
1133	if 'b' in fld: lattSym += 'b'
1134	if 'c' in fld: lattSym += 'c'
1135	if len(lattSym) != 3:
1136	lattSym = 'abc'
1137	# lattSym = Fives[lattSym]
1138	return lattSym
1139
1140
1141	def Text2MT(mcifOpr,CIF=True):
1142	"From space group cif text returns matrix/translation"
1143	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],
1144	'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],
1145	'+x-y':[1,-1,0],'+y-x':[-1,1,0]}
1146	ops = mcifOpr.split(",")
1147	M = []
1148	T = []
1149	for op in ops[:3]:
1150	ip = len(op)
1151	if '/' in op:
1152	try: #mcif format
1153	nP = op.count('+')
1154	opMT = op.split('+')
1155	T.append(eval(opMT[nP]))
1156	if nP == 2:
1157	opMT[0] = '+'.join(opMT[0:2])
1158	except NameError: #normal cif format
1159	ip = op.index('/')
1160	T.append(eval(op[:ip+2]))
1161	opMT = [op[ip+2:],'']
1162	else:
1163	opMT = [op,'']
1164	T.append(0.)
1165	M.append(XYZ[opMT[0].lower()])
1166	return np.array(M),np.array(T)
1167
1168	def MagText2MTS(mcifOpr,CIF=True):
1169	"From magnetic space group cif text returns matrix/translation + spin flip"
1170	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],
1171	'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],
1172	'+x-y':[1,-1,0],'+y-x':[-1,1,0]}
1173	ops = mcifOpr.split(",")
1174	M = []
1175	T = []
1176	for op in ops[:3]:
1177	ip = len(op)
1178	if '/' in op:
1179	try: #mcif format
1180	nP = op.count('+')
1181	opMT = op.split('+')
1182	T.append(eval(opMT[nP]))
1183	if nP == 2:
1184	opMT[0] = '+'.join(opMT[0:2])
1185	except NameError: #normal cif format
1186	ip = op.index('/')
1187	T.append(eval(op[:ip+2]))
1188	opMT = [op[ip+2:],'']
1189	else:
1190	opMT = [op,'']
1191	T.append(0.)
1192	M.append(XYZ[opMT[0].lower()])
1193	spnflp = 1
1194	if '-1' in ops[3]:
1195	spnflp = -1
1196	return np.array(M),np.array(T),spnflp
1197
1198	def MagSSText2MTS(mcifOpr):
1199	"From magnetic super space group cif text returns matrix/translation + spin flip"
1200	XYZ = {'x1':[1,0,0,0],'-x1':[-1,0,0,0],
1201	'x2':[0,1,0,0],'-x2':[0,-1,0,0],
1202	'x3':[0,0,1,0],'-x3':[0,0,-1,0],
1203	'x4':[0,0,0,1],'-x4':[0,0,0,-1],
1204	'x1-x2':[1,-1,0,0],'-x1+x2':[-1,1,0,0],
1205	'x1-x4':[1,0,0,-1],'-x1+x4':[-1,0,0,1],
1206	'x2-x4':[0,1,0,-1],'-x2+x4':[0,-1,0,1],
1207	'-x1-x2+x4':[-1,-1,0,1],'x1+x2-x4':[1,1,0,-1]}
1208	ops = mcifOpr.split(",")
1209	M = []
1210	T = []
1211	for op in ops[:4]:
1212	ip = len(op)
1213	if '/' in op:
1214	ip = op.index('/')-2
1215	T.append(eval(op[ip:]))
1216	else:
1217	T.append(0.)
1218	M.append(XYZ[op[:ip]])
1219	spnflp = 1
1220	if '-1' in ops[4]:
1221	spnflp = -1
1222	return np.array(M),np.array(T),spnflp
1223
1224	def GetSGSpin(SGData,MSgSym):
1225	'get spin generators from magnetic space group symbol'
1226	SpGrp = SGData['SpGrp']
1227	mSgSym = MSgSym+' '
1228	Flds = SpGrp.split()
1229	iB = 0
1230	Spn = [1,] #for identity generator
1231	if len(Flds) == 2: #-1, 2/m, 4/m & 6/m; 1 or 2 generators
1232	fld = Flds[1]
1233	iF = mSgSym[iB:].index(fld[0])+iB
1234	jF = mSgSym[iF:].index(fld[-1])+iF
1235	if '/' in mSgSym[iF:jF]:
1236	if "'" in mSgSym[iF:jF]:
1237	Spn.append(-1)
1238	else:
1239	Spn.append(1)
1240	if "'" == mSgSym[jF+1]:
1241	Spn.append(-1)
1242	else:
1243	Spn.append(1)
1244	elif len(Flds) == 3: # 3m & m3; 1 or 2 generator
1245	if SGData['SGPtGrp'] == '-3m':
1246	if not mSgSym.count("'"):
1247	Spn += [1,1,]
1248	elif mSgSym.count("'") == 2:
1249	Spn += [-1,1,]
1250	elif "3'" in mSgSym:
1251	Spn += [1,-1,]
1252	else:
1253	Spn += [-1,-1,]
1254	else:
1255	if "'" in mSgSym: #could be 1 or 2 '; doesn't matter.
1256	Spn.append(-1)
1257	else:
1258	Spn.append(1)
1259	else: #the rest; 3 generators. NB: any ' before / in 1st field ignored
1260	for fld in Flds[1:]:
1261	iF = mSgSym[iB:].index(fld[0])+iB
1262	jF = mSgSym[iF:].index(fld[-1])+iF
1263	if "'" == mSgSym[jF+1]:
1264	Spn.append(-1)
1265	iB = jF+2
1266	else:
1267	Spn.append(1)
1268	iB = jF+1
1269	return Spn
1270
1271	def GenMagOps(SGData):
1272	FlpSpn = SGData['SGSpin']
1273	Nsym = len(SGData['SGOps'])
1274	Ncv = len(SGData['SGCen'])
1275	sgOp = [M for M,T in SGData['SGOps']]
1276	oprName = [GetOprPtrName(str(irtx)) for irtx in PackRot(SGData['SGOps'])]
1277	if SGData['SGInv'] and not SGData['SGFixed'] and not SGData['SGGray']:
1278	Nsym *= 2
1279	sgOp += [-M for M,T in SGData['SGOps']]
1280	oprName += [GetOprPtrName(str(-irtx)) for irtx in PackRot(SGData['SGOps'])]
1281	Nsyms = 0
1282	sgOps = []
1283	OprNames = []
1284	for incv in range(Ncv):
1285	Nsyms += Nsym
1286	sgOps += sgOp
1287	OprNames += oprName
1288	if SGData['SGFixed']:
1289	SpnFlp = SGData['SpnFlp']
1290	else:
1291	SpnFlp = np.ones(Nsym,dtype=np.int)
1292	GenFlg = SGData.get('GenFlg',[0])
1293	Ngen = len(SGData['SGGen'])
1294	# print ('GenFlg:',SGData['GenFlg'])
1295	# print ('GenSym:',SGData['GenSym'])
1296	Nfl = len(GenFlg)
1297	for ieqv in range(Nsym):
1298	for iunq in range(Nfl):
1299	if SGData['SGGen'][ieqv%Ngen] & GenFlg[iunq]:
1300	SpnFlp[ieqv] *= FlpSpn[iunq]
1301	# print ('\nMagSpGrp:',SGData['MagSpGrp'],Ncv)
1302	# print ('FlpSpn:',Nfl,FlpSpn)
1303	for incv in range(Ncv):
1304	if incv:
1305	try:
1306	SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))
1307	except IndexError:
1308	FlpSpn = [1,]+FlpSpn
1309	SpnFlp = np.concatenate((SpnFlp,SpnFlp[:Nsym]*FlpSpn[Nfl+incv-1]))
1310	detM = [nl.det(M) for M in sgOp]
1311	MagMom = SpnFlpnp.array(NcvdetM) #duplicate for no. centerings
1312	SGData['MagMom'] = MagMom
1313	# print ('SgOps:',OprNames)
1314	# print ('SGGen:',SGData['SGGen'])
1315	# print ('SpnFlp:',SpnFlp)
1316	# print ('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 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,8,9]):
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	if SGData['SGGray'] and SSymbol[-1] == 's':
1641	SSymbol = SSymbol[:-1]
1642	try:
1643	modsym,gensym = splitSSsym(SSymbol)
1644	except ValueError:
1645	return 'Error in superspace symbol '+SSymbol,None
1646	modQ = [Fracs[mod] for mod in modsym]
1647	SSGKl = SGData['SSGKl'][:]
1648	if SGData['SGLaue'] in ['2/m','mmm']:
1649	SSGKl = fixMonoOrtho()
1650	Ngen = len(gensym)
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	SSgSpc += 's'
1661	SSGData = {'SSpGrp':SSgSpc,'modQ':modQ,'modSymb':modsym,'SSGKl':SSGKl}
1662	SSCen = np.zeros((len(SGData['SGCen']),4))
1663	for icen,cen in enumerate(SGData['SGCen']):
1664	SSCen[icen,0:3] = cen
1665	SSCen[0] = np.zeros(4)
1666	SSGData['SSGCen'] = SSCen
1667	SSGData['SSGOps'] = []
1668	for iop,op in enumerate(SGData['SGOps']):
1669	T = np.zeros(4)
1670	ssop = np.zeros((4,4))
1671	ssop[:3,:3] = op[0]
1672	T[:3] = op[1]
1673	SSGData['SSGOps'].append([ssop,T])
1674	E,Result = genSSGOps()
1675	if E:
1676	SSGData['SSGOps'] = Result
1677	if DEBUG:
1678	print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1679	for Op in Result:
1680	print (SSMT2text(Op).replace(' ',''))
1681	if SGData['SGInv']:
1682	for Op in Result:
1683	Op = [-Op[0],-Op[1]%1.]
1684	print (SSMT2text(Op).replace(' ',''))
1685	return None,SSGData
1686	else:
1687	return Result+'\nOperator conflict - incorrect superspace symbol',None
1688
1689	def SSChoice(SGData):
1690	'''
1691	Gets the unique set of possible super space groups for a given space group
1692	'''
1693	laueSS = {'-1':['(abg)',],
1694	'2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1695	'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1696	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1697	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1698	'4/m':['(00g)','(1/21/2g)'],
1699	'4/mmm':['(00g)','(1/21/2g)'],
1700	'3':['(00g)','(1/31/3g)'],
1701	'3m1':['(00g)'],
1702	'31m':['(00g)','(1/31/3g)'],
1703	'6/m':['(00g)',],
1704	'6/mmm':['(00g)',]}
1705
1706	laueTS = {'-1':['',],
1707	'2/m':['','s','0s','ss','s0'],
1708	'mmm':['','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1709	'4/m':['','q','s','s0',],
1710	'4/mmm':['','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1711	'3':['','t'],
1712	'3m1':['','t0','0s','t00','0s0'],
1713	'31m':['','t00','0ss'],
1714	'6/m':['','h','t','s','s0'],
1715	'6/mmm':['','h00','t00','s00','ss0','0ss','s0s','s0s0','00ss','s00s','ss00','0ss0','0s0s']}
1716	laue = SGData['SGLaue']
1717	SSChoice = []
1718	for ax in laueSS[laue]:
1719	for sx in laueTS[laue]:
1720	SSChoice.append(ax+sx)
1721	ssChoice = []
1722	ssHash = []
1723	for item in SSChoice:
1724	E,SSG = SSpcGroup(SGData,item)
1725	if SSG:
1726	sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1727	if sshash not in ssHash:
1728	ssHash.append(sshash)
1729	ssChoice.append(item)
1730	return ssChoice
1731
1732	def fixGray(SGData,SSymbol):
1733	modsym,gensym = SSymbol.replace(' ','').split(')')
1734	modsym += ')'
1735	sgPtGp = SGData['SGPtGrp']
1736	if gensym:
1737	if sgPtGp in ['1','2','m','3','4','6'] and len(gensym) == 1:
1738	gensym += 's'
1739	elif sgPtGp in ['2/m','4/m','6/m'] and len(gensym) == 2:
1740	gensym += 's'
1741	elif sgPtGp in ['4/mmm','6/mmm'] and len(gensym) == 4:
1742	gensym += 's'
1743	elif len(gensym) == 3:
1744	gensym += 's'
1745	else:
1746	if sgPtGp in ['1','2','m','3','4','6']:
1747	gensym += '0s'
1748	elif sgPtGp in ['2/m','4/m','6/m']:
1749	gensym += '00s'
1750	elif sgPtGp in ['4/mmm','6/mmm']:
1751	gensym += '0000s'
1752	else:
1753	gensym += '000s'
1754	return modsym+gensym
1755
1756	def splitSSsym(SSymbol):
1757	'''
1758	Splits supersymmetry symbol into two lists of strings
1759	'''
1760	modsym,gensym = SSymbol.replace(' ','').split(')')
1761	modsym = modsym.replace(',','')
1762	if "1'" in modsym:
1763	gensym = gensym[:-1]
1764	modsym = modsym.replace("1'",'')
1765	if gensym in ['0','00','000','0000']: #get rid of extraneous symbols
1766	gensym = ''
1767	nfrac = modsym.count('/')
1768	modsym = modsym.lstrip('(')
1769	if nfrac == 0:
1770	modsym = list(modsym)
1771	elif nfrac == 1:
1772	pos = modsym.find('/')
1773	if pos == 1:
1774	modsym = [modsym[:3],modsym[3],modsym[4]]
1775	elif pos == 2:
1776	modsym = [modsym[0],modsym[1:4],modsym[4]]
1777	else:
1778	modsym = [modsym[0],modsym[1],modsym[2:]]
1779	else:
1780	lpos = modsym.find('/')
1781	rpos = modsym.rfind('/')
1782	if lpos == 1 and rpos == 4:
1783	modsym = [modsym[:3],modsym[3:6],modsym[6]]
1784	elif lpos == 1 and rpos == 5:
1785	modsym = [modsym[:3],modsym[3],modsym[4:]]
1786	else:
1787	modsym = [modsym[0],modsym[1:4],modsym[4:]]
1788	gensym = list(gensym)
1789	return modsym,gensym
1790
1791	def SSGPrint(SGData,SSGData):
1792	'''
1793	Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1794
1795	:param SGData: space group data structure as defined in SpcGroup above.
1796	:param SSGData: from :func:`SSpcGroup`
1797	:returns:
1798	SSGText - list of strings with the superspace group details
1799	SGTable - list of strings for each of the operations
1800	'''
1801	Mult = len(SSGData['SSGCen'])len(SSGData['SSGOps'])(int(SGData['SGInv'])+1)
1802	SSGText = []
1803	SSGText.append(' Superspace Group: '+SSGData['SSpGrp'])
1804	CentStr = 'centrosymmetric'
1805	if not SGData['SGInv']:
1806	CentStr = 'non'+CentStr
1807	if SGData['SGLatt'] in 'ABCIFR':
1808	SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1809	else:
1810	SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1811	SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1812	SGptGp = SGData['SGPtGrp']
1813	if SGData['SGGray']:
1814	SGptGp += "1'"
1815	SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1816	SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1817	SSGText.append(' Multiplicity of a general site is '+str(Mult))
1818	if SGData['SGUniq'] in ['a','b','c']:
1819	SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1820	if SGData['SGInv']:
1821	SSGText.append(' The inversion center is located at 0,0,0')
1822	if SGData['SGPolax']:
1823	SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1824	SSGText.append(' ')
1825	if len(SSGData['SSGCen']) > 1:
1826	SSGText.append(' The equivalent positions are:')
1827	SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
1828	else:
1829	SSGText.append(' The equivalent positions are:\n')
1830	SSGTable = []
1831	for i,Opr in enumerate(SSGData['SSGOps']):
1832	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1833	if SGData['SGGray']:
1834	SSGTable.append(" for 1'")
1835	for i,Opr in enumerate(SSGData['SSGOps']):
1836	Opr2 = [Opr[0],Opr[1]+np.array([0,0,0,.5])]
1837	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr2)))
1838	return SSGText,SSGTable
1839
1840	def SSGModCheck(Vec,modSymb,newMod=True):
1841	''' Checks modulation vector compatibility with supersymmetry space group symbol.
1842	if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1843	'''
1844	Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1845	modQ = [Fracs[mod] for mod in modSymb]
1846	if newMod:
1847	newVec = Vec
1848	if not np.any(Vec):
1849	newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1850	return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)], \
1851	[True if mod in ['a','b','g'] else False for mod in modSymb]
1852	else:
1853	return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1854
1855	def SSMT2text(Opr):
1856	"From superspace group matrix/translation operator returns text version"
1857	XYZS = ('x','y','z','t') #Stokes, Campbell & van Smaalen notation
1858	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1859	Fld = ''
1860	M,T = Opr
1861	for j in range(4):
1862	IJ = ''
1863	for k in range(4):
1864	txt = str(int(round(M[j][k])))
1865	txt = txt.replace('1',XYZS[k]).replace('0','')
1866	if '2' in txt:
1867	txt += XYZS[k]
1868	if IJ and M[j][k] > 0:
1869	IJ += '+'+txt
1870	else:
1871	IJ += txt
1872	IK = int(round(T[j]*12))%12
1873	if IK:
1874	if not IJ:
1875	break
1876	if IJ[0] == '-':
1877	Fld += (TRA[IK]+IJ).rjust(8)
1878	else:
1879	Fld += (TRA[IK]+'+'+IJ).rjust(8)
1880	else:
1881	Fld += IJ.rjust(8)
1882	if j != 3: Fld += ', '
1883	return Fld
1884
1885	def SSLatt2text(SSGCen):
1886	"Lattice centering vectors to text"
1887	lattTxt = ''
1888	lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1889	for vec in SSGCen:
1890	lattTxt += ' '
1891	for item in vec:
1892	lattTxt += '%s,'%(lattDir[int(item*12)])
1893	lattTxt = lattTxt.rstrip(',')
1894	lattTxt += ';'
1895	lattTxt = lattTxt.rstrip(';').lstrip(' ')
1896	return lattTxt
1897
1898	def SSpaceGroup(SGSymbol,SSymbol):
1899	'''
1900	Print the output of SSpcGroup in a nicely formatted way.
1901
1902	:param SGSymbol: space group symbol with spaces between axial fields.
1903	:param SSymbol: superspace group symbol extension (string).
1904	:returns: nothing
1905	'''
1906
1907	E,A = SpcGroup(SGSymbol)
1908	if E > 0:
1909	print (SGErrors(E))
1910	return
1911	E,B = SSpcGroup(A,SSymbol)
1912	if E > 0:
1913	print (E)
1914	return
1915	for l in SSGPrint(B):
1916	print (l)
1917
1918	def SGProd(OpA,OpB):
1919	'''
1920	Form space group operator product. OpA & OpB are [M,V] pairs;
1921	both must be of same dimension (3 or 4). Returns [M,V] pair
1922	'''
1923	A,U = OpA
1924	B,V = OpB
1925	M = np.inner(B,A.T)
1926	W = np.inner(B,U)+V
1927	return M,W
1928
1929	def GetLittleGrpOps(SGData,vec):
1930	''' Find rotation part of operators that leave vec unchanged
1931
1932	:param SGData: space group data structure as defined in SpcGroup above.
1933	:param vec: a numpy array of fractional vector coordinates
1934	:returns: Little - list of operators [M,T] that form the little gropu
1935	'''
1936	Little = []
1937	Ops = SGData['SGOps'][:]
1938	if SGData['SGInv']:
1939	Ops += [[-M,-T] for [M,T] in Ops]
1940	for [M,T] in Ops:
1941	tvec = np.inner(M,vec)%1.
1942	if np.allclose(tvec,vec%1.):
1943	Little.append([M,T])
1944	return Little
1945
1946	def MoveToUnitCell(xyz):
1947	'''
1948	Translates a set of coordinates so that all values are >=0 and < 1
1949
1950	:param xyz: a list or numpy array of fractional coordinates
1951	:returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
1952	'''
1953	XYZ = (np.array(xyz)+10.)%1.
1954	cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
1955	return XYZ,cell
1956
1957	def Opposite(XYZ,toler=0.0002):
1958	'''
1959	Gives opposite corner, edge or face of unit cell for position within tolerance.
1960	Result may be just outside the cell within tolerance
1961
1962	:param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1963	:param toler: unit cell fraction tolerance making opposite
1964	:returns:
1965	XYZ: dict of opposite positions; key=unit cell & always contains XYZ
1966	'''
1967	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]]
1968	TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1969	perm = TB*perm3
1970	cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
1971	D = dict(zip(cperm,perm))
1972	new = {}
1973	for key in D:
1974	new[key] = np.array(D[key])+np.array(XYZ)
1975	return new
1976
1977	def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1978	'''
1979	Generates the equivalent positions for a specified coordinate and space group
1980
1981	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
1982	:param SGData: from :func:`SpcGroup`
1983	:param All: True return all equivalent positions including duplicates;
1984	False return only unique positions
1985	:param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1986	:param Move: True move generated atom positions to be inside cell
1987	False do not move atoms
1988	:return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
1989
1990	* [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1991	* Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1992	* is centering operator number (1-4) and - is for inversion
1993	Cell = unit cell translations needed to put new positions inside cell
1994	[UijEquiv] - equivalent Uij; absent if no Uij given
1995	* +1/-1 for spin inversion of operator - empty if not magnetic
1996
1997	'''
1998	XYZEquiv = []
1999	UijEquiv = []
2000	Idup = []
2001	Cell = []
2002	inv = int(SGData['SGInv']+1)
2003	icen = SGData['SGCen']
2004	if SGData['SGFixed']:
2005	inv = 1
2006	SpnFlp = SGData.get('SpnFlp',[])
2007	spnflp = []
2008	X = np.array(XYZ)
2009	mj = 0
2010	cell0 = np.zeros(3,dtype=np.int32)
2011	if Move:
2012	X,cell0 = MoveToUnitCell(X)
2013	for ic,cen in enumerate(icen):
2014	C = np.array(cen)
2015	for invers in range(inv):
2016	for io,[M,T] in enumerate(SGData['SGOps']):
2017	idup = ((io+1)+100ic)(1-2*invers)
2018	XT = np.inner(M,X)+T
2019	if len(Uij):
2020	U = Uij2U(Uij)
2021	U = np.inner(M,np.inner(U,M).T)
2022	newUij = U2Uij(U)
2023	if invers:
2024	XT = -XT
2025	XT += C
2026	cell = np.zeros(3,dtype=np.int32)+cell0
2027	cellj = np.zeros(3,dtype=np.int32)
2028	if Move:
2029	newX,cellj = MoveToUnitCell(XT)
2030	else:
2031	newX = XT
2032	cell += cellj
2033	if All:
2034	if np.allclose(newX,X,atol=0.0002):
2035	idup = False
2036	else:
2037	if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
2038	idup = False
2039	if All or idup:
2040	XYZEquiv.append(newX)
2041	Idup.append(idup)
2042	Cell.append(cell)
2043	if len(Uij):
2044	UijEquiv.append(newUij)
2045	if len(SpnFlp):
2046	spnflp.append(SpnFlp[mj])
2047	else:
2048	spnflp.append(1)
2049	mj += 1
2050	if len(Uij):
2051	return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2052	else:
2053	return zip(XYZEquiv,Idup,Cell,spnflp)
2054
2055	def GenHKL(HKL,SGData):
2056	''' Generates all equivlent reflections including Friedel pairs
2057	:param HKL: [h,k,l] must be integral values
2058	:param SGData: space group data obtained from SpcGroup
2059	:returns: array Uniq: equivalent reflections
2060	'''
2061
2062	Ops = SGData['SGOps']
2063	OpM = np.array([op[0] for op in Ops])
2064	Uniq = np.inner(OpM,HKL)
2065	Uniq = list(Uniq)+list(-1*Uniq)
2066	return np.array(Uniq)
2067
2068	def GenHKLf(HKL,SGData):
2069	'''
2070	Uses old GSAS Fortran routine genhkl.for
2071
2072	:param HKL: [h,k,l] must be integral values for genhkl.for to work
2073	:param SGData: space group data obtained from SpcGroup
2074	:returns: iabsnt,mulp,Uniq,phi
2075
2076	* iabsnt = True if reflection is forbidden by symmetry
2077	* mulp = reflection multiplicity including Friedel pairs
2078	* Uniq = numpy array of equivalent hkl in descending order of h,k,l
2079	* phi = phase offset for each equivalent h,k,l
2080
2081	'''
2082	hklf = list(HKL)+[0,] #could be numpy array!
2083	Ops = SGData['SGOps']
2084	OpM = np.array([op[0] for op in Ops],order='F')
2085	OpT = np.array([op[1] for op in Ops])
2086	Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2087
2088	import pyspg
2089	Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2090	h,k,l,f = Uniq
2091	Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2092	phi = f[:Nuniq]
2093	return iabsnt,mulp,Uniq,phi
2094
2095	def checkSSLaue(HKL,SGData,SSGData):
2096	#Laue check here - Toss HKL if outside unique Laue part
2097	h,k,l,m = HKL
2098	if SGData['SGLaue'] == '2/m':
2099	if SGData['SGUniq'] == 'a':
2100	if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2101	return False
2102	elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2103	return False
2104	else:
2105	return True
2106	elif SGData['SGUniq'] == 'b':
2107	if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2108	return False
2109	elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2110	return False
2111	else:
2112	return True
2113	elif SGData['SGUniq'] == 'c':
2114	if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2115	return False
2116	elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2117	return False
2118	else:
2119	return True
2120	elif SGData['SGLaue'] == 'mmm':
2121	if 'a' in SSGData['modSymb']:
2122	if h == 0 and m < 0:
2123	return False
2124	else:
2125	return True
2126	elif 'b' in SSGData['modSymb']:
2127	if k == 0 and m < 0:
2128	return False
2129	else:
2130	return True
2131	elif 'g' in SSGData['modSymb']:
2132	if l == 0 and m < 0:
2133	return False
2134	else:
2135	return True
2136	else: #tetragonal, trigonal, hexagonal (& triclinic?)
2137	if l == 0 and m < 0:
2138	return False
2139	else:
2140	return True
2141
2142	def checkHKLextc(HKL,SGData):
2143	'''
2144	Checks if reflection extinct - does not check centering
2145
2146	:param HKL: [h,k,l]
2147	:param SGData: space group data obtained from SpcGroup
2148	:returns: True if extinct; False if allowed
2149
2150	'''
2151	Ops = SGData['SGOps']
2152	OpM = np.array([op[0] for op in Ops])
2153	OpT = np.array([op[1] for op in Ops])
2154	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2155	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2156	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2157	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
2158	if dhkl.any():
2159	continue
2160	else:
2161	if phkl%1.:
2162	return True
2163	return False
2164
2165	def checkMagextc(HKL,SGData):
2166	'''
2167	Checks if reflection magnetically extinct; does fullcheck (centering, too)
2168	uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2169
2170	:param HKL: [h,k,l]
2171	:param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2172	:returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2173
2174	'''
2175	Ops = SGData['SGOps']
2176	Ncen = len(SGData['SGCen'])
2177	OpM = np.array([op[0] for op in Ops])
2178	OpT = np.array([op[1] for op in Ops])
2179	if SGData['SGInv'] and not SGData['SGFixed']:
2180	OpM = np.vstack((OpM,-OpM))
2181	OpT = np.vstack((OpT,-OpT))%1.
2182	OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2183	OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2184	Spn = SGData['SpnFlp'][:len(OpM)]
2185	Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2186	DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2187	PHKL = np.reshape(np.cos(twopinp.inner(HKL,OpT))Mag,(-1,))[:,nxs,nxs]*OpM #compute T(R,theta) eq(7)
2188	Ftest = np.random.rand(3) #random magnetic moment
2189	Psum = np.zeros(3)
2190	nsum = 0.
2191	nA = 0
2192	for dhkl,phkl in zip(DHKL,PHKL):
2193	if not np.allclose(dhkl,HKL): #test for eq(5)
2194	continue
2195	else:
2196	nA += 1
2197	nsum += np.trace(phkl) #eq(8)
2198	pterm = np.inner(Ftest,phkl) #eq(9)
2199	Psum += pterm
2200	if nsum/nA > 1.: #only need to look at nA=1 frok eq(8)
2201	return False
2202	if np.allclose(Psum,np.zeros(3)):
2203	return True
2204	else:
2205	if np.inner(HKL,Psum):
2206	return True
2207	return False
2208
2209	def checkSSextc(HKL,SSGData):
2210	Ops = SSGData['SSGOps']
2211	OpM = np.array([op[0] for op in Ops])
2212	OpT = np.array([op[1] for op in Ops])
2213	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2214	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2215	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2216	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
2217	if dhkl.any():
2218	continue
2219	else:
2220	if phkl%1.:
2221	return False
2222	return True
2223
2224	################################################################################
2225	#### Site symmetry tables
2226	################################################################################
2227
2228	OprName = {
2229	'-6643': ['-1',1],'6479' : ['2(z)',2],'-6479': ['m(z)',3],
2230	'6481' : ['m(y)',4],'-6481': ['2(y)',5],'6641' : ['m(x)',6],
2231	'-6641': ['2(x)',7],'6591' : ['m(+-0)',8],'-6591': ['2(+-0)',9],
2232	'6531' : ['m(110)',10],'-6531': ['2(110)',11],'6537' : ['4(z)',12],
2233	'-6537': ['-4(z)',13],'975' : ['3(111)',14],'6456' : ['3',15],
2234	'-489' : ['3(+--)',16],'483' : ['3(-+-)',17],'-969' : ['3(--+)',18],
2235	'819' : ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' : ['m(0+-)',21],
2236	'-2431': ['2(0+-)',22],'-657' : ['m(xz)',23],'657' : ['2(xz)',24],
2237	'1943' : ['-4(x)',25],'-1943': ['4(x)',26],'-2429': ['m(yz)',27],
2238	'2429' : ['2(yz)',28],'639' : ['-4(y)',29],'-639' : ['4(y)',30],
2239	'-6484': ['2(010)',4],'6484' : ['m(010)',5],'-6668': ['2(100)',6],
2240	'6668' : ['m(100)',7],'-6454': ['2(120)',18],'6454' : ['m(120)',19],
2241	'-6638': ['2(210)',20],'6638' : ['m(210)',21], #search in SytSym ends at m(210)
2242	'2223' : ['3(+++)2',39],
2243	'6538' : ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2244	'2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489' :['-3(+--)1',45],
2245	'801' : ['4(y)1',46],'1945' : ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2246	'6585' : ['4(z)3',49],'6584' : ['3(z)2',50],'6666' : ['6(z)5 ',51],
2247	'6643' : ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2248	'-6666': ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2249	'-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2250	'969' :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2251	'-2151':['-3(+--)2',64], }
2252
2253	KNsym = {
2254	'0' :' 1 ','1' :' -1 ','64' :' 2(x)','32' :' m(x)',
2255	'97' :' 2/m(x)','16' :' 2(y)','8' :' m(y)','25' :' 2/m(y)',
2256	'2' :' 2(z)','4' :' m(z)','7' :' 2/m(z)','134217728' :' 2(yz)',
2257	'67108864' :' m(yz)','201326593' :' 2/m(yz)','2097152' :' 2(0+-)','1048576' :' m(0+-)',
2258	'3145729' :'2/m(0+-)','8388608' :' 2(xz)','4194304' :' m(xz)','12582913' :' 2/m(xz)',
2259	'524288' :' 2(+0-)','262144' :' m(+0-)','796433' :'2/m(+0-)','1024' :' 2(xy)',
2260	'512' :' m(xy)','1537' :' 2/m(xy)','256' :' 2(+-0)','128' :' m(+-0)',
2261	'385' :'2/m(+-0)','76' :' mm2(x)','52' :' mm2(y)','42' :' mm2(z)',
2262	'135266336' :' mm2(yz)','69206048' :'mm2(0+-)','8650760' :' mm2(xz)','4718600' :'mm2(+0-)',
2263	'1156' :' mm2(xy)','772' :'mm2(+-0)','82' :' 222 ','136314944' :' 222(x)',
2264	'8912912' :' 222(y)','1282' :' 222(z)','127' :' mmm ','204472417' :' mmm(x)',
2265	'13369369' :' mmm(y)','1927' :' mmm(z)','33554496' :' 4(x)','16777280' :' -4(x)',
2266	'50331745' :'4/m(x)' ,'169869394' :'422(x)','84934738' :'-42m(x)','101711948' :'4mm(x)',
2267	'254804095' :'4/mmm(x)','536870928 ':' 4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2268	'545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2269	'2050' :' 4(z)','4098' :' -4(z)','6151' :'4/m(z)','3410' :'422(z)',
2270	'4818' :'-42m(z)','2730' :'4mm(z)','8191' :'4/mmm(z)','8192' :' 3(111)',
2271	'8193' :' -3(111)','2629888' :' 32(111)','1319040' :' 3m(111)','3940737' :'-3m(111)',
2272	'32768' :' 3(+--)','32769' :' -3(+--)','10519552' :' 32(+--)','5276160' :' 3m(+--)',
2273	'15762945' :'-3m(+--)','65536' :' 3(-+-)','65537' :' -3(-+-)','134808576' :' 32(-+-)',
2274	'67437056' :' 3m(-+-)','202180097' :'-3m(-+-)','131072' :' 3(--+)','131073' :' -3(--+)',
2275	'142737664' :' 32(--+)','71434368' :' 3m(--+)','214040961' :'-3m(--+)','237650' :' 23 ',
2276	'237695' :' m3 ','715894098' :' 432 ','358068946' :' -43m ','1073725439':' m3m ',
2277	'68157504' :' mm2(d100)','4456464' :' mm2(d010)','642' :' mm2(d001)','153092172' :'-4m2(x)',
2278	'277348404' :'-4m2(y)','5418' :'-4m2(z)','1075726335':' 6/mmm ','1074414420':'-6m2(100)',
2279	'1075070124':'-6m2(120)','1075069650':' 6mm ','1074414890':' 622 ','1073758215':' 6/m ',
2280	'1073758212':' -6 ','1073758210':' 6 ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2281	'1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2282	'1073758209':' -3 ','1073758208':' 3 ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2283	'1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2284	'1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2285	'1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2286	'1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2287	'1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2288	'1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2289	'1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2290	'1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2291	'1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':' -1 ',
2292	}
2293
2294	NXUPQsym = {
2295	'1' :(28,29,28,28),'-1' :( 1,29,28, 0),'2(x)' :(12,18,12,25),'m(x)' :(25,18,12,25),
2296	'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),
2297	'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),
2298	'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),
2299	'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),
2300	'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),
2301	'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),
2302	'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),
2303	'mm2(yz)' :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)' :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2304	'mm2(xy)' :( 6,11, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222' :( 1,10, 0,-1),'222(x)' :( 1,13, 0,-1),
2305	'222(y)' :( 1,12, 0,-1),'222(z)' :( 1,11, 0,-1),'mmm' :( 1,10, 0,-1),'mmm(x)' :( 1,13, 0,-1),
2306	'mmm(y)' :( 1,12, 0,-1),'mmm(z)' :( 1,11, 0,-1),'4(x)' :(12, 4,12, 0),'-4(x)' :( 1, 4,12, 0),
2307	'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),
2308	'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),
2309	'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,),
2310	'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),
2311	'-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),
2312	'-3(111)' :( 1, 5, 2, 0),'32(111)' :( 1, 5, 0, 2),'3m(111)' :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2313	'3(+--)' :( 5, 8, 5, 0),'-3(+--)' :( 1, 8, 5, 0),'32(+--)' :( 1, 8, 0, 5),'3m(+--)' :( 5, 8, 0, 5),
2314	'-3m(+--)' :( 1, 8, 0,-1),'3(-+-)' :( 4, 7, 4, 0),'-3(-+-)' :( 1, 7, 4, 0),'32(-+-)' :( 1, 7, 0, 4),
2315	'3m(-+-)' :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)' :( 3, 6, 3, 0),'-3(--+)' :( 1, 6, 3, 0),
2316	'32(--+)' :( 1, 6, 0, 3),'3m(--+)' :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23' :( 1, 1, 0, 0),
2317	'm3' :( 1, 1, 0, 0),'432' :( 1, 1, 0, 0),'-43m' :( 1, 1, 0, 0),'m3m' :( 1, 1, 0, 0),
2318	'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)' :( 1, 4, 0,-1),
2319	'-4m2(y)' :( 1, 3, 0,-1),'-4m2(z)' :( 1, 2, 0,-1),'6/mmm' :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2320	'-6m2(120)':( 1, 9, 0,-1),'6mm' :(14, 9, 0,-1),'622' :( 1, 9, 0,-1),'6/m' :( 1, 9,14,-1),
2321	'-6' :( 1, 9,14, 0),'6' :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2322	'3m(100)' :(14, 9, 0,14),'3m(120)' :(14, 9, 0,14),'32(100)' :( 1, 9, 0,14),'32(120)' :( 1, 9, 0,14),
2323	'-3' :( 1, 9,14, 0),'3' :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2324	'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2325	'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2326	'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2327	'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),
2328	'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),
2329	'm(001)' :(23,16,14,23),'m(100)' :(26,25,12,26),'m(010)' :(27,28,13,27),'m(110)' :(18,19, 6,18),
2330	'm(120)' :(24,27,15,24),'m(210)' :(25,26,16,25),'m(+-0)' :(17,20, 7,17),'2(001)' :(14,16,14,23),
2331	'2(100)' :(12,25,12,26),'2(010)' :(13,28,13,27),'2(110)' :( 6,19, 6,18),'2(120)' :(15,27,15,24),
2332	'2(210)' :(16,26,16,25),'2(+-0)' :( 7,20, 7,17),'-1' :( 1,29,28, 0)
2333	}
2334
2335	CSxinel = [[], # 0th empty - indices are Fortran style
2336	[[0,0,0],[ 0.0, 0.0, 0.0]], #1 0 0 0
2337	[[1,1,1],[ 1.0, 1.0, 1.0]], #2 X X X
2338	[[1,1,1],[ 1.0, 1.0,-1.0]], #3 X X -X
2339	[[1,1,1],[ 1.0,-1.0, 1.0]], #4 X -X X
2340	[[1,1,1],[ 1.0,-1.0,-1.0]], #5 -X X X
2341	[[1,1,0],[ 1.0, 1.0, 0.0]], #6 X X 0
2342	[[1,1,0],[ 1.0,-1.0, 0.0]], #7 X -X 0
2343	[[1,0,1],[ 1.0, 0.0, 1.0]], #8 X 0 X
2344	[[1,0,1],[ 1.0, 0.0,-1.0]], #9 X 0 -X
2345	[[0,1,1],[ 0.0, 1.0, 1.0]], #10 0 Y Y
2346	[[0,1,1],[ 0.0, 1.0,-1.0]], #11 0 Y -Y
2347	[[1,0,0],[ 1.0, 0.0, 0.0]], #12 X 0 0
2348	[[0,1,0],[ 0.0, 1.0, 0.0]], #13 0 Y 0
2349	[[0,0,1],[ 0.0, 0.0, 1.0]], #14 0 0 Z
2350	[[1,1,0],[ 1.0, 2.0, 0.0]], #15 X 2X 0
2351	[[1,1,0],[ 2.0, 1.0, 0.0]], #16 2X X 0
2352	[[1,1,2],[ 1.0, 1.0, 1.0]], #17 X X Z
2353	[[1,1,2],[ 1.0,-1.0, 1.0]], #18 X -X Z
2354	[[1,2,1],[ 1.0, 1.0, 1.0]], #19 X Y X
2355	[[1,2,1],[ 1.0, 1.0,-1.0]], #20 X Y -X
2356	[[1,2,2],[ 1.0, 1.0, 1.0]], #21 X Y Y
2357	[[1,2,2],[ 1.0, 1.0,-1.0]], #22 X Y -Y
2358	[[1,2,0],[ 1.0, 1.0, 0.0]], #23 X Y 0
2359	[[1,0,2],[ 1.0, 0.0, 1.0]], #24 X 0 Z
2360	[[0,1,2],[ 0.0, 1.0, 1.0]], #25 0 Y Z
2361	[[1,1,2],[ 1.0, 2.0, 1.0]], #26 X 2X Z
2362	[[1,1,2],[ 2.0, 1.0, 1.0]], #27 2X X Z
2363	[[1,2,3],[ 1.0, 1.0, 1.0]], #28 X Y Z
2364	]
2365
2366	CSuinel = [[], # 0th empty - indices are Fortran style
2367	[[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
2368	[[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
2369	[[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
2370	[[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
2371	[[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
2372	[[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
2373	[[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
2374	[[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
2375	[[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
2376	[[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
2377	[[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
2378	[[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
2379	[[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
2380	[[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
2381	[[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
2382	[[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
2383	[[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
2384	[[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
2385	[[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
2386	[[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
2387	[[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
2388	[[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
2389	[[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
2390	[[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
2391	[[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
2392	[[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
2393	[[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
2394	[[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
2395	[[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
2396	]
2397
2398	################################################################################
2399	#### Site symmetry routines
2400	################################################################################
2401
2402	def GetOprPtrName(key):
2403	'Needs a doc string'
2404	try:
2405	oprName = OprName[key][0]
2406	except KeyError:
2407	return key
2408	return oprName.replace('(','').replace(')','')
2409
2410	def GetOprPtrNumber(key):
2411	'Needs a doc string'
2412	try:
2413	return OprName[key][1]
2414	except KeyError:
2415	return key
2416
2417	def GetOprName(key):
2418	'Needs a doc string'
2419	return OprName[key][0]
2420
2421	def GetKNsym(key):
2422	'Needs a doc string'
2423	try:
2424	return KNsym[key].strip()
2425	except KeyError:
2426	return 'sp'
2427
2428	def GetNXUPQsym(siteSym):
2429	'''
2430	The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2431	'''
2432	return NXUPQsym[siteSym]
2433
2434	def GetCSxinel(siteSym):
2435	"returns Xyz terms, multipliers, GUI flags"
2436	indx = GetNXUPQsym(siteSym.strip())
2437	return CSxinel[indx[0]]
2438
2439	def GetCSuinel(siteSym):
2440	"returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2441	indx = GetNXUPQsym(siteSym.strip())
2442	return CSuinel[indx[1]]
2443
2444	def GetCSpqinel(SpnFlp,dupDir):
2445	"returns Mxyz terms, multipliers, GUI flags"
2446	CSI = [[1,2,3],[1.0,1.0,1.0]]
2447	for sopr in dupDir:
2448	# print (sopr,dupDir[sopr])
2449	opr = sopr.replace("'",'')
2450	indx = GetNXUPQsym(opr)
2451	if SpnFlp[dupDir[sopr]] > 0:
2452	csi = CSxinel[indx[2]] #P
2453	else:
2454	csi = CSxinel[indx[3]] #Q
2455	# print(opr,indx,csi,CSI)
2456	if not len(csi):
2457	return [[0,0,0],[0.,0.,0.]]
2458	for kcs in [0,1,2]:
2459	if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2460	jcs = CSI[0][kcs]
2461	for ics in [0,1,2]:
2462	if CSI[0][ics] == jcs:
2463	CSI[0][ics] = 0
2464	CSI[1][ics] = 0.
2465	elif CSI[0][ics] > jcs:
2466	CSI[0][ics] = CSI[0][ics]-1
2467	elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2468	CSI[1][kcs] = csi[1][kcs]
2469	elif CSI[0][kcs] >= csi[0][kcs]:
2470	CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2471	if CSI[1][kcs] != csi[1][kcs]:
2472	if CSI[1][kcs] == 1.:
2473	CSI[1][kcs] = csi[1][kcs]
2474	# print(CSI)
2475	return CSI
2476
2477	def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2478	phase = np.sum(XYZ*wave)
2479	ssopinv = nl.inv(ssop[0])
2480	mst = ssopinv[3][:3]
2481	epsinv = ssopinv[3][3]
2482	sdet = nl.det(sop[0])
2483	ssdet = nl.det(ssop[0])
2484	dtau = mst(XYZ-sop[1])-epsinvssop[1][3]
2485	dT = 1.0
2486	if np.any(dtau%.5):
2487	sumdtau = np.sum(dtau%.5)
2488	dT = 0.
2489	if np.abs(sumdtau-.5) > 1.e-4:
2490	dT = np.tan(np.pi*sumdtau)
2491	tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2492	return sdet,ssdet,dtau,dT,tauT
2493
2494	def OpsfromStringOps(A,SGData,SSGData):
2495	SGOps = SGData['SGOps']
2496	SSGOps = SSGData['SSGOps']
2497	Ax = A.split('+')
2498	Ax[0] = int(Ax[0])
2499	iC = 1
2500	if Ax[0] < 0:
2501	iC = -1
2502	Ax[0] = abs(Ax[0])
2503	nA = Ax[0]%100-1
2504	nC = Ax[0]//100
2505	unit = [0,0,0]
2506	if len(Ax) > 1:
2507	unit = eval('['+Ax[1]+']')
2508	return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2509
2510	def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2511
2512	def orderParms(CSI):
2513	parms = [0,]
2514	for csi in CSI:
2515	for i in [0,1,2]:
2516	if csi[i] not in parms:
2517	parms.append(csi[i])
2518	for csi in CSI:
2519	for i in [0,1,2]:
2520	csi[i] = parms.index(csi[i])
2521	return CSI
2522
2523	def fracCrenel(tau,Toff,Twid):
2524	Tau = (tau-Toff[:,nxs])%1.
2525	A = np.where(Tau<Twid[:,nxs],1.,0.)
2526	return A
2527
2528	def fracFourier(tau,nH,fsin,fcos):
2529	SA = np.sin(2.nHnp.pi*tau)
2530	CB = np.cos(2.nHnp.pi*tau)
2531	A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2532	B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2533	return A+B
2534
2535	def posFourier(tau,nH,psin,pcos):
2536	SA = np.sin(2nHnp.pi*tau)
2537	CB = np.cos(2nHnp.pi*tau)
2538	A = SA[nxs,nxs,:]*psin[:,:,nxs]
2539	B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2540	return A+B
2541
2542	def posSawtooth(tau,Toff,slopes):
2543	Tau = (tau-Toff)%1.
2544	A = slopes[:,nxs]*Tau
2545	return A
2546
2547	def posZigZag(tau,Tmm,XYZmax):
2548	DT = Tmm[1]-Tmm[0]
2549	slopeUp = 2.*XYZmax/DT
2550	slopeDn = 2.*XYZmax/(1.-DT)
2551	A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp((t-Tmm[0])%1.),XYZmax-slopeDn((t-Tmm[1])%1.)) for t in tau])
2552	return A
2553
2554	def posBlock(tau,Tmm,XYZmax):
2555	A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2556	return A
2557
2558	def DoFrac():
2559	delt2 = np.eye(2)*0.001
2560	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2561	dFTP = []
2562	if siteSym == '1':
2563	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2564	elif siteSym == '-1':
2565	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2566	else:
2567	FSC = np.ones(2,dtype='i')
2568	CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2569	if 'Crenel' in waveType:
2570	dF = np.zeros_like(tau)
2571	else:
2572	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2573	dFT = np.zeros_like(dF)
2574	dFTP = []
2575	for i in SdIndx:
2576	sop = Sop[i]
2577	ssop = SSop[i]
2578	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2579	fsc = np.ones(2,dtype='i')
2580	if 'Crenel' in waveType:
2581	dFT = np.zeros_like(tau)
2582	fsc = [1,1]
2583	else: #Fourier
2584	dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2585	dFT = nl.det(sop[0])*dFT
2586	dFT = dFT[:,np.argsort(tauT)]
2587	dFT[0] *= ssdet
2588	dFT[1] *= sdet
2589	dFTP.append(dFT)
2590
2591	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2592	fsc = [1,1]
2593	if dT:
2594	CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2595	else:
2596	CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2597	FSC = np.zeros(2,dtype='i')
2598	return CSI,dF,dFTP
2599	else:
2600	for i in range(2):
2601	if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2602	fsc[i] = 1
2603	else:
2604	fsc[i] = 0
2605	FSC &= fsc
2606	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2607	n = -1
2608	for i,F in enumerate(FSC):
2609	if F:
2610	n += 1
2611	CSI[0][i] = n+1
2612	CSI[1][i] = 1.0
2613
2614	return CSI,dF,dFTP
2615
2616	def DoXYZ():
2617	delt4 = np.ones(4)*0.001
2618	delt5 = np.ones(5)*0.001
2619	delt6 = np.eye(6)*0.001
2620	if 'Fourier' in waveType:
2621	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2622	#3x6x12 modulated position array (X,Spos,tau)& force positive
2623	elif waveType == 'Sawtooth':
2624	dX = posSawtooth(tau,delt4[0],delt4[1:])
2625	elif waveType in ['ZigZag','Block']:
2626	if waveType == 'ZigZag':
2627	dX = posZigZag(tau,delt5[:2],delt5[2:])
2628	else:
2629	dX = posBlock(tau,delt5[:2],delt5[2:])
2630	dXTP = []
2631	if siteSym == '1':
2632	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2633	elif siteSym == '-1':
2634	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.],]
2635	else:
2636	if 'Fourier' in waveType:
2637	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2638	elif waveType == 'Sawtooth':
2639	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2640	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2641	elif waveType in ['ZigZag','Block']:
2642	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2643	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2644	XSC = np.ones(6,dtype='i')
2645	dXTP = []
2646	for i in SdIndx:
2647	sop = Sop[i]
2648	ssop = SSop[i]
2649	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2650	xsc = np.ones(6,dtype='i')
2651	if 'Fourier' in waveType:
2652	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2653	elif waveType == 'Sawtooth':
2654	dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,nxs,nxs]
2655	elif waveType == 'ZigZag':
2656	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2657	elif waveType == 'Block':
2658	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2659	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2660	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2661	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2662	dXTP.append(dXT)
2663	if waveType == 'Fourier':
2664	for i in range(3):
2665	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2666	xsc[i] = 0
2667	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2668	xsc[i+3] = 0
2669	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2670	xsc[3:6] = 0
2671	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2672	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2673	if dT:
2674	if '(x)' in siteSym:
2675	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2676	if 'm' in siteSym and len(SdIndx) == 1:
2677	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2678	elif '(y)' in siteSym:
2679	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2680	if 'm' in siteSym and len(SdIndx) == 1:
2681	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2682	elif '(z)' in siteSym:
2683	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2684	if 'm' in siteSym and len(SdIndx) == 1:
2685	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2686	else:
2687	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2688	if '4/mmm' in laue:
2689	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2690	if '(xy)' in siteSym:
2691	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2692	if dT:
2693	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2694	else:
2695	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2696	if '(xy)' in siteSym or '(+-0)' in siteSym:
2697	mul = 1
2698	if '(+-0)' in siteSym:
2699	mul = -1
2700	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2701	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2702	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2703	xsc[3:5] = 0
2704	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2705	CSI[0][:2] = [[12,0,0],[12,0,0]]
2706	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2707	xsc[:2] = 0
2708	XSC &= xsc
2709	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2710	if waveType == 'Fourier':
2711	n = -1
2712	if debug: print (XSC)
2713	for i,X in enumerate(XSC):
2714	if X:
2715	n += 1
2716	CSI[0][i][0] = n+1
2717	CSI[1][i][0] = 1.0
2718
2719	return list(CSI),dX,dXTP
2720
2721	def DoUij():
2722	delt12 = np.eye(12)*0.0001
2723	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2724	dUTP = []
2725	if siteSym == '1':
2726	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2727	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2728	elif siteSym == '-1':
2729	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2730	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2731	else:
2732	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2733	USC = np.ones(12,dtype='i')
2734	dUTP = []
2735	dtau = 0.
2736	for i in SdIndx:
2737	sop = Sop[i]
2738	ssop = SSop[i]
2739	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2740	usc = np.ones(12,dtype='i')
2741	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2742	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2743	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2744	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2745	dUT = dUT[:,:,np.argsort(tauT)]
2746	dUT[:,:6,:] =(ssdetsdet)
2747	dUTP.append(dUT)
2748	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2749	if dT:
2750	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2751	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2752	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2753	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2754	else:
2755	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2756	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2757	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2758	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2759	if 'mm2(x)' in siteSym and dT:
2760	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2761	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2762	elif '(xy)' in siteSym and dT:
2763	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2764	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2765	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2766	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2767	elif '(x)' in siteSym and dT:
2768	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2769	elif '(y)' in siteSym and dT:
2770	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2771	elif '(z)' in siteSym and dT:
2772	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2773	for i in range(6):
2774	if not USC[i]:
2775	CSI[0][i] = [0,0,0]
2776	CSI[1][i] = [0.,0.,0.]
2777	CSI[0][i+6] = [0,0,0]
2778	CSI[1][i+6] = [0.,0.,0.]
2779	else:
2780	for i in range(6):
2781	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2782	usc[i] = 0
2783	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2784	usc[i+6] = 0
2785	if np.any(dUT[1,0,:]):
2786	if '4/m' in siteSym:
2787	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2788	if ssop[1][3]:
2789	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2790	usc[9] = 1
2791	else:
2792	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2793	usc[9] = 0
2794	elif '4' in siteSym:
2795	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2796	CSI[0][:2] = [[11,0,0],[11,0,0]]
2797	if ssop[1][3]:
2798	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2799	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2800	usc[2] = 0
2801	usc[8] = 0
2802	usc[3] = 1
2803	usc[9] = 1
2804	else:
2805	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2806	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2807	usc[2] = 1
2808	usc[8] = 1
2809	usc[3] = 0
2810	usc[9] = 0
2811	elif 'xy' in siteSym or '+-0' in siteSym:
2812	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2813	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2814	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2815	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2816	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2817	usc[4:6] = 0
2818	usc[6:8] = 0
2819
2820	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2821	USC &= usc
2822	if debug: print (USC)
2823	if not np.any(dtau%.5):
2824	n = -1
2825	for i,U in enumerate(USC):
2826	if U:
2827	n += 1
2828	CSI[0][i][0] = n+1
2829	CSI[1][i][0] = 1.0
2830
2831	return list(CSI),dU,dUTP
2832
2833	def DoMag():
2834	delt6 = np.eye(6)*0.001
2835	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2836	dMTP = []
2837	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2838	if siteSym == '1':
2839	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2840	elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2841	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2842	else:
2843	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2844	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2845	MSC = np.ones(6,dtype='i')
2846	dMTP = []
2847	for i in SdIndx:
2848	sop = Sop[i]
2849	ssop = SSop[i]
2850	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2851	msc = np.ones(6,dtype='i')
2852	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2853	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2854	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2855	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2856	dMTP.append(dMT)
2857	for i in range(3):
2858	if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2859	msc[i] = 0
2860	if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2861	msc[i+3] = 0
2862	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2863	msc[3:6] = 0
2864	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2865	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2866	if dT:
2867	if '(x)' in siteSym:
2868	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2869	if 'm' in siteSym and len(SdIndx) == 1:
2870	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2871	elif '(y)' in siteSym:
2872	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2873	if 'm' in siteSym and len(SdIndx) == 1:
2874	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2875	elif '(z)' in siteSym:
2876	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2877	if 'm' in siteSym and len(SdIndx) == 1:
2878	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2879	else:
2880	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2881	if '4/mmm' in laue:
2882	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2883	if '(xy)' in siteSym:
2884	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2885	if dT:
2886	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2887	else:
2888	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2889	if '(xy)' in siteSym or '(+-0)' in siteSym:
2890	mul = 1
2891	if '(+-0)' in siteSym:
2892	mul = -1
2893	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2894	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2895	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2896	msc[3:5] = 0
2897	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2898	CSI[0][:2] = [[12,0,0],[12,0,0]]
2899	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2900	msc[:2] = 0
2901	MSC &= msc
2902	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2903	n = -1
2904	if debug: print (MSC)
2905	for i,M in enumerate(MSC):
2906	if M:
2907	n += 1
2908	CSI[0][i][0] = n+1
2909	CSI[1][i][0] = 1.0
2910
2911	return list(CSI),dM,dMTP
2912
2913	if debug: print ('super space group: '+SSGData['SSpGrp'])
2914	xyz = np.array(XYZ)%1.
2915	SGOps = copy.deepcopy(SGData['SGOps'])
2916	laue = SGData['SGLaue']
2917	siteSym = SytSym(XYZ,SGData)[0].strip()
2918	if debug: print ('siteSym: '+siteSym)
2919	SSGOps = copy.deepcopy(SSGData['SSGOps'])
2920	#expand ops to include inversions if any
2921	if SGData['SGInv']:
2922	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2923	SGOps.append([-op[0],-op[1]%1.])
2924	SSGOps.append([-sop[0],-sop[1]%1.])
2925	#build set of sym ops around special position
2926	SSop = []
2927	Sop = []
2928	Sdtau = []
2929	for iop,Op in enumerate(SGOps):
2930	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2931	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2932	SSop.append(SSGOps[iop])
2933	Sop.append(SGOps[iop])
2934	ssopinv = nl.inv(SSGOps[iop][0])
2935	mst = ssopinv[3][:3]
2936	epsinv = ssopinv[3][3]
2937	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
2938	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
2939	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2940	#setup displacement arrays
2941	tau = np.linspace(-1,1,49,True)
2942	#make modulation arrays - one parameter at a time
2943	if Stype == 'Sfrac':
2944	CSI,dF,dFTP = DoFrac()
2945	elif Stype == 'Spos':
2946	CSI,dF,dFTP = DoXYZ()
2947	CSI[0] = orderParms(CSI[0])
2948	elif Stype == 'Sadp':
2949	CSI,dF,dFTP = DoUij()
2950	CSI[0] = orderParms(CSI[0])
2951	elif Stype == 'Smag':
2952	CSI,dF,dFTP = DoMag()
2953
2954	if debug:
2955	return CSI,dF,dFTP
2956	else:
2957	return CSI,[],[]
2958
2959	def MustrainNames(SGData):
2960	'Needs a doc string'
2961	laue = SGData['SGLaue']
2962	uniq = SGData['SGUniq']
2963	if laue in ['m3','m3m']:
2964	return ['S400','S220']
2965	elif laue in ['6/m','6/mmm','3m1']:
2966	return ['S400','S004','S202']
2967	elif laue in ['31m','3']:
2968	return ['S400','S004','S202','S211']
2969	elif laue in ['3R','3mR']:
2970	return ['S400','S220','S310','S211']
2971	elif laue in ['4/m','4/mmm']:
2972	return ['S400','S004','S220','S022']
2973	elif laue in ['mmm']:
2974	return ['S400','S040','S004','S220','S202','S022']
2975	elif laue in ['2/m']:
2976	SHKL = ['S400','S040','S004','S220','S202','S022']
2977	if uniq == 'a':
2978	SHKL += ['S013','S031','S211']
2979	elif uniq == 'b':
2980	SHKL += ['S301','S103','S121']
2981	elif uniq == 'c':
2982	SHKL += ['S130','S310','S112']
2983	return SHKL
2984	else:
2985	SHKL = ['S400','S040','S004','S220','S202','S022']
2986	SHKL += ['S310','S103','S031','S130','S301','S013']
2987	SHKL += ['S211','S121','S112']
2988	return SHKL
2989
2990	def HStrainVals(HSvals,SGData):
2991	laue = SGData['SGLaue']
2992	uniq = SGData['SGUniq']
2993	DIJ = np.zeros(6)
2994	if laue in ['m3','m3m']:
2995	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2996	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2997	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2998	elif laue in ['3R','3mR']:
2999	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3000	elif laue in ['4/m','4/mmm']:
3001	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3002	elif laue in ['mmm']:
3003	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3004	elif laue in ['2/m']:
3005	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3006	if uniq == 'a':
3007	DIJ[5] = HSvals[3]
3008	elif uniq == 'b':
3009	DIJ[4] = HSvals[3]
3010	elif uniq == 'c':
3011	DIJ[3] = HSvals[3]
3012	else:
3013	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3014	return DIJ
3015
3016	def HStrainNames(SGData):
3017	'Needs a doc string'
3018	laue = SGData['SGLaue']
3019	uniq = SGData['SGUniq']
3020	if laue in ['m3','m3m']:
3021	return ['D11','eA'] #add cubic strain term
3022	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3023	return ['D11','D33']
3024	elif laue in ['3R','3mR']:
3025	return ['D11','D12']
3026	elif laue in ['4/m','4/mmm']:
3027	return ['D11','D33']
3028	elif laue in ['mmm']:
3029	return ['D11','D22','D33']
3030	elif laue in ['2/m']:
3031	Dij = ['D11','D22','D33']
3032	if uniq == 'a':
3033	Dij += ['D23']
3034	elif uniq == 'b':
3035	Dij += ['D13']
3036	elif uniq == 'c':
3037	Dij += ['D12']
3038	return Dij
3039	else:
3040	Dij = ['D11','D22','D33','D12','D13','D23']
3041	return Dij
3042
3043	def MustrainCoeff(HKL,SGData):
3044	'Needs a doc string'
3045	#NB: order of terms is the same as returned by MustrainNames
3046	laue = SGData['SGLaue']
3047	uniq = SGData['SGUniq']
3048	h,k,l = HKL
3049	Strm = []
3050	if laue in ['m3','m3m']:
3051	Strm.append(h4+k4+l**4)
3052	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3053	elif laue in ['6/m','6/mmm','3m1']:
3054	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3055	Strm.append(l**4)
3056	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3057	elif laue in ['31m','3']:
3058	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3059	Strm.append(l**4)
3060	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3061	Strm.append(4.0hkl(h+k))
3062	elif laue in ['3R','3mR']:
3063	Strm.append(h4+k4+l**4)
3064	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3065	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
3066	Strm.append(4.0(klh2+hlk2+hkl*2))
3067	elif laue in ['4/m','4/mmm']:
3068	Strm.append(h4+k4)
3069	Strm.append(l**4)
3070	Strm.append(3.0(hk)**2)
3071	Strm.append(3.0((hl)*2+(kl)**2))
3072	elif laue in ['mmm']:
3073	Strm.append(h**4)
3074	Strm.append(k**4)
3075	Strm.append(l**4)
3076	Strm.append(3.0(hk)**2)
3077	Strm.append(3.0(hl)**2)
3078	Strm.append(3.0(kl)**2)
3079	elif laue in ['2/m']:
3080	Strm.append(h**4)
3081	Strm.append(k**4)
3082	Strm.append(l**4)
3083	Strm.append(3.0(hk)**2)
3084	Strm.append(3.0(hl)**2)
3085	Strm.append(3.0(kl)**2)
3086	if uniq == 'a':
3087	Strm.append(2.0kl**3)
3088	Strm.append(2.0lk**3)
3089	Strm.append(4.0klh*2)
3090	elif uniq == 'b':
3091	Strm.append(2.0lh**3)
3092	Strm.append(2.0hl**3)
3093	Strm.append(4.0hlk*2)
3094	elif uniq == 'c':
3095	Strm.append(2.0hk**3)
3096	Strm.append(2.0kh**3)
3097	Strm.append(4.0hkl*2)
3098	else:
3099	Strm.append(h**4)
3100	Strm.append(k**4)
3101	Strm.append(l**4)
3102	Strm.append(3.0(hk)**2)
3103	Strm.append(3.0(hl)**2)
3104	Strm.append(3.0(kl)**2)
3105	Strm.append(2.0kh**3)
3106	Strm.append(2.0hl**3)
3107	Strm.append(2.0lk**3)
3108	Strm.append(2.0hk**3)
3109	Strm.append(2.0lh**3)
3110	Strm.append(2.0kl**3)
3111	Strm.append(4.0klh*2)
3112	Strm.append(4.0hlk*2)
3113	Strm.append(4.0khl*2)
3114	return Strm
3115
3116	def MuShklMean(SGData,Amat,Shkl):
3117
3118	def genMustrain(xyz,Shkl):
3119	uvw = np.inner(Amat.T,xyz)
3120	Strm = np.array(MustrainCoeff(uvw,SGData))
3121	Sum = np.sum(np.multiply(Shkl,Strm))
3122	Sum = np.where(Sum > 0.01,Sum,0.01)
3123	Sum = np.sqrt(Sum)
3124	return Sum*xyz
3125
3126	PHI = np.linspace(0.,360.,30,True)
3127	PSI = np.linspace(0.,180.,30,True)
3128	X = np.outer(npcosd(PHI),npsind(PSI))
3129	Y = np.outer(npsind(PHI),npsind(PSI))
3130	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3131	XYZ = np.dstack((X,Y,Z))
3132	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3133	return np.sqrt(np.sum(XYZ**2)/900.)
3134
3135	def Muiso2Shkl(muiso,SGData,cell):
3136	"this is to convert isotropic mustrain to generalized Shkls"
3137	import GSASIIlattice as G2lat
3138	A = G2lat.cell2AB(cell)[0]
3139
3140	def minMus(Shkl,muiso,H,SGData,A):
3141	U = np.inner(A.T,H)
3142	S = np.array(MustrainCoeff(U,SGData))
3143	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,nxs]),axis=0))
3144	rad = np.sqrt(np.sum((Sum[:,nxs]H)*2,axis=1))
3145	return (muiso-rad)**2
3146
3147	laue = SGData['SGLaue']
3148	PHI = np.linspace(0.,360.,60,True)
3149	PSI = np.linspace(0.,180.,60,True)
3150	X = np.outer(npsind(PHI),npsind(PSI))
3151	Y = np.outer(npcosd(PHI),npsind(PSI))
3152	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3153	HKL = np.dstack((X,Y,Z))
3154	if laue in ['m3','m3m']:
3155	S0 = [1000.,1000.]
3156	elif laue in ['6/m','6/mmm','3m1']:
3157	S0 = [1000.,1000.,1000.]
3158	elif laue in ['31m','3']:
3159	S0 = [1000.,1000.,1000.,1000.]
3160	elif laue in ['3R','3mR']:
3161	S0 = [1000.,1000.,1000.,1000.]
3162	elif laue in ['4/m','4/mmm']:
3163	S0 = [1000.,1000.,1000.,1000.]
3164	elif laue in ['mmm']:
3165	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3166	elif laue in ['2/m']:
3167	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3168	else:
3169	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
3170	1000.,1000.,0.,0.,0.]
3171	S0 = np.array(S0)
3172	HKL = np.reshape(HKL,(-1,3))
3173	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3174	return result[0]
3175
3176	def PackRot(SGOps):
3177	IRT = []
3178	for ops in SGOps:
3179	M = ops[0]
3180	irt = 0
3181	for j in range(2,-1,-1):
3182	for k in range(2,-1,-1):
3183	irt *= 3
3184	irt += M[k][j]
3185	IRT.append(int(irt))
3186	return IRT
3187
3188	def SytSym(XYZ,SGData):
3189	'''
3190	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3191
3192	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
3193	:param SGData: from SpcGroup
3194	:Returns: a four element tuple:
3195
3196	* The 1st element is a code for the site symmetry (see GetKNsym)
3197	* The 2nd element is the site multiplicity
3198	* Ndup number of overlapping operators
3199	* dupDir Dict - dictionary of overlapping operators
3200
3201	'''
3202	Mult = 1
3203	Isym = 0
3204	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3205	Isym = 1073741824
3206	Jdup = 0
3207	Ndup = 0
3208	dupDir = {}
3209	inv = SGData['SGInv']+1
3210	icen = SGData['SGCen']
3211	if SGData['SGFixed']: #already in list of operators
3212	inv = 1
3213	Xeqv = list(GenAtom(XYZ,SGData,True))
3214	# for xeqv in Xeqv: print(xeqv)
3215	IRT = PackRot(SGData['SGOps'])
3216	L = -1
3217	for ic,cen in enumerate(icen):
3218	for invers in range(int(inv)):
3219	for io,ops in enumerate(SGData['SGOps']):
3220	irtx = (1-2invers)IRT[io]
3221	L += 1
3222	if not Xeqv[L][1]:
3223	Ndup = io
3224	Jdup += 1
3225	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
3226	if jx < 39:
3227	px = GetOprName(str(irtx))
3228	if Xeqv[L][-1] < 0:
3229	if '(' in px:
3230	px = px.split('(')
3231	px[0] += "'"
3232	px = '('.join(px)
3233	else:
3234	px += "'"
3235	dupDir[px] = L
3236	Isym += 2**(jx-1)
3237	if Isym == 1073741824: Isym = 0
3238	Mult = len(SGData['SGOps'])len(SGData['SGCen'])inv//Jdup
3239
3240	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3241
3242	def MagSytSym(SytSym,dupDir,SGData):
3243	'''
3244	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3245	'''
3246	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3247	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3248	# print('dupDir',dupDir)
3249	SplitSytSym = SytSym.split('(')
3250	if SGData['SGGray']:
3251	return SytSym+"1'"
3252	if SytSym == '1': #genersl position
3253	return SytSym
3254	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
3255	try:
3256	MagSytSym = SGData['MagSpGrp'].split()[1]
3257	except IndexError:
3258	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3259	if len(SplitSytSym) > 1:
3260	MagSytSym += '('+SplitSytSym[1]
3261	return MagSytSym
3262	if len(dupDir) == 1:
3263	return list(dupDir.keys())[0]
3264
3265
3266	if '2/m' in SytSym: #done I think; last 2wo might be not needed
3267	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3268	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3269	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3270	'(110)':['2(110)','m(110)']}
3271
3272	elif '4/mmm' in SytSym:
3273	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
3274	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
3275	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
3276	elif '4mm' in SytSym:
3277	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3278	elif '422' in SytSym:
3279	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3280	elif '-4m2' in SytSym:
3281	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3282	elif '-42m' in SytSym:
3283	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3284	elif '-4' in SytSym:
3285	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3286	elif '4' in SytSym:
3287	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3288
3289	elif '222' in SytSym:
3290	ops = {'':['2(x)','2(y)','2(z)'],
3291	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3292	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3293	'(110)':['2(z)','2(110)','2(+-0)',],}
3294	elif 'mm2' in SytSym:
3295	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3296	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3297	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3298	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3299	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3300	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3301	'(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3302	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3303	'(110)':['m(z)','m(+-0)','2(110)',],}
3304	elif 'mmm' in SytSym:
3305	ops = {'':['m(x)','m(y)','m(z)'],
3306	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3307	'(110)':['m(z)','m(110)','m(+-0)',],
3308	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3309
3310	elif '32' in SytSym:
3311	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3312	elif '23' in SytSym:
3313	ops = {'':['2(x)','3(111)']}
3314	elif 'm3' in SytSym:
3315	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3316	elif '3m' in SytSym:
3317	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3318	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3319	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3320
3321	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]: #not simple cases
3322	MagSytSym = SytSym
3323	if "-1'" in dupDir:
3324	if '-6' in SytSym:
3325	MagSytSym = MagSytSym.replace('-6',"-6'")
3326	elif '-3m' in SytSym:
3327	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3328	elif '-3' in SytSym:
3329	MagSytSym = MagSytSym.replace('-3',"-3'")
3330	elif '-6m2' in SytSym:
3331	if "m'(110)" in dupDir:
3332	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3333	elif '6/m' in SytSym:
3334	if "m'(z)" in dupDir:
3335	MagSytSym = "6'/m'"
3336	elif '6mm' in SytSym:
3337	if "m'(110)" in dupDir:
3338	MagSytSym = "6'm'm"
3339	return MagSytSym
3340	try:
3341	axis = '('+SytSym.split('(')[1]
3342	except IndexError:
3343	axis = ''
3344	MagSytSym = ''
3345	for m in ops[axis]:
3346	if m in dupDir:
3347	MagSytSym += m.split('(')[0]
3348	else:
3349	MagSytSym += m.split('(')[0]+"'"
3350	if '2/m' in SytSym and '2' in m:
3351	MagSytSym += '/'
3352	if '-3/m' in SytSym:
3353	MagSytSym = '-'+MagSytSym
3354
3355	MagSytSym += axis
3356	# some exceptions & special rules
3357	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3358	return MagSytSym
3359
3360	# if len(GenSym) == 3:
3361	# if SGSpin[1] < 0:
3362	# if 'mm2' in SytSym:
3363	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3364	# else: #bad rule for I41/a
3365	# MagSytSym = SplitSytSym[0]+"'"
3366	# if len(SplitSytSym) > 1:
3367	# MagSytSym += '('+SplitSytSym[1]
3368	# else:
3369	# MagSytSym = SytSym
3370	# if len(SplitSytSym) >1:
3371	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3372	# MagSytSym = MagSytSym.replace('-4',"-4'")
3373	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3374	# MagSytSym = MagSytSym.replace('-6',"-6'")
3375	# return MagSytSym
3376	#
3377	return SytSym
3378
3379	def UpdateSytSym(Phase):
3380	''' Update site symmetry/site multiplicity after space group/BNS lattice change
3381	'''
3382	generalData = Phase['General']
3383	SGData = generalData['SGData']
3384	Atoms = Phase['Atoms']
3385	cx,ct,cs,cia = generalData['AtomPtrs']
3386	for atom in Atoms:
3387	XYZ = atom[cx:cx+3]
3388	sytsym,Mult = SytSym(XYZ,SGData)[:2]
3389	sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3390	atom[cs] = sytsym
3391	if generalData['Type'] == 'magnetic':
3392	magSytSym = MagSytSym(sytSym,dupDir,SGData)
3393	atom[cs] = magSytSym
3394	atom[cs+1] = Mult
3395	return
3396
3397	def ElemPosition(SGData):
3398	''' Under development.
3399	Object here is to return a list of symmetry element types and locations suitable
3400	for say drawing them.
3401	So far I have the element type... getting all possible locations without lookup may be impossible!
3402	'''
3403	Inv = SGData['SGInv']
3404	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3405	# get operators & expand if centrosymmetric
3406	Ops = SGData['SGOps']
3407	opM = np.array([op[0].T for op in Ops])
3408	opT = np.array([op[1] for op in Ops])
3409	if Inv:
3410	opM = np.concatenate((opM,-opM))
3411	opT = np.concatenate((opT,-opT))
3412	opMT = list(zip(opM,opT))
3413	for M,T in opMT[1:]: #skip I
3414	Dt = int(nl.det(M))
3415	Tr = int(np.trace(M))
3416	Dt = -(Dt-1)//2
3417	Es = eleSym[Tr][Dt]
3418	if Dt: #rotation-inversion
3419	I = np.eye(3)
3420	if Tr == 1: #mirrors/glides
3421	if np.any(T): #glide
3422	M2 = np.inner(M,M)
3423	MT = np.inner(M,T)+T
3424	print ('glide',Es,MT)
3425	print (M2)
3426	else: #mirror
3427	print ('mirror',Es,T)
3428	print (I-M)
3429	X = [-1,-1,-1]
3430	elif Tr == -3: # pure inversion
3431	X = np.inner(nl.inv(I-M),T)
3432	print ('inversion',Es,X)
3433	else: #other rotation-inversion
3434	M2 = np.inner(M,M)
3435	MT = np.inner(M,T)+T
3436	print ('rot-inv',Es,MT)
3437	print (M2)
3438	X = [-1,-1,-1]
3439	else: #rotations
3440	print ('rotation',Es)
3441	X = [-1,-1,-1]
3442	#SymElements.append([Es,X])
3443
3444	return #SymElements
3445
3446	def ApplyStringOps(A,SGData,X,Uij=[]):
3447	'Needs a doc string'
3448	SGOps = SGData['SGOps']
3449	SGCen = SGData['SGCen']
3450	Ax = A.split('+')
3451	Ax[0] = int(Ax[0])
3452	iC = 1
3453	if Ax[0] < 0:
3454	iC = -1
3455	Ax[0] = abs(Ax[0])
3456	nA = Ax[0]%100-1
3457	cA = Ax[0]//100
3458	Cen = SGCen[cA]
3459	M,T = SGOps[nA]
3460	if len(Ax)>1:
3461	cellA = Ax[1].split(',')
3462	cellA = np.array([int(a) for a in cellA])
3463	else:
3464	cellA = np.zeros(3)
3465	newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3466	if len(Uij):
3467	U = Uij2U(Uij)
3468	U = np.inner(M,np.inner(U,M).T)
3469	newUij = U2Uij(U)
3470	return [newX,newUij]
3471	else:
3472	return newX
3473
3474	def ApplyStringOpsMom(A,SGData,Mom):
3475	'Needs a doc string'
3476	SGOps = SGData['SGOps']
3477	Ax = A.split('+')
3478	Ax[0] = int(Ax[0])
3479	iAx = abs(Ax[0])
3480	nA = iAx%100-1
3481	nC = len(SGOps)*(iAx//100)
3482	NA = nA
3483	if Ax[0] < 0:
3484	NA += len(SGOps)
3485	M,T = SGOps[nA]
3486	if SGData['SGGray']:
3487	newMom = -np.inner(Mom,M).Tnl.det(M)SGData['SpnFlp'][NA+nC]
3488	else:
3489	newMom = np.inner(Mom,M).Tnl.det(M)SGData['SpnFlp'][NA+nC]
3490	# print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom)
3491	return newMom
3492
3493	def StringOpsProd(A,B,SGData):
3494	"""
3495	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3496	where '-' indicates inversion, c(>0) is the cell centering operator,
3497	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3498	Should return resultant string - C. SGData - dictionary using entries:
3499
3500	* 'SGCen': cell centering vectors [0,0,0] at least
3501	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3502
3503	"""
3504	SGOps = SGData['SGOps']
3505	SGCen = SGData['SGCen']
3506	#1st split out the cell translation part & work on the operator parts
3507	Ax = A.split('+'); Bx = B.split('+')
3508	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3509	iC = 0
3510	if Ax[0]*Bx[0] < 0:
3511	iC = 1
3512	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3513	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3514	cA = Ax[0]//100; cB = Bx[0]//100
3515	Cen = (SGCen[cA]+SGCen[cB])%1.0
3516	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3517	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3518	Mc = np.inner(Ma,Mb.T)
3519	# print Ma,Mb,Mc
3520	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3521	# print Ta,Tb,Tc
3522	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3523	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3524	#now the cell translation part
3525	if len(Ax)>1:
3526	cellA = Ax[1].split(',')
3527	cellA = [int(a) for a in cellA]
3528	else:
3529	cellA = [0,0,0]
3530	if len(Bx)>1:
3531	cellB = Bx[1].split(',')
3532	cellB = [int(b) for b in cellB]
3533	else:
3534	cellB = [0,0,0]
3535	cellC = np.add(cellA,cellB)
3536	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3537	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3538	return C
3539
3540	def U2Uij(U):
3541	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3542	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3543
3544	def Uij2U(Uij):
3545	#returns the thermal motion tensor U from Uij as numpy array
3546	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3547
3548	def StandardizeSpcName(spcgroup):
3549	'''Accept a spacegroup name where spaces may have not been used
3550	in the names according to the GSAS convention (spaces between symmetry
3551	for each axis) and return the space group name as used in GSAS
3552	'''
3553	rspc = spcgroup.replace(' ','').upper()
3554	# deal with rhombohedral and hexagonal setting designations
3555	rhomb = ''
3556	if rspc[-1:] == 'R':
3557	rspc = rspc[:-1]
3558	rhomb = ' R'
3559	gray = ''
3560	if "1'" in rspc:
3561	gray = " 1'"
3562	rspc = rspc.replace("1'",'')
3563	rspc = rspc.replace("'",'')
3564	if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3565	rspc = rspc[:-1]
3566	if rspc[1:3] in ['M3','N3','A3','D3']: #fix cubic old style
3567	rspc.replace('3','-3')
3568	bns = -1
3569	try:
3570	bns = rspc.index('_')
3571	rspc = rspc.replace(rspc[bns:bns+2],'')
3572	except ValueError:
3573	pass
3574	# look for a match in the spacegroup lists
3575	for i in spglist.values():
3576	for spc in i:
3577	if rspc == spc.replace(' ','').upper():
3578	return spc+gray+rhomb
3579	# how about the post-2002 orthorhombic names?
3580	if rspc in sgequiv_2002_orthorhombic:
3581	return sgequiv_2002_orthorhombic[rspc]+gray
3582	else:
3583	# not found
3584	return ''
3585
3586	def SpaceGroupNumber(spcgroup):
3587	SGNo = -1
3588	SpcGp = StandardizeSpcName(spcgroup)
3589	if not SpcGp:
3590	return SGNo
3591	try:
3592	SGNo = spgbyNum.index(SpcGp)
3593	except ValueError:
3594	pass
3595	return SGNo
3596
3597
3598	spgbyNum = []
3599	'''Space groups indexed by number'''
3600	spgbyNum = [None,
3601	'P 1','P -1', #1-2
3602	'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3603	'C 2/m','P 2/c','P 21/c','C 2/c', #3-15
3604	'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3605	'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3606	'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3607	'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3608	'C m m 2','C m c 21','C c c 2',
3609	'A m m 2','A b m 2','A m a 2','A b a 2',
3610	'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3611	'P m m m','P n n n','P c c m','P b a n',
3612	'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3613	'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3614	'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3615	'F m m m', 'F d d d',
3616	'I m m m','I b a m','I b c a','I m m a', #16-74
3617	'P 4','P 41','P 42','P 43',
3618	'I 4','I 41',
3619	'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3620	'I 4/m','I 41/a',
3621	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3622	'P 42 21 2','P 43 2 2','P 43 21 2',
3623	'I 4 2 2','I 41 2 2',
3624	'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3625	'P 42 m c','P 42 b c',
3626	'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3627	'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3628	'P -4 c 2','P -4 b 2','P -4 n 2',
3629	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3630	'P 4/m m m','P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m',
3631	'P 4/m n c','P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
3632	'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m','P 42/n m c',
3633	'P 42/n c m',
3634	'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3635	'P 3','P 31','P 32','R 3','P -3','R -3',
3636	'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3637	'R 3 2',
3638	'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3639	'R 3 m','R 3 c',
3640	'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3641	'R -3 m','R -3 c', #75-167
3642	'P 6','P 61',
3643	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3644	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3645	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3646	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c', #168-194
3647	'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3648	'F m -3','F d -3','I m -3',
3649	'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3650	'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3651	'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3652	'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3653	'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3654	'I m -3 m','I a -3 d',] #195-230
3655	altSettingOrtho = {}
3656	''' A dictionary of alternate settings for orthorhombic unit cells
3657	'''
3658	altSettingOrtho = {
3659	'P 2 2 2' :{'abc':'P 2 2 2','cab':'P 2 2 2','bca':'P 2 2 2','acb':'P 2 2 2','bac':'P 2 2 2','cba':'P 2 2 2'},
3660	'P 2 2 21' :{'abc':'P 2 2 21','cab':'P 21 2 2','bca':'P 2 21 2','acb':'P 2 21 2','bac':'P 2 2 21','cba':'P 21 2 2'},
3661	'P 21 21 2':{'abc':'P 21 21 2','cab':'P 2 21 21','bca':'P 21 2 21','acb':'P 21 2 21','bac':'P 21 21 2','cba':'P 2 21 21'},
3662	'P 21 21 21':{'abc':'P 21 21 21','cab':'P 21 21 21','bca':'P 21 21 21','acb':'P 21 21 21','bac':'P 21 21 21','cba':'P 21 21 21'},
3663	'C 2 2 21':{'abc':'C 2 2 21','cab':'A 21 2 2','bca':'B 2 21 2','acb':'B 2 21 2','bac':'C 2 2 21','cba':'A 21 2 2'},
3664	'C 2 2 2':{'abc':'C 2 2 2','cab':'A 2 2 2','bca':'B 2 2 2','acb':'B 2 2 2','bac':'C 2 2 2','cba':'A 2 2 2'},
3665	'F 2 2 2':{'abc':'F 2 2 2','cab':'F 2 2 2','bca':'F 2 2 2','acb':'F 2 2 2','bac':'F 2 2 2','cba':'F 2 2 2'},
3666	'I 2 2 2':{'abc':'I 2 2 2','cab':'I 2 2 2','bca':'I 2 2 2','acb':'I 2 2 2','bac':'I 2 2 2','cba':'I 2 2 2'},
3667	'I 21 21 21':{'abc':'I 21 21 21','cab':'I 21 21 21','bca':'I 21 21 21','acb':'I 21 21 21','bac':'I 21 21 21','cba':'I 21 21 21'},
3668	'P m m 2':{'abc':'P m m 2','cab':'P 2 m m','bca':'P m 2 m','acb':'P m 2 m','bac':'P m m 2','cba':'P 2 m m'},
3669	'P m c 21':{'abc':'P m c 21','cab':'P 21 m a','bca':'P b 21 m','acb':'P m 21 b','bac':'P c m 21','cba':'P 21 a m'},
3670	'P c c 2':{'abc':'P c c 2','cab':'P 2 a a','bca':'P b 2 b','acb':'P b 2 b','bac':'P c c 2','cba':'P 2 a a'},
3671	'P m a 2':{'abc':'P m a 2','cab':'P 2 m b','bca':'P c 2 m','acb':'P m 2 a','bac':'P b m 2','cba':'P 2 c m'},
3672	'P c a 21':{'abc':'P c a 21','cab':'P 21 a b','bca':'P c 21 b','acb':'P b 21 a','bac':'P b c 21','cba':'P 21 c a'},
3673	'P n c 2':{'abc':'P n c 2','cab':'P 2 n a','bca':'P b 2 n','acb':'P n 2 b','bac':'P c n 2','cba':'P 2 a n'},
3674	'P m n 21':{'abc':'P m n 21','cab':'P 21 m n','bca':'P n 21 m','acb':'P m 21 n','bac':'P n m 21','cba':'P 21 n m'},
3675	'P b a 2':{'abc':'P b a 2','cab':'P 2 c b','bca':'P c 2 a','acb':'P c 2 a','bac':'P b a 2','cba':'P 2 c b'},
3676	'P n a 21':{'abc':'P n a 21','cab':'P 21 n b','bca':'P c 21 n','acb':'P n 21 a','bac':'P b n 21','cba':'P 21 c n'},
3677	'P n n 2':{'abc':'P n n 2','cab':'P 2 n n','bca':'P n 2 n','acb':'P n 2 n','bac':'P n n 2','cba':'P 2 n n'},
3678	'C m m 2':{'abc':'C m m 2','cab':'A 2 m m','bca':'B m 2 m','acb':'B m 2 m','bac':'C m m 2','cba':'A 2 m m'},
3679	'C m c 21':{'abc':'C m c 21','cab':'A 21 m a','bca':'B b 21 m','acb':'B m 21 b','bac':'C c m 21','cba':'A 21 a m'},
3680	'C c c 2':{'abc':'C c c 2','cab':'A 2 a a','bca':'B b 2 b','acb':'B b 2 b','bac':'C c c 2','cba':'A 2 a a'},
3681	'A m m 2':{'abc':'A m m 2','cab':'B 2 m m','bca':'C m 2 m','acb':'A m 2 m','bac':'B m m 2','cba':'C 2 m m'},
3682	'A b m 2':{'abc':'A b m 2','cab':'B 2 c m','bca':'C m 2 a','acb':'A c 2 m','bac':'B m a 2','cba':'C 2 m b'},
3683	'A m a 2':{'abc':'A m a 2','cab':'B 2 m b','bca':'C c 2 m','acb':'A m 2 a','bac':'B b m 2','cba':'C 2 c m'},
3684	'A b a 2':{'abc':'A b a 2','cab':'B 2 c b','bca':'C c 2 a','acb':'A c 2 a','bac':'B b a 2','cba':'C 2 c b'},
3685	'F m m 2':{'abc':'F m m 2','cab':'F 2 m m','bca':'F m 2 m','acb':'F m 2 m','bac':'F m m 2','cba':'F 2 m m'},
3686	'F d d 2':{'abc':'F d d 2','cab':'F 2 d d','bca':'F d 2 d','acb':'F d 2 d','bac':'F d d 2','cba':'F 2 d d'},
3687	'I m m 2':{'abc':'I m m 2','cab':'I 2 m m','bca':'I m 2 m','acb':'I m 2 m','bac':'I m m 2','cba':'I 2 m m'},
3688	'I b a 2':{'abc':'I b a 2','cab':'I 2 c b','bca':'I c 2 a','acb':'I c 2 a','bac':'I b a 2','cba':'I 2 c b'},
3689	'I m a 2':{'abc':'I m a 2','cab':'I 2 m b','bca':'I c 2 m','acb':'I m 2 a','bac':'I b m 2','cba':'I 2 c m'},
3690	'P m m m':{'abc':'P m m m','cab':'P m m m','bca':'P m m m','acb':'P m m m','bac':'P m m m','cba':'P m m m'},
3691	'P n n n':{'abc':'P n n n','cab':'P n n n','bca':'P n n n','acb':'P n n n','bac':'P n n n','cba':'P n n n'},
3692	'P c c m':{'abc':'P c c m','cab':'P m a a','bca':'P b m b','acb':'P b m b','bac':'P c c m','cba':'P m a a'},
3693	'P b a n':{'abc':'P b a n','cab':'P n c b','bca':'P c n a','acb':'P c n a','bac':'P b a n','cba':'P n c b'},
3694	'P m m a':{'abc':'P m m a','cab':'P b m m','bca':'P m c m','acb':'P m a m','bac':'P m m b','cba':'P c m m'},
3695	'P n n a':{'abc':'P n n a','cab':'P b n n','bca':'P n c n','acb':'P n a n','bac':'P n n b','cba':'P c n n'},
3696	'P m n a':{'abc':'P m n a','cab':'P b m n','bca':'P n c m','acb':'P m a n','bac':'P n m b','cba':'P c n m'},
3697	'P c c a':{'abc':'P c c a','cab':'P b a a','bca':'P b c b','acb':'P b a b','bac':'P c c b','cba':'P c a a'},
3698	'P b a m':{'abc':'P b a m','cab':'P m c b','bca':'P c m a','acb':'P c m a','bac':'P b a m','cba':'P m c b'},
3699	'P c c n':{'abc':'P c c n','cab':'P n a a','bca':'P b n b','acb':'P b n b','bac':'P c c n','cba':'P n a a'},
3700	'P b c m':{'abc':'P b c m','cab':'P m c a','bca':'P b m a','acb':'P c m b','bac':'P c a m','cba':'P m a b'},
3701	'P n n m':{'abc':'P n n m','cab':'P m n n','bca':'P n m n','acb':'P n m n','bac':'P n n m','cba':'P m n n'},
3702	'P m m n':{'abc':'P m m n','cab':'P n m m','bca':'P m n m','acb':'P m n m','bac':'P m m n','cba':'P n m m'},
3703	'P b c n':{'abc':'P b c n','cab':'P n c a','bca':'P b n a','acb':'P c n b','bac':'P c a n','cba':'P n a b'},
3704	'P b c a':{'abc':'P b c a','cab':'P b c a','bca':'P b c a','acb':'P c a b','bac':'P c a b','cba':'P c a b'},
3705	'P n m a':{'abc':'P n m a','cab':'P b n m','bca':'P m c n','acb':'P n a m','bac':'P m n b','cba':'P c m n'},
3706	'C m c m':{'abc':'C m c m','cab':'A m m a','bca':'B b m m','acb':'B m m b','bac':'C c m m','cba':'A m a m'},
3707	'C m c a':{'abc':'C m c a','cab':'A b m a','bca':'B b c m','acb':'B m a b','bac':'C c m b','cba':'A c a m'},
3708	'C m m m':{'abc':'C m m m','cab':'A m m m','bca':'B m m m','acb':'B m m m','bac':'C m m m','cba':'A m m m'},
3709	'C c c m':{'abc':'C c c m','cab':'A m a a','bca':'B b m b','acb':'B b m b','bac':'C c c m','cba':'A m a a'},
3710	'C m m a':{'abc':'C m m a','cab':'A b m m','bca':'B m c m','acb':'B m a m','bac':'C m m b','cba':'A c m m'},
3711	'C c c a':{'abc':'C c a a','cab':'A b a a','bca':'B b c b','acb':'B b a b','bac':'C c c b','cba':'A c a a'},
3712	'F m m m':{'abc':'F m m m','cab':'F m m m','bca':'F m m m','acb':'F m m m','bac':'F m m m','cba':'F m m m'},
3713	'F d d d':{'abc':'F d d d','cab':'F d d d','bca':'F d d d','acb':'F d d d','bac':'F d d d','cba':'F d d d'},
3714	'I m m m':{'abc':'I m m m','cab':'I m m m','bca':'I m m m','acb':'I m m m','bac':'I m m m','cba':'I m m m'},
3715	'I b a m':{'abc':'I b a m','cab':'I m c b','bca':'I c m a','acb':'I c m a','bac':'I b a m','cba':'I m c b'},
3716	'I b c a':{'abc':'I b c a','cab':'I b c a','bca':'I b c a','acb':'I c a b','bac':'I c a b','cba':'I c a b'},
3717	'I m m a':{'abc':'I m m a','cab':'I b m m','bca':'I m c m','acb':'I m a m','bac':'I m m b','cba':'I c m m'},
3718	}
3719	spg2origins = {}
3720	''' A dictionary of all spacegroups that have 2nd settings; the value is the
3721	1st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3722	'''
3723	spg2origins = {
3724	'P n n n':[-.25,-.25,-.25],
3725	'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3726	'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3727	'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3728	'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3729	'F d d d':[-.125,-.125,-.125],
3730	'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3731	'P 4/n b m':[-.25,-.25,0],'P 4/n n c':[-.25,-.25,-.25],'P 4/n m m':[-.25,-.25,0],'P 4/n c c':[-.25,-.25,0],
3732	'P 42/n b c':[-.25,-.25,-.25],'P 42/n n m':[-.25,.25,-.25],'P 42/n m c':[-.25,.25,-.25],'P 42/n c m':[-.25,.25,-.25],
3733	'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3734	'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3735	'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3736	'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3737	'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3738	spglist = {}
3739	'''A dictionary of space groups as ordered and named in the pre-2002 International
3740	Tables Volume A, except that spaces are used following the GSAS convention to
3741	separate the different crystallographic directions.
3742	Note that the symmetry codes here will recognize many non-standard space group
3743	symbols with different settings. They are ordered by Laue group
3744	'''
3745	spglist = {
3746	'P1' : ('P 1','P -1',), # 1-2
3747	'C1' : ('C 1','C -1',),
3748	'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3749	'P 2/m','P 21/m','P 2/c','P 2/a','P 2/n','P 21/c','P 21/a','P 21/n',), #3-15
3750	'C2/m':('C 2','C m','C c','C n',
3751	'C 2/m','C 2/c','C 2/n',),
3752	'A2/m':('A 2','A m','A a','A n',
3753	'A 2/m','A 2/a','A 2/n',),
3754	'I2/m':('I 2','I m','I a','I n','I c',
3755	'I 2/m','I 2/a','I 2/c','I 2/n',),
3756	'Pmmm':('P 2 2 2',
3757	'P 2 2 21','P 21 2 2','P 2 21 2',
3758	'P 21 21 2','P 2 21 21','P 21 2 21',
3759	'P 21 21 21',
3760	'P m m 2','P 2 m m','P m 2 m',
3761	'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3762	'P c c 2','P 2 a a','P b 2 b',
3763	'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3764	'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3765	'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3766	'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3767	'P b a 2','P 2 c b','P c 2 a',
3768	'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3769	'P n n 2','P 2 n n','P n 2 n',
3770	'P m m m','P n n n',
3771	'P c c m','P m a a','P b m b',
3772	'P b a n','P n c b','P c n a',
3773	'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3774	'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3775	'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3776	'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3777	'P b a m','P m c b','P c m a',
3778	'P c c n','P n a a','P b n b',
3779	'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3780	'P n n m','P m n n','P n m n',
3781	'P m m n','P n m m','P m n m',
3782	'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3783	'P b c a','P c a b',
3784	'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3785	),
3786	'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3787	'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3788	'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3789	'C m c m','C c m m','C m c a','C c m b',
3790	'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3791	'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3792	'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3793	'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3794	'A m m a','A m a m','A b m a','A c a m',
3795	'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3796	'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3797	'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3798	'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3799	'B b m m','B m m b','B b c m','B m a b',
3800	'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3801	'Immm':('I 2 2 2','I 21 21 21',
3802	'I m m 2','I m 2 m','I 2 m m',
3803	'I b a 2','I 2 c b','I c 2 a',
3804	'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3805	'I m m m','I b a m','I m c b','I c m a',
3806	'I b c a','I c a b',
3807	'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3808	'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3809	'F m m 2','F m 2 m','F 2 m m',
3810	'F d d 2','F d 2 d','F 2 d d',),
3811	'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3812	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3813	'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3814	'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3815	'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3816	'P -4 n 2','P 4/m m m','P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m',
3817	'P 4/m n c','P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
3818	'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m','P 42/n m c',
3819	'P 42/n c m',),
3820	'I4/mmm':('I 4','I 41','I -4','I 4/m','I 41/a','I 4 2 2','I 41 2 2','I 4 m m',
3821	'I 4 c m','I 41 m d','I 41 c d',
3822	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d','I 4/m m m','I 4/m c m',
3823	'I 41/a m d','I 41/a c d'),
3824	'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3825	'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3826	'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3827	'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3828	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3829	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3830	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3831	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3832	'Pm3m': ('P 2 3','P 21 3','P m 3','P m -3','P n 3','P n -3','P a 3','P a -3',
3833	'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3834	'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3835	'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3836	'Im3m':('I 2 3','I 21 3','I m 3','I m -3','I a 3','I a -3', 'I 4 3 2','I 41 3 2',
3837	'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a 3 d','I a -3 d','I n 3 n','I n -3 n'),
3838	'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3839	'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3840	'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3841	'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3842	}
3843	sgequiv_2002_orthorhombic = {}
3844	''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3845	along with the pre-2002 name. The e designates a double glide-plane
3846	'''
3847	sgequiv_2002_orthorhombic = {
3848	'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3849	'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3850	'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3851	'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3852	'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3853	'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3854	'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3855	'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3856
3857	#'A few non-standard space groups for test use'
3858	nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h',
3859	'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3860	'R 3 c r','R -3 c r','R -3 m r',),
3861
3862	#Use the space groups types in this order to list the symbols in the
3863	#order they are listed in the International Tables, vol. A'''
3864	symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal',
3865	'trigonal', 'hexagonal', 'cubic')
3866
3867	# self-test materials follow. Requires files in directory testinp
3868	selftestlist = []
3869	'''Defines a list of self-tests'''
3870	selftestquiet = True
3871	def _ReportTest():
3872	'Report name and doc string of current routine when ``selftestquiet`` is False'
3873	if not selftestquiet:
3874	import inspect
3875	caller = inspect.stack()[1][3]
3876	doc = eval(caller).__doc__
3877	if doc is not None:
3878	print('testing '+__file__+' with '+caller+' ('+doc+')')
3879	else:
3880	print('testing '+__file__()+" with "+caller)
3881	def test0():
3882	'''self-test #0: exercise MoveToUnitCell'''
3883	_ReportTest()
3884	msg = "MoveToUnitCell failed"
3885	assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3886	assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3887	assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3888	assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3889	selftestlist.append(test0)
3890
3891	def test1():
3892	'''self-test #1: SpcGroup against previous results'''
3893	#'''self-test #1: SpcGroup and SGPrint against previous results'''
3894	_ReportTest()
3895	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3896	if ospath.exists(testdir):
3897	if testdir not in sys.path: sys.path.insert(0,testdir)
3898	import spctestinp
3899	def CompareSpcGroup(spc, referr, refdict, reflist):
3900	'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3901	# if an error is reported, the dictionary can be ignored
3902	msg0 = "CompareSpcGroup failed on space group %s" % spc
3903	result = SpcGroup(spc)
3904	if result[0] == referr and referr > 0: return True
3905	# #print result[1]['SpGrp']
3906	#msg = msg0 + " in list lengths"
3907	#assert len(keys) == len(refdict.keys()), msg
3908	for key in refdict.keys():
3909	if key == 'SGOps' or key == 'SGCen':
3910	msg = msg0 + (" in key %s length" % key)
3911	assert len(refdict[key]) == len(result[1][key]), msg
3912	for i in range(len(refdict[key])):
3913	msg = msg0 + (" in key %s level 0" % key)
3914	assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
3915	msg = msg0 + (" in key %s level 1" % key)
3916	assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
3917	else:
3918	msg = msg0 + (" in key %s" % key)
3919	assert result[1][key] == refdict[key], msg
3920	msg = msg0 + (" in key %s reflist" % key)
3921	#for (l1,l2) in zip(reflist, SGPrint(result[1])):
3922	# assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
3923	# for now disable SGPrint testing, output has changed
3924	#assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
3925	for spc in spctestinp.SGdat:
3926	CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
3927	selftestlist.append(test1)
3928
3929	def test2():
3930	'''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
3931	_ReportTest()
3932	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3933	if ospath.exists(testdir):
3934	if testdir not in sys.path: sys.path.insert(0,testdir)
3935	import sgtbxtestinp
3936	def CompareWcctbx(spcname, cctbx_in, debug=0):
3937	'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
3938	cctbx = cctbx_in[:] # make copy so we don't delete from the original
3939	spc = (SpcGroup(spcname))[1]
3940	if debug: print (spc['SpGrp'])
3941	if debug: print (spc['SGCen'])
3942	latticetype = spcname.strip().upper()[0]
3943	# lattice type of R implies Hexagonal centering", fix the rhombohedral settings
3944	if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
3945	assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
3946	onebar = [1]
3947	if spc['SGInv']: onebar.append(-1)
3948	for (op,off) in spc['SGOps']:
3949	for inv in onebar:
3950	for cen in spc['SGCen']:
3951	noff = off + cen
3952	noff = MoveToUnitCell(noff)[0]
3953	mult = tuple((op*inv).ravel().tolist())
3954	if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
3955	for refop in cctbx:
3956	if debug: print (refop)
3957	# check the transform
3958	if refop[:9] != mult: continue
3959	if debug: print ("mult match")
3960	# check the translation
3961	reftrans = list(refop[-3:])
3962	reftrans = MoveToUnitCell(reftrans)[0]
3963	if all(abs(noff - reftrans) < 1.e-5):
3964	cctbx.remove(refop)
3965	break
3966	else:
3967	assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
3968	for key in sgtbxtestinp.sgtbx:
3969	CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
3970	selftestlist.append(test2)
3971
3972	def test3():
3973	'''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
3974	for selected space groups against info in IT Volume A '''
3975	_ReportTest()
3976	def ExerciseSiteSym (spc, crdlist):
3977	'compare site symmetries and multiplicities for a specified space group'
3978	msg = "failed on site sym test for %s" % spc
3979	(E,S) = SpcGroup(spc)
3980	assert not E, msg
3981	for t in crdlist:
3982	symb, m, n, od = SytSym(t[0],S)
3983	if symb.strip() != t[2].strip() or m != t[1]:
3984	print (spc,t[0],m,n,symb,t[2],od)
3985	assert m == t[1]
3986	#assert symb.strip() == t[2].strip()
3987
3988	ExerciseSiteSym('p 1',[
3989	((0.13,0.22,0.31),1,'1'),
3990	((0,0,0),1,'1'),
3991	])
3992	ExerciseSiteSym('p -1',[
3993	((0.13,0.22,0.31),2,'1'),
3994	((0,0.5,0),1,'-1'),
3995	])
3996	ExerciseSiteSym('C 2/c',[
3997	((0.13,0.22,0.31),8,'1'),
3998	((0.0,.31,0.25),4,'2(y)'),
3999	((0.25,.25,0.5),4,'-1'),
4000	((0,0.5,0),4,'-1'),
4001	])
4002	ExerciseSiteSym('p 2 2 2',[
4003	((0.13,0.22,0.31),4,'1'),
4004	((0,0.5,.31),2,'2(z)'),
4005	((0.5,.31,0.5),2,'2(y)'),
4006	((.11,0,0),2,'2(x)'),
4007	((0,0.5,0),1,'222'),
4008	])
4009	ExerciseSiteSym('p 4/n',[
4010	((0.13,0.22,0.31),8,'1'),
4011	((0.25,0.75,.31),4,'2(z)'),
4012	((0.5,0.5,0.5),4,'-1'),
4013	((0,0.5,0),4,'-1'),
4014	((0.25,0.25,.31),2,'4(001)'),
4015	((0.25,.75,0.5),2,'-4(001)'),
4016	((0.25,.75,0.0),2,'-4(001)'),
4017	])
4018	ExerciseSiteSym('p 31 2 1',[
4019	((0.13,0.22,0.31),6,'1'),
4020	((0.13,0.0,0.833333333),3,'2(100)'),
4021	((0.13,0.13,0.),3,'2(110)'),
4022	])
4023	ExerciseSiteSym('R 3 c',[
4024	((0.13,0.22,0.31),18,'1'),
4025	((0.0,0.0,0.31),6,'3'),
4026	])
4027	ExerciseSiteSym('R 3 c R',[
4028	((0.13,0.22,0.31),6,'1'),
4029	((0.31,0.31,0.31),2,'3(111)'),
4030	])
4031	ExerciseSiteSym('P 63 m c',[
4032	((0.13,0.22,0.31),12,'1'),
4033	((0.11,0.22,0.31),6,'m(100)'),
4034	((0.333333,0.6666667,0.31),2,'3m(100)'),
4035	((0,0,0.31),2,'3m(100)'),
4036	])
4037	ExerciseSiteSym('I a -3',[
4038	((0.13,0.22,0.31),48,'1'),
4039	((0.11,0,0.25),24,'2(x)'),
4040	((0.11,0.11,0.11),16,'3(111)'),
4041	((0,0,0),8,'-3(111)'),
4042	])
4043	selftestlist.append(test3)
4044
4045	if __name__ == '__main__':
4046	# run self-tests
4047	selftestquiet = False
4048	for test in selftestlist:
4049	test()
4050	print ("OK")

Note: See TracBrowser for help on using the repository browser.

Download in other formats: