Context Navigation

source: trunk/GSASIIspc.py @ 3625

Last change on this file since 3625 was 3625, checked in by vondreele, 7 years ago
fix 2 bugs in mag. space group stuff correct Mag Structure tutorial name
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 170.2 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-09-27 17:30:16 +0000 (Thu, 27 Sep 2018) $
12	# $Author: vondreele $
13	# $Revision: 3625 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 3625 2018-09-27 17:30:16Z 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: 3625 $")
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]}
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]})
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	dF = None
2560	dFTP = None
2561	if siteSym == '1':
2562	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2563	elif siteSym == '-1':
2564	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2565	else:
2566	delt2 = np.eye(2)*0.001
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	dX,dXTP = None,None
2618	if siteSym == '1':
2619	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2620	elif siteSym == '-1':
2621	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.],]
2622	else:
2623	delt4 = np.ones(4)*0.001
2624	delt5 = np.ones(5)*0.001
2625	delt6 = np.eye(6)*0.001
2626	if 'Fourier' in waveType:
2627	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2628	#3x6x12 modulated position array (X,Spos,tau)& force positive
2629	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2630	elif waveType == 'Sawtooth':
2631	dX = posSawtooth(tau,delt4[0],delt4[1:])
2632	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0]]),
2633	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2634	elif waveType in ['ZigZag','Block']:
2635	if waveType == 'ZigZag':
2636	dX = posZigZag(tau,delt5[:2],delt5[2:])
2637	else:
2638	dX = posBlock(tau,delt5[:2],delt5[2:])
2639	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2640	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2641	XSC = np.ones(6,dtype='i')
2642	dXTP = []
2643	for i in SdIndx:
2644	sop = Sop[i]
2645	ssop = SSop[i]
2646	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2647	xsc = np.ones(6,dtype='i')
2648	if 'Fourier' in waveType:
2649	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2650	elif waveType == 'Sawtooth':
2651	dXT = posSawtooth(tauT,delt4[0],delt4[1:])+np.array(XYZ)[:,nxs,nxs]
2652	elif waveType == 'ZigZag':
2653	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2654	elif waveType == 'Block':
2655	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2656	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2657	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2658	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2659	dXTP.append(dXT)
2660	if waveType == 'Fourier':
2661	for i in range(3):
2662	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2663	xsc[i] = 0
2664	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2665	xsc[i+3] = 0
2666	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2667	xsc[3:6] = 0
2668	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2669	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2670	if dT:
2671	if '(x)' in siteSym:
2672	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2673	if 'm' in siteSym and len(SdIndx) == 1:
2674	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2675	elif '(y)' in siteSym:
2676	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2677	if 'm' in siteSym and len(SdIndx) == 1:
2678	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2679	elif '(z)' in siteSym:
2680	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2681	if 'm' in siteSym and len(SdIndx) == 1:
2682	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2683	else:
2684	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2685	if '4/mmm' in laue:
2686	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2687	if '(xy)' in siteSym:
2688	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2689	if dT:
2690	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2691	else:
2692	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2693	if '(xy)' in siteSym or '(+-0)' in siteSym:
2694	mul = 1
2695	if '(+-0)' in siteSym:
2696	mul = -1
2697	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2698	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2699	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2700	xsc[3:5] = 0
2701	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2702	CSI[0][:2] = [[12,0,0],[12,0,0]]
2703	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2704	xsc[:2] = 0
2705	XSC &= xsc
2706	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2707	if waveType == 'Fourier':
2708	n = -1
2709	if debug: print (XSC)
2710	for i,X in enumerate(XSC):
2711	if X:
2712	n += 1
2713	CSI[0][i][0] = n+1
2714	CSI[1][i][0] = 1.0
2715
2716	return list(CSI),dX,dXTP
2717
2718	def DoUij():
2719	dU,dUTP = None,None
2720	if siteSym == '1':
2721	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2722	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2723	elif siteSym == '-1':
2724	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2725	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2726	else:
2727	tau = np.linspace(0,1,49,True)
2728	delt12 = np.eye(12)*0.0001
2729	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2730	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2731	USC = np.ones(12,dtype='i')
2732	dUTP = []
2733	dtau = 0.
2734	for i in SdIndx:
2735	sop = Sop[i]
2736	ssop = SSop[i]
2737	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2738	usc = np.ones(12,dtype='i')
2739	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2740	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2741	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2742	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2743	dUT = dUT[:,:,np.argsort(tauT)]
2744	dUT[:,:6,:] =(ssdetsdet)
2745	dUTP.append(dUT)
2746	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2747	if dT:
2748	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2749	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2750	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2751	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2752	else:
2753	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2754	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2755	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2756	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2757	if 'mm2(x)' in siteSym and dT:
2758	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2759	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2760	elif '(xy)' in siteSym and dT:
2761	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2762	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2763	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2764	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2765	elif '(x)' in siteSym and dT:
2766	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2767	elif '(y)' in siteSym and dT:
2768	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2769	elif '(z)' in siteSym and dT:
2770	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2771	for i in range(6):
2772	if not USC[i]:
2773	CSI[0][i] = [0,0,0]
2774	CSI[1][i] = [0.,0.,0.]
2775	CSI[0][i+6] = [0,0,0]
2776	CSI[1][i+6] = [0.,0.,0.]
2777	else:
2778	for i in range(6):
2779	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2780	usc[i] = 0
2781	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2782	usc[i+6] = 0
2783	if np.any(dUT[1,0,:]):
2784	if '4/m' in siteSym:
2785	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2786	if ssop[1][3]:
2787	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2788	usc[9] = 1
2789	else:
2790	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2791	usc[9] = 0
2792	elif '4' in siteSym:
2793	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2794	CSI[0][:2] = [[11,0,0],[11,0,0]]
2795	if ssop[1][3]:
2796	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2797	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2798	usc[2] = 0
2799	usc[8] = 0
2800	usc[3] = 1
2801	usc[9] = 1
2802	else:
2803	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2804	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2805	usc[2] = 1
2806	usc[8] = 1
2807	usc[3] = 0
2808	usc[9] = 0
2809	elif 'xy' in siteSym or '+-0' in siteSym:
2810	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2811	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2812	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2813	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2814	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2815	usc[4:6] = 0
2816	usc[6:8] = 0
2817
2818	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2819	USC &= usc
2820	if debug: print (USC)
2821	if not np.any(dtau%.5):
2822	n = -1
2823	for i,U in enumerate(USC):
2824	if U:
2825	n += 1
2826	CSI[0][i][0] = n+1
2827	CSI[1][i][0] = 1.0
2828
2829	return list(CSI),dU,dUTP
2830
2831	def DoMag():
2832	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2833	if siteSym == '1':
2834	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2835	elif siteSym in ['-1','mmm','2/m(x)','2/m(y)','2/m(z)','4/mmm001']:
2836	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2837	else:
2838	delt6 = np.eye(6)*0.001
2839	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2840	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2841	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2842	MSC = np.ones(6,dtype='i')
2843	dMTP = []
2844	for i in SdIndx:
2845	sop = Sop[i]
2846	ssop = SSop[i]
2847	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2848	msc = np.ones(6,dtype='i')
2849	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2850	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2851	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2852	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2853	dMTP.append(dMT)
2854	for i in range(3):
2855	if not np.allclose(dM[i,i,:],dMT[i,i,:]):
2856	msc[i] = 0
2857	if not np.allclose(dM[i,i+3,:],dMT[i,i+3,:]):
2858	msc[i+3] = 0
2859	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2860	msc[3:6] = 0
2861	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2862	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2863	if dT:
2864	if '(x)' in siteSym:
2865	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2866	if 'm' in siteSym and len(SdIndx) == 1:
2867	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2868	elif '(y)' in siteSym:
2869	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2870	if 'm' in siteSym and len(SdIndx) == 1:
2871	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2872	elif '(z)' in siteSym:
2873	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2874	if 'm' in siteSym and len(SdIndx) == 1:
2875	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2876	else:
2877	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2878	if '4/mmm' in laue:
2879	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2880	if '(xy)' in siteSym:
2881	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2882	if dT:
2883	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2884	else:
2885	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2886	if '(xy)' in siteSym or '(+-0)' in siteSym:
2887	mul = 1
2888	if '(+-0)' in siteSym:
2889	mul = -1
2890	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2891	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2892	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2893	msc[3:5] = 0
2894	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2895	CSI[0][:2] = [[12,0,0],[12,0,0]]
2896	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2897	msc[:2] = 0
2898	MSC &= msc
2899	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2900	n = -1
2901	if debug: print (MSC)
2902	for i,M in enumerate(MSC):
2903	if M:
2904	n += 1
2905	CSI[0][i][0] = n+1
2906	CSI[1][i][0] = 1.0
2907
2908	return CSI,None,None
2909
2910	if debug: print ('super space group: '+SSGData['SSpGrp'])
2911	xyz = np.array(XYZ)%1.
2912	SGOps = copy.deepcopy(SGData['SGOps'])
2913	laue = SGData['SGLaue']
2914	siteSym = SytSym(XYZ,SGData)[0].strip()
2915	if debug: print ('siteSym: '+siteSym)
2916	SSGOps = copy.deepcopy(SSGData['SSGOps'])
2917	#expand ops to include inversions if any
2918	if SGData['SGInv']:
2919	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
2920	SGOps.append([-op[0],-op[1]%1.])
2921	SSGOps.append([-sop[0],-sop[1]%1.])
2922	#build set of sym ops around special position
2923	SSop = []
2924	Sop = []
2925	Sdtau = []
2926	for iop,Op in enumerate(SGOps):
2927	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
2928	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
2929	SSop.append(SSGOps[iop])
2930	Sop.append(SGOps[iop])
2931	ssopinv = nl.inv(SSGOps[iop][0])
2932	mst = ssopinv[3][:3]
2933	epsinv = ssopinv[3][3]
2934	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
2935	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
2936	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
2937	#setup displacement arrays
2938	tau = np.linspace(-1,1,49,True)
2939	#make modulation arrays - one parameter at a time
2940	if Stype == 'Sfrac':
2941	CSI,dF,dFTP = DoFrac()
2942	elif Stype == 'Spos':
2943	CSI,dF,dFTP = DoXYZ()
2944	CSI[0] = orderParms(CSI[0])
2945	elif Stype == 'Sadp':
2946	CSI,dF,dFTP = DoUij()
2947	CSI[0] = orderParms(CSI[0])
2948	elif Stype == 'Smag':
2949	CSI,dF,dFTP = DoMag()
2950
2951	if debug:
2952	return CSI,dF,dFTP
2953	else:
2954	return CSI
2955
2956	def MustrainNames(SGData):
2957	'Needs a doc string'
2958	laue = SGData['SGLaue']
2959	uniq = SGData['SGUniq']
2960	if laue in ['m3','m3m']:
2961	return ['S400','S220']
2962	elif laue in ['6/m','6/mmm','3m1']:
2963	return ['S400','S004','S202']
2964	elif laue in ['31m','3']:
2965	return ['S400','S004','S202','S211']
2966	elif laue in ['3R','3mR']:
2967	return ['S400','S220','S310','S211']
2968	elif laue in ['4/m','4/mmm']:
2969	return ['S400','S004','S220','S022']
2970	elif laue in ['mmm']:
2971	return ['S400','S040','S004','S220','S202','S022']
2972	elif laue in ['2/m']:
2973	SHKL = ['S400','S040','S004','S220','S202','S022']
2974	if uniq == 'a':
2975	SHKL += ['S013','S031','S211']
2976	elif uniq == 'b':
2977	SHKL += ['S301','S103','S121']
2978	elif uniq == 'c':
2979	SHKL += ['S130','S310','S112']
2980	return SHKL
2981	else:
2982	SHKL = ['S400','S040','S004','S220','S202','S022']
2983	SHKL += ['S310','S103','S031','S130','S301','S013']
2984	SHKL += ['S211','S121','S112']
2985	return SHKL
2986
2987	def HStrainVals(HSvals,SGData):
2988	laue = SGData['SGLaue']
2989	uniq = SGData['SGUniq']
2990	DIJ = np.zeros(6)
2991	if laue in ['m3','m3m']:
2992	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
2993	elif laue in ['6/m','6/mmm','3m1','31m','3']:
2994	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
2995	elif laue in ['3R','3mR']:
2996	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
2997	elif laue in ['4/m','4/mmm']:
2998	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
2999	elif laue in ['mmm']:
3000	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3001	elif laue in ['2/m']:
3002	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3003	if uniq == 'a':
3004	DIJ[5] = HSvals[3]
3005	elif uniq == 'b':
3006	DIJ[4] = HSvals[3]
3007	elif uniq == 'c':
3008	DIJ[3] = HSvals[3]
3009	else:
3010	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3011	return DIJ
3012
3013	def HStrainNames(SGData):
3014	'Needs a doc string'
3015	laue = SGData['SGLaue']
3016	uniq = SGData['SGUniq']
3017	if laue in ['m3','m3m']:
3018	return ['D11','eA'] #add cubic strain term
3019	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3020	return ['D11','D33']
3021	elif laue in ['3R','3mR']:
3022	return ['D11','D12']
3023	elif laue in ['4/m','4/mmm']:
3024	return ['D11','D33']
3025	elif laue in ['mmm']:
3026	return ['D11','D22','D33']
3027	elif laue in ['2/m']:
3028	Dij = ['D11','D22','D33']
3029	if uniq == 'a':
3030	Dij += ['D23']
3031	elif uniq == 'b':
3032	Dij += ['D13']
3033	elif uniq == 'c':
3034	Dij += ['D12']
3035	return Dij
3036	else:
3037	Dij = ['D11','D22','D33','D12','D13','D23']
3038	return Dij
3039
3040	def MustrainCoeff(HKL,SGData):
3041	'Needs a doc string'
3042	#NB: order of terms is the same as returned by MustrainNames
3043	laue = SGData['SGLaue']
3044	uniq = SGData['SGUniq']
3045	h,k,l = HKL
3046	Strm = []
3047	if laue in ['m3','m3m']:
3048	Strm.append(h4+k4+l**4)
3049	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3050	elif laue in ['6/m','6/mmm','3m1']:
3051	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3052	Strm.append(l**4)
3053	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3054	elif laue in ['31m','3']:
3055	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3056	Strm.append(l**4)
3057	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3058	Strm.append(4.0hkl(h+k))
3059	elif laue in ['3R','3mR']:
3060	Strm.append(h4+k4+l**4)
3061	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3062	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
3063	Strm.append(4.0(klh2+hlk2+hkl*2))
3064	elif laue in ['4/m','4/mmm']:
3065	Strm.append(h4+k4)
3066	Strm.append(l**4)
3067	Strm.append(3.0(hk)**2)
3068	Strm.append(3.0((hl)*2+(kl)**2))
3069	elif laue in ['mmm']:
3070	Strm.append(h**4)
3071	Strm.append(k**4)
3072	Strm.append(l**4)
3073	Strm.append(3.0(hk)**2)
3074	Strm.append(3.0(hl)**2)
3075	Strm.append(3.0(kl)**2)
3076	elif laue in ['2/m']:
3077	Strm.append(h**4)
3078	Strm.append(k**4)
3079	Strm.append(l**4)
3080	Strm.append(3.0(hk)**2)
3081	Strm.append(3.0(hl)**2)
3082	Strm.append(3.0(kl)**2)
3083	if uniq == 'a':
3084	Strm.append(2.0kl**3)
3085	Strm.append(2.0lk**3)
3086	Strm.append(4.0klh*2)
3087	elif uniq == 'b':
3088	Strm.append(2.0lh**3)
3089	Strm.append(2.0hl**3)
3090	Strm.append(4.0hlk*2)
3091	elif uniq == 'c':
3092	Strm.append(2.0hk**3)
3093	Strm.append(2.0kh**3)
3094	Strm.append(4.0hkl*2)
3095	else:
3096	Strm.append(h**4)
3097	Strm.append(k**4)
3098	Strm.append(l**4)
3099	Strm.append(3.0(hk)**2)
3100	Strm.append(3.0(hl)**2)
3101	Strm.append(3.0(kl)**2)
3102	Strm.append(2.0kh**3)
3103	Strm.append(2.0hl**3)
3104	Strm.append(2.0lk**3)
3105	Strm.append(2.0hk**3)
3106	Strm.append(2.0lh**3)
3107	Strm.append(2.0kl**3)
3108	Strm.append(4.0klh*2)
3109	Strm.append(4.0hlk*2)
3110	Strm.append(4.0khl*2)
3111	return Strm
3112
3113	def MuShklMean(SGData,Amat,Shkl):
3114
3115	def genMustrain(xyz,Shkl):
3116	uvw = np.inner(Amat.T,xyz)
3117	Strm = np.array(MustrainCoeff(uvw,SGData))
3118	Sum = np.sum(np.multiply(Shkl,Strm))
3119	Sum = np.where(Sum > 0.01,Sum,0.01)
3120	Sum = np.sqrt(Sum)
3121	return Sum*xyz
3122
3123	PHI = np.linspace(0.,360.,30,True)
3124	PSI = np.linspace(0.,180.,30,True)
3125	X = np.outer(npcosd(PHI),npsind(PSI))
3126	Y = np.outer(npsind(PHI),npsind(PSI))
3127	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3128	XYZ = np.dstack((X,Y,Z))
3129	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3130	return np.sqrt(np.sum(XYZ**2)/900.)
3131
3132	def Muiso2Shkl(muiso,SGData,cell):
3133	"this is to convert isotropic mustrain to generalized Shkls"
3134	import GSASIIlattice as G2lat
3135	A = G2lat.cell2AB(cell)[0]
3136
3137	def minMus(Shkl,muiso,H,SGData,A):
3138	U = np.inner(A.T,H)
3139	S = np.array(MustrainCoeff(U,SGData))
3140	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,nxs]),axis=0))
3141	rad = np.sqrt(np.sum((Sum[:,nxs]H)*2,axis=1))
3142	return (muiso-rad)**2
3143
3144	laue = SGData['SGLaue']
3145	PHI = np.linspace(0.,360.,60,True)
3146	PSI = np.linspace(0.,180.,60,True)
3147	X = np.outer(npsind(PHI),npsind(PSI))
3148	Y = np.outer(npcosd(PHI),npsind(PSI))
3149	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3150	HKL = np.dstack((X,Y,Z))
3151	if laue in ['m3','m3m']:
3152	S0 = [1000.,1000.]
3153	elif laue in ['6/m','6/mmm','3m1']:
3154	S0 = [1000.,1000.,1000.]
3155	elif laue in ['31m','3']:
3156	S0 = [1000.,1000.,1000.,1000.]
3157	elif laue in ['3R','3mR']:
3158	S0 = [1000.,1000.,1000.,1000.]
3159	elif laue in ['4/m','4/mmm']:
3160	S0 = [1000.,1000.,1000.,1000.]
3161	elif laue in ['mmm']:
3162	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3163	elif laue in ['2/m']:
3164	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3165	else:
3166	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
3167	1000.,1000.,0.,0.,0.]
3168	S0 = np.array(S0)
3169	HKL = np.reshape(HKL,(-1,3))
3170	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3171	return result[0]
3172
3173	def PackRot(SGOps):
3174	IRT = []
3175	for ops in SGOps:
3176	M = ops[0]
3177	irt = 0
3178	for j in range(2,-1,-1):
3179	for k in range(2,-1,-1):
3180	irt *= 3
3181	irt += M[k][j]
3182	IRT.append(int(irt))
3183	return IRT
3184
3185	def SytSym(XYZ,SGData):
3186	'''
3187	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3188
3189	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
3190	:param SGData: from SpcGroup
3191	:Returns: a four element tuple:
3192
3193	* The 1st element is a code for the site symmetry (see GetKNsym)
3194	* The 2nd element is the site multiplicity
3195	* Ndup number of overlapping operators
3196	* dupDir Dict - dictionary of overlapping operators
3197
3198	'''
3199	Mult = 1
3200	Isym = 0
3201	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3202	Isym = 1073741824
3203	Jdup = 0
3204	Ndup = 0
3205	dupDir = {}
3206	inv = SGData['SGInv']+1
3207	icen = SGData['SGCen']
3208	if SGData['SGFixed']: #already in list of operators
3209	inv = 1
3210	Xeqv = list(GenAtom(XYZ,SGData,True))
3211	# for xeqv in Xeqv: print(xeqv)
3212	IRT = PackRot(SGData['SGOps'])
3213	L = -1
3214	for ic,cen in enumerate(icen):
3215	for invers in range(int(inv)):
3216	for io,ops in enumerate(SGData['SGOps']):
3217	irtx = (1-2invers)IRT[io]
3218	L += 1
3219	if not Xeqv[L][1]:
3220	Ndup = io
3221	Jdup += 1
3222	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
3223	if jx < 39:
3224	px = GetOprName(str(irtx))
3225	if Xeqv[L][-1] < 0:
3226	if '(' in px:
3227	px = px.split('(')
3228	px[0] += "'"
3229	px = '('.join(px)
3230	else:
3231	px += "'"
3232	dupDir[px] = L
3233	Isym += 2**(jx-1)
3234	if Isym == 1073741824: Isym = 0
3235	Mult = len(SGData['SGOps'])len(SGData['SGCen'])inv//Jdup
3236
3237	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3238
3239	def MagSytSym(SytSym,dupDir,SGData):
3240	'''
3241	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3242	'''
3243	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3244	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3245	# print('dupDir',dupDir)
3246	SplitSytSym = SytSym.split('(')
3247	if SGData['SGGray']:
3248	return SytSym+"1'"
3249	if SytSym == '1': #genersl position
3250	return SytSym
3251	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
3252	try:
3253	MagSytSym = SGData['MagSpGrp'].split()[1]
3254	except IndexError:
3255	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3256	if len(SplitSytSym) > 1:
3257	MagSytSym += '('+SplitSytSym[1]
3258	return MagSytSym
3259	if len(dupDir) == 1:
3260	return list(dupDir.keys())[0]
3261
3262
3263	if '2/m' in SytSym: #done I think; last 2wo might be not needed
3264	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3265	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3266	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3267	'(110)':['2(110)','m(110)']}
3268
3269	elif '4/mmm' in SytSym:
3270	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
3271	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
3272	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
3273	elif '4mm' in SytSym:
3274	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3275	elif '422' in SytSym:
3276	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3277	elif '-4m2' in SytSym:
3278	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3279	elif '-42m' in SytSym:
3280	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3281	elif '-4' in SytSym:
3282	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3283	elif '4' in SytSym:
3284	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3285
3286	elif '222' in SytSym:
3287	ops = {'':['2(x)','2(y)','2(z)'],
3288	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3289	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3290	'(110)':['2(z)','2(110)','2(+-0)',],}
3291	elif 'mm2' in SytSym:
3292	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3293	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3294	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3295	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3296	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3297	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3298	'(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3299	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3300	'(110)':['m(z)','m(+-0)','2(110)',],}
3301	elif 'mmm' in SytSym:
3302	ops = {'':['m(x)','m(y)','m(z)'],
3303	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3304	'(110)':['m(z)','m(110)','m(+-0)',],
3305	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3306
3307	elif '32' in SytSym:
3308	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)']}
3309	elif '23' in SytSym:
3310	ops = {'':['2(x)','3(111)']}
3311	elif 'm3' in SytSym:
3312	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3313	elif '3m' in SytSym:
3314	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3315	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3316	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3317
3318	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m',]: #not simple cases
3319	MagSytSym = SytSym
3320	if "-1'" in dupDir:
3321	if '-6' in SytSym:
3322	MagSytSym = MagSytSym.replace('-6',"-6'")
3323	elif '-3m' in SytSym:
3324	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3325	elif '-3' in SytSym:
3326	MagSytSym = MagSytSym.replace('-3',"-3'")
3327	elif '-6m2' in SytSym:
3328	if "m'(110)" in dupDir:
3329	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3330	elif '6/m' in SytSym:
3331	if "m'(z)" in dupDir:
3332	MagSytSym = "6'/m'"
3333	elif '6mm' in SytSym:
3334	if "m'(110)" in dupDir:
3335	MagSytSym = "6'm'm"
3336	return MagSytSym
3337	try:
3338	axis = '('+SytSym.split('(')[1]
3339	except IndexError:
3340	axis = ''
3341	MagSytSym = ''
3342	for m in ops[axis]:
3343	if m in dupDir:
3344	MagSytSym += m.split('(')[0]
3345	else:
3346	MagSytSym += m.split('(')[0]+"'"
3347	if '2/m' in SytSym and '2' in m:
3348	MagSytSym += '/'
3349	if '-3/m' in SytSym:
3350	MagSytSym = '-'+MagSytSym
3351
3352	MagSytSym += axis
3353	# some exceptions & special rules
3354	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3355	return MagSytSym
3356
3357	# if len(GenSym) == 3:
3358	# if SGSpin[1] < 0:
3359	# if 'mm2' in SytSym:
3360	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3361	# else: #bad rule for I41/a
3362	# MagSytSym = SplitSytSym[0]+"'"
3363	# if len(SplitSytSym) > 1:
3364	# MagSytSym += '('+SplitSytSym[1]
3365	# else:
3366	# MagSytSym = SytSym
3367	# if len(SplitSytSym) >1:
3368	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3369	# MagSytSym = MagSytSym.replace('-4',"-4'")
3370	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3371	# MagSytSym = MagSytSym.replace('-6',"-6'")
3372	# return MagSytSym
3373	#
3374	return SytSym
3375
3376	def UpdateSytSym(Phase):
3377	''' Update site symmetry/site multiplicity after space group/BNS lattice change
3378	'''
3379	generalData = Phase['General']
3380	SGData = generalData['SGData']
3381	Atoms = Phase['Atoms']
3382	cx,ct,cs,cia = generalData['AtomPtrs']
3383	for atom in Atoms:
3384	XYZ = atom[cx:cx+3]
3385	sytsym,Mult = SytSym(XYZ,SGData)[:2]
3386	sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3387	atom[cs] = sytsym
3388	if generalData['Type'] == 'magnetic':
3389	magSytSym = MagSytSym(sytSym,dupDir,SGData)
3390	atom[cs] = magSytSym
3391	atom[cs+1] = Mult
3392	return
3393
3394	def ElemPosition(SGData):
3395	''' Under development.
3396	Object here is to return a list of symmetry element types and locations suitable
3397	for say drawing them.
3398	So far I have the element type... getting all possible locations without lookup may be impossible!
3399	'''
3400	Inv = SGData['SGInv']
3401	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3402	# get operators & expand if centrosymmetric
3403	Ops = SGData['SGOps']
3404	opM = np.array([op[0].T for op in Ops])
3405	opT = np.array([op[1] for op in Ops])
3406	if Inv:
3407	opM = np.concatenate((opM,-opM))
3408	opT = np.concatenate((opT,-opT))
3409	opMT = list(zip(opM,opT))
3410	for M,T in opMT[1:]: #skip I
3411	Dt = int(nl.det(M))
3412	Tr = int(np.trace(M))
3413	Dt = -(Dt-1)//2
3414	Es = eleSym[Tr][Dt]
3415	if Dt: #rotation-inversion
3416	I = np.eye(3)
3417	if Tr == 1: #mirrors/glides
3418	if np.any(T): #glide
3419	M2 = np.inner(M,M)
3420	MT = np.inner(M,T)+T
3421	print ('glide',Es,MT)
3422	print (M2)
3423	else: #mirror
3424	print ('mirror',Es,T)
3425	print (I-M)
3426	X = [-1,-1,-1]
3427	elif Tr == -3: # pure inversion
3428	X = np.inner(nl.inv(I-M),T)
3429	print ('inversion',Es,X)
3430	else: #other rotation-inversion
3431	M2 = np.inner(M,M)
3432	MT = np.inner(M,T)+T
3433	print ('rot-inv',Es,MT)
3434	print (M2)
3435	X = [-1,-1,-1]
3436	else: #rotations
3437	print ('rotation',Es)
3438	X = [-1,-1,-1]
3439	#SymElements.append([Es,X])
3440
3441	return #SymElements
3442
3443	def ApplyStringOps(A,SGData,X,Uij=[]):
3444	'Needs a doc string'
3445	SGOps = SGData['SGOps']
3446	SGCen = SGData['SGCen']
3447	Ax = A.split('+')
3448	Ax[0] = int(Ax[0])
3449	iC = 1
3450	if Ax[0] < 0:
3451	iC = -1
3452	Ax[0] = abs(Ax[0])
3453	nA = Ax[0]%100-1
3454	cA = Ax[0]//100
3455	Cen = SGCen[cA]
3456	M,T = SGOps[nA]
3457	if len(Ax)>1:
3458	cellA = Ax[1].split(',')
3459	cellA = np.array([int(a) for a in cellA])
3460	else:
3461	cellA = np.zeros(3)
3462	newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3463	if len(Uij):
3464	U = Uij2U(Uij)
3465	U = np.inner(M,np.inner(U,M).T)
3466	newUij = U2Uij(U)
3467	return [newX,newUij]
3468	else:
3469	return newX
3470
3471	def ApplyStringOpsMom(A,SGData,Mom):
3472	'Needs a doc string'
3473	SGOps = SGData['SGOps']
3474	Ax = A.split('+')
3475	Ax[0] = int(Ax[0])
3476	iAx = abs(Ax[0])
3477	nA = iAx%100-1
3478	nC = len(SGOps)*(iAx//100)
3479	NA = nA
3480	if Ax[0] < 0:
3481	NA += len(SGOps)
3482	M,T = SGOps[nA]
3483	if SGData['SGGray']: #no nonzero moments for gray groups!
3484	newMom = [0.,0.,0.]
3485	else:
3486	newMom = np.inner(Mom,M).Tnl.det(M)SGData['SpnFlp'][NA+nC]
3487	# print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom)
3488	return newMom
3489
3490	def StringOpsProd(A,B,SGData):
3491	"""
3492	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3493	where '-' indicates inversion, c(>0) is the cell centering operator,
3494	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3495	Should return resultant string - C. SGData - dictionary using entries:
3496
3497	* 'SGCen': cell centering vectors [0,0,0] at least
3498	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3499
3500	"""
3501	SGOps = SGData['SGOps']
3502	SGCen = SGData['SGCen']
3503	#1st split out the cell translation part & work on the operator parts
3504	Ax = A.split('+'); Bx = B.split('+')
3505	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3506	iC = 0
3507	if Ax[0]*Bx[0] < 0:
3508	iC = 1
3509	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3510	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3511	cA = Ax[0]//100; cB = Bx[0]//100
3512	Cen = (SGCen[cA]+SGCen[cB])%1.0
3513	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3514	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3515	Mc = np.inner(Ma,Mb.T)
3516	# print Ma,Mb,Mc
3517	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3518	# print Ta,Tb,Tc
3519	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3520	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3521	#now the cell translation part
3522	if len(Ax)>1:
3523	cellA = Ax[1].split(',')
3524	cellA = [int(a) for a in cellA]
3525	else:
3526	cellA = [0,0,0]
3527	if len(Bx)>1:
3528	cellB = Bx[1].split(',')
3529	cellB = [int(b) for b in cellB]
3530	else:
3531	cellB = [0,0,0]
3532	cellC = np.add(cellA,cellB)
3533	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3534	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3535	return C
3536
3537	def U2Uij(U):
3538	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3539	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3540
3541	def Uij2U(Uij):
3542	#returns the thermal motion tensor U from Uij as numpy array
3543	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3544
3545	def StandardizeSpcName(spcgroup):
3546	'''Accept a spacegroup name where spaces may have not been used
3547	in the names according to the GSAS convention (spaces between symmetry
3548	for each axis) and return the space group name as used in GSAS
3549	'''
3550	rspc = spcgroup.replace(' ','').upper()
3551	# deal with rhombohedral and hexagonal setting designations
3552	rhomb = ''
3553	if rspc[-1:] == 'R':
3554	rspc = rspc[:-1]
3555	rhomb = ' R'
3556	gray = ''
3557	if "1'" in rspc:
3558	gray = " 1'"
3559	rspc = rspc.replace("1'",'')
3560	rspc = rspc.replace("'",'')
3561	if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3562	rspc = rspc[:-1]
3563	if rspc[1:3] in ['M3','N3','A3','D3']: #fix cubic old style
3564	rspc.replace('3','-3')
3565	bns = -1
3566	try:
3567	bns = rspc.index('_')
3568	rspc = rspc.replace(rspc[bns:bns+2],'')
3569	except ValueError:
3570	pass
3571	# look for a match in the spacegroup lists
3572	for i in spglist.values():
3573	for spc in i:
3574	if rspc == spc.replace(' ','').upper():
3575	return spc+gray+rhomb
3576	# how about the post-2002 orthorhombic names?
3577	if rspc in sgequiv_2002_orthorhombic:
3578	return sgequiv_2002_orthorhombic[rspc]+gray
3579	else:
3580	# not found
3581	return ''
3582
3583	def SpaceGroupNumber(spcgroup):
3584	SGNo = -1
3585	SpcGp = StandardizeSpcName(spcgroup)
3586	if not SpcGp:
3587	return SGNo
3588	try:
3589	SGNo = spgbyNum.index(SpcGp)
3590	except ValueError:
3591	pass
3592	return SGNo
3593
3594
3595	spgbyNum = []
3596	'''Space groups indexed by number'''
3597	spgbyNum = [None,
3598	'P 1','P -1', #1-2
3599	'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3600	'C 2/m','P 2/c','P 21/c','C 2/c', #3-15
3601	'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3602	'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3603	'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3604	'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3605	'C m m 2','C m c 21','C c c 2',
3606	'A m m 2','A b m 2','A m a 2','A b a 2',
3607	'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3608	'P m m m','P n n n','P c c m','P b a n',
3609	'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3610	'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3611	'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3612	'F m m m', 'F d d d',
3613	'I m m m','I b a m','I b c a','I m m a', #16-74
3614	'P 4','P 41','P 42','P 43',
3615	'I 4','I 41',
3616	'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3617	'I 4/m','I 41/a',
3618	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3619	'P 42 21 2','P 43 2 2','P 43 21 2',
3620	'I 4 2 2','I 41 2 2',
3621	'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3622	'P 42 m c','P 42 b c',
3623	'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3624	'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3625	'P -4 c 2','P -4 b 2','P -4 n 2',
3626	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3627	'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',
3628	'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',
3629	'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',
3630	'P 42/n c m',
3631	'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3632	'P 3','P 31','P 32','R 3','P -3','R -3',
3633	'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3634	'R 3 2',
3635	'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3636	'R 3 m','R 3 c',
3637	'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3638	'R -3 m','R -3 c', #75-167
3639	'P 6','P 61',
3640	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3641	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3642	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3643	'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
3644	'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3645	'F m -3','F d -3','I m -3',
3646	'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3647	'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3648	'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3649	'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3650	'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3651	'I m -3 m','I a -3 d',] #195-230
3652	altSettingOrtho = {}
3653	''' A dictionary of alternate settings for orthorhombic unit cells
3654	'''
3655	altSettingOrtho = {
3656	'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'},
3657	'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'},
3658	'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'},
3659	'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'},
3660	'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'},
3661	'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'},
3662	'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'},
3663	'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'},
3664	'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'},
3665	'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'},
3666	'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'},
3667	'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'},
3668	'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'},
3669	'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'},
3670	'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'},
3671	'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'},
3672	'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'},
3673	'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'},
3674	'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'},
3675	'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'},
3676	'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'},
3677	'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'},
3678	'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'},
3679	'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'},
3680	'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'},
3681	'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'},
3682	'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'},
3683	'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'},
3684	'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'},
3685	'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'},
3686	'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'},
3687	'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'},
3688	'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'},
3689	'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'},
3690	'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'},
3691	'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'},
3692	'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'},
3693	'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'},
3694	'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'},
3695	'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'},
3696	'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'},
3697	'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'},
3698	'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'},
3699	'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'},
3700	'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'},
3701	'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'},
3702	'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'},
3703	'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'},
3704	'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'},
3705	'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'},
3706	'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'},
3707	'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'},
3708	'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'},
3709	'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'},
3710	'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'},
3711	'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'},
3712	'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'},
3713	'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'},
3714	'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'},
3715	}
3716	spg2origins = {}
3717	''' A dictionary of all spacegroups that have 2nd settings; the value is the
3718	1st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3719	'''
3720	spg2origins = {
3721	'P n n n':[-.25,-.25,-.25],
3722	'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3723	'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3724	'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3725	'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3726	'F d d d':[-.125,-.125,-.125],
3727	'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3728	'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],
3729	'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],
3730	'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3731	'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3732	'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3733	'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3734	'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3735	spglist = {}
3736	'''A dictionary of space groups as ordered and named in the pre-2002 International
3737	Tables Volume A, except that spaces are used following the GSAS convention to
3738	separate the different crystallographic directions.
3739	Note that the symmetry codes here will recognize many non-standard space group
3740	symbols with different settings. They are ordered by Laue group
3741	'''
3742	spglist = {
3743	'P1' : ('P 1','P -1',), # 1-2
3744	'C1' : ('C 1','C -1',),
3745	'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3746	'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
3747	'C2/m':('C 2','C m','C c','C n',
3748	'C 2/m','C 2/c','C 2/n',),
3749	'A2/m':('A 2','A m','A a','A n',
3750	'A 2/m','A 2/a','A 2/n',),
3751	'Pmmm':('P 2 2 2',
3752	'P 2 2 21','P 21 2 2','P 2 21 2',
3753	'P 21 21 2','P 2 21 21','P 21 2 21',
3754	'P 21 21 21',
3755	'P m m 2','P 2 m m','P m 2 m',
3756	'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3757	'P c c 2','P 2 a a','P b 2 b',
3758	'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3759	'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3760	'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3761	'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3762	'P b a 2','P 2 c b','P c 2 a',
3763	'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3764	'P n n 2','P 2 n n','P n 2 n',
3765	'P m m m','P n n n',
3766	'P c c m','P m a a','P b m b',
3767	'P b a n','P n c b','P c n a',
3768	'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3769	'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3770	'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3771	'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3772	'P b a m','P m c b','P c m a',
3773	'P c c n','P n a a','P b n b',
3774	'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3775	'P n n m','P m n n','P n m n',
3776	'P m m n','P n m m','P m n m',
3777	'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3778	'P b c a','P c a b',
3779	'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3780	),
3781	'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3782	'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3783	'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3784	'C m c m','C c m m','C m c a','C c m b',
3785	'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3786	'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3787	'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3788	'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3789	'A m m a','A m a m','A b m a','A c a m',
3790	'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3791	'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3792	'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3793	'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3794	'B b m m','B m m b','B b c m','B m a b',
3795	'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3796	'Immm':('I 2 2 2','I 21 21 21',
3797	'I m m 2','I m 2 m','I 2 m m',
3798	'I b a 2','I 2 c b','I c 2 a',
3799	'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3800	'I m m m','I b a m','I m c b','I c m a',
3801	'I b c a','I c a b',
3802	'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3803	'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3804	'F m m 2','F m 2 m','F 2 m m',
3805	'F d d 2','F d 2 d','F 2 d d',),
3806	'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3807	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3808	'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3809	'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3810	'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3811	'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',
3812	'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',
3813	'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',
3814	'P 42/n c m',),
3815	'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',
3816	'I 4 c m','I 41 m d','I 41 c d',
3817	'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',
3818	'I 41/a m d','I 41/a c d'),
3819	'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3820	'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3821	'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3822	'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3823	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3824	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3825	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3826	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3827	'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',
3828	'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3829	'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3830	'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3831	'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',
3832	'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'),
3833	'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3834	'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3835	'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3836	'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3837	}
3838	sgequiv_2002_orthorhombic = {}
3839	''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3840	along with the pre-2002 name. The e designates a double glide-plane
3841	'''
3842	sgequiv_2002_orthorhombic = {
3843	'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3844	'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3845	'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3846	'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3847	'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3848	'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3849	'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3850	'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3851
3852	#'A few non-standard space groups for test use'
3853	nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h',
3854	'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3855	'R 3 c r','R -3 c r','R -3 m r',),
3856
3857	#Use the space groups types in this order to list the symbols in the
3858	#order they are listed in the International Tables, vol. A'''
3859	symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal',
3860	'trigonal', 'hexagonal', 'cubic')
3861
3862	# self-test materials follow. Requires files in directory testinp
3863	selftestlist = []
3864	'''Defines a list of self-tests'''
3865	selftestquiet = True
3866	def _ReportTest():
3867	'Report name and doc string of current routine when ``selftestquiet`` is False'
3868	if not selftestquiet:
3869	import inspect
3870	caller = inspect.stack()[1][3]
3871	doc = eval(caller).__doc__
3872	if doc is not None:
3873	print('testing '+__file__+' with '+caller+' ('+doc+')')
3874	else:
3875	print('testing '+__file__()+" with "+caller)
3876	def test0():
3877	'''self-test #0: exercise MoveToUnitCell'''
3878	_ReportTest()
3879	msg = "MoveToUnitCell failed"
3880	assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3881	assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3882	assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3883	assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3884	selftestlist.append(test0)
3885
3886	def test1():
3887	'''self-test #1: SpcGroup against previous results'''
3888	#'''self-test #1: SpcGroup and SGPrint against previous results'''
3889	_ReportTest()
3890	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3891	if ospath.exists(testdir):
3892	if testdir not in sys.path: sys.path.insert(0,testdir)
3893	import spctestinp
3894	def CompareSpcGroup(spc, referr, refdict, reflist):
3895	'Compare output from GSASIIspc.SpcGroup with results from a previous run'
3896	# if an error is reported, the dictionary can be ignored
3897	msg0 = "CompareSpcGroup failed on space group %s" % spc
3898	result = SpcGroup(spc)
3899	if result[0] == referr and referr > 0: return True
3900	# #print result[1]['SpGrp']
3901	#msg = msg0 + " in list lengths"
3902	#assert len(keys) == len(refdict.keys()), msg
3903	for key in refdict.keys():
3904	if key == 'SGOps' or key == 'SGCen':
3905	msg = msg0 + (" in key %s length" % key)
3906	assert len(refdict[key]) == len(result[1][key]), msg
3907	for i in range(len(refdict[key])):
3908	msg = msg0 + (" in key %s level 0" % key)
3909	assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
3910	msg = msg0 + (" in key %s level 1" % key)
3911	assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
3912	else:
3913	msg = msg0 + (" in key %s" % key)
3914	assert result[1][key] == refdict[key], msg
3915	msg = msg0 + (" in key %s reflist" % key)
3916	#for (l1,l2) in zip(reflist, SGPrint(result[1])):
3917	# assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
3918	# for now disable SGPrint testing, output has changed
3919	#assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
3920	for spc in spctestinp.SGdat:
3921	CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
3922	selftestlist.append(test1)
3923
3924	def test2():
3925	'''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
3926	_ReportTest()
3927	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3928	if ospath.exists(testdir):
3929	if testdir not in sys.path: sys.path.insert(0,testdir)
3930	import sgtbxtestinp
3931	def CompareWcctbx(spcname, cctbx_in, debug=0):
3932	'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
3933	cctbx = cctbx_in[:] # make copy so we don't delete from the original
3934	spc = (SpcGroup(spcname))[1]
3935	if debug: print (spc['SpGrp'])
3936	if debug: print (spc['SGCen'])
3937	latticetype = spcname.strip().upper()[0]
3938	# lattice type of R implies Hexagonal centering", fix the rhombohedral settings
3939	if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
3940	assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
3941	onebar = [1]
3942	if spc['SGInv']: onebar.append(-1)
3943	for (op,off) in spc['SGOps']:
3944	for inv in onebar:
3945	for cen in spc['SGCen']:
3946	noff = off + cen
3947	noff = MoveToUnitCell(noff)[0]
3948	mult = tuple((op*inv).ravel().tolist())
3949	if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
3950	for refop in cctbx:
3951	if debug: print (refop)
3952	# check the transform
3953	if refop[:9] != mult: continue
3954	if debug: print ("mult match")
3955	# check the translation
3956	reftrans = list(refop[-3:])
3957	reftrans = MoveToUnitCell(reftrans)[0]
3958	if all(abs(noff - reftrans) < 1.e-5):
3959	cctbx.remove(refop)
3960	break
3961	else:
3962	assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
3963	for key in sgtbxtestinp.sgtbx:
3964	CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
3965	selftestlist.append(test2)
3966
3967	def test3():
3968	'''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
3969	for selected space groups against info in IT Volume A '''
3970	_ReportTest()
3971	def ExerciseSiteSym (spc, crdlist):
3972	'compare site symmetries and multiplicities for a specified space group'
3973	msg = "failed on site sym test for %s" % spc
3974	(E,S) = SpcGroup(spc)
3975	assert not E, msg
3976	for t in crdlist:
3977	symb, m, n, od = SytSym(t[0],S)
3978	if symb.strip() != t[2].strip() or m != t[1]:
3979	print (spc,t[0],m,n,symb,t[2],od)
3980	assert m == t[1]
3981	#assert symb.strip() == t[2].strip()
3982
3983	ExerciseSiteSym('p 1',[
3984	((0.13,0.22,0.31),1,'1'),
3985	((0,0,0),1,'1'),
3986	])
3987	ExerciseSiteSym('p -1',[
3988	((0.13,0.22,0.31),2,'1'),
3989	((0,0.5,0),1,'-1'),
3990	])
3991	ExerciseSiteSym('C 2/c',[
3992	((0.13,0.22,0.31),8,'1'),
3993	((0.0,.31,0.25),4,'2(y)'),
3994	((0.25,.25,0.5),4,'-1'),
3995	((0,0.5,0),4,'-1'),
3996	])
3997	ExerciseSiteSym('p 2 2 2',[
3998	((0.13,0.22,0.31),4,'1'),
3999	((0,0.5,.31),2,'2(z)'),
4000	((0.5,.31,0.5),2,'2(y)'),
4001	((.11,0,0),2,'2(x)'),
4002	((0,0.5,0),1,'222'),
4003	])
4004	ExerciseSiteSym('p 4/n',[
4005	((0.13,0.22,0.31),8,'1'),
4006	((0.25,0.75,.31),4,'2(z)'),
4007	((0.5,0.5,0.5),4,'-1'),
4008	((0,0.5,0),4,'-1'),
4009	((0.25,0.25,.31),2,'4(001)'),
4010	((0.25,.75,0.5),2,'-4(001)'),
4011	((0.25,.75,0.0),2,'-4(001)'),
4012	])
4013	ExerciseSiteSym('p 31 2 1',[
4014	((0.13,0.22,0.31),6,'1'),
4015	((0.13,0.0,0.833333333),3,'2(100)'),
4016	((0.13,0.13,0.),3,'2(110)'),
4017	])
4018	ExerciseSiteSym('R 3 c',[
4019	((0.13,0.22,0.31),18,'1'),
4020	((0.0,0.0,0.31),6,'3'),
4021	])
4022	ExerciseSiteSym('R 3 c R',[
4023	((0.13,0.22,0.31),6,'1'),
4024	((0.31,0.31,0.31),2,'3(111)'),
4025	])
4026	ExerciseSiteSym('P 63 m c',[
4027	((0.13,0.22,0.31),12,'1'),
4028	((0.11,0.22,0.31),6,'m(100)'),
4029	((0.333333,0.6666667,0.31),2,'3m(100)'),
4030	((0,0,0.31),2,'3m(100)'),
4031	])
4032	ExerciseSiteSym('I a -3',[
4033	((0.13,0.22,0.31),48,'1'),
4034	((0.11,0,0.25),24,'2(x)'),
4035	((0.11,0.11,0.11),16,'3(111)'),
4036	((0,0,0),8,'-3(111)'),
4037	])
4038	selftestlist.append(test3)
4039
4040	if __name__ == '__main__':
4041	# run self-tests
4042	selftestquiet = False
4043	for test in selftestlist:
4044	test()
4045	print ("OK")

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

Download in other formats: