Context Navigation

source: trunk/GSASIIspc.py @ 4475

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

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

Download in other formats: