Context Navigation

source: trunk/GSASIIspc.py @ 3286

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

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

Download in other formats: