Context Navigation

source: trunk/GSASIIspc.py @ 3247

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

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

Download in other formats: