Context Navigation

source: trunk/GSASIIspc.py @ 4860

Last change on this file since 4860 was 4810, checked in by vondreele, 4 years ago
fix depreciation warning in G2indx.findMV fix numpy round bug in G2lattice.MaxIndex? fix bug in G2pwdGUI.OnFindMV fix bug in G2spc.SSpcGroup
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 174.9 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: 2021-02-12 15:51:39 +0000 (Fri, 12 Feb 2021) $
12	# $Author: vondreele $
13	# $Revision: 4810 $
14	# $URL: trunk/GSASIIspc.py $
15	# $Id: GSASIIspc.py 4810 2021-02-12 15:51:39Z 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: 4810 $")
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	if mod:
1405	return [-SSGKl[i] if mod[i] in ['a','b','g'] else SSGKl[i] for i in range(3)]
1406	else:
1407	return [SSGKl[i] for i in range(3)]
1408
1409	def extendSSGOps(SSGOps):
1410	for OpA in SSGOps:
1411	OpAtxt = SSMT2text(OpA)
1412	if 't' not in OpAtxt:
1413	continue
1414	for OpB in SSGOps:
1415	OpBtxt = SSMT2text(OpB)
1416	if 't' not in OpBtxt:
1417	continue
1418	OpC = list(SGProd(OpB,OpA))
1419	OpC[1] %= 1.
1420	OpCtxt = SSMT2text(OpC)
1421	# print OpAtxt.replace(' ','')+' * '+OpBtxt.replace(' ','')+' = '+OpCtxt.replace(' ','')
1422	for k,OpD in enumerate(SSGOps):
1423	OpDtxt = SSMT2text(OpD)
1424	OpDtxt2 = ''
1425	if SGData['SGGray']:
1426	OpDtxt2 = SSMT2text([OpD[0],OpD[1]+np.array([0.,0.,0.,.5])])
1427	# print ' ('+OpCtxt.replace(' ','')+' = ? '+OpDtxt.replace(' ','')+')'
1428	if OpCtxt == OpDtxt:
1429	continue
1430	elif OpCtxt == OpDtxt2:
1431	continue
1432	elif OpCtxt.split(',')[:3] == OpDtxt.split(',')[:3]:
1433	if 't' not in OpDtxt:
1434	SSGOps[k] = OpC
1435	# print k,' new:',OpCtxt.replace(' ','')
1436	break
1437	else:
1438	OpCtxt = OpCtxt.replace(' ','')
1439	OpDtxt = OpDtxt.replace(' ','')
1440	Txt = OpCtxt+' conflicts with '+OpDtxt
1441	# print (Txt)
1442	return False,Txt
1443	return True,SSGOps
1444
1445	def findMod(modSym):
1446	for a in ['a','b','g']:
1447	if a in modSym:
1448	return a
1449
1450	def genSSGOps():
1451	SSGOps = SSGData['SSGOps'][:]
1452	iFrac = {}
1453	for i,frac in enumerate(SSGData['modSymb']):
1454	if frac in ['1/2','1/3','1/4','1/6','1']:
1455	iFrac[i] = frac+'.'
1456	# print SGData['SpGrp']+SSymbol
1457	# print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac,'modSymb',SSGData['modSymb']
1458	# set identity & 1,-1; triclinic
1459	SSGOps[0][0][3,3] = 1.
1460	## expand if centrosymmetric
1461	# if SGData['SGInv']:
1462	# SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
1463	# monoclinic - all done & all checked
1464	if SGData['SGPtGrp'] in ['2','m']: #OK
1465	SSGOps[1][0][3,3] = SSGKl[0]
1466	SSGOps[1][1][3] = genQ[0]
1467	for i in iFrac:
1468	SSGOps[1][0][3,i] = -SSGKl[0]
1469	elif SGData['SGPtGrp'] == '2/m': #OK
1470	SSGOps[1][0][3,3] = SSGKl[1]
1471	if 's' in gensym:
1472	SSGOps[1][1][3] = 0.5
1473	for i in iFrac:
1474	SSGOps[1][0][3,i] = SSGKl[0]
1475
1476	# orthorhombic - all OK not fully checked
1477	elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm']: #OK
1478	if SGData['SGPtGrp'] == '222':
1479	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[1,2],2:[1,3]},'b':{2:[3,2],0:[1,2]}} #OK
1480	elif SGData['SGPtGrp'] == 'mm2':
1481	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}} #OK
1482	elif SGData['SGPtGrp'] == 'm2m':
1483	OrOps = {'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]}} #OK
1484	elif SGData['SGPtGrp'] == '2mm':
1485	OrOps = {'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]},'g':{0:[1,3],1:[2,3]}} #OK
1486	a = findMod(SSGData['modSymb'])
1487	OrFrac = OrOps[a]
1488	for j in iFrac:
1489	for i in OrFrac[j]:
1490	SSGOps[i][0][3,j] = -2.eval(iFrac[j])SSGKl[i-1]
1491	for i in [0,1,2]:
1492	SSGOps[i+1][0][3,3] = SSGKl[i]
1493	SSGOps[i+1][1][3] = genQ[i]
1494	E,SSGOps = extendSSGOps(SSGOps)
1495	if not E:
1496	return E,SSGOps
1497	elif SGData['SGPtGrp'] == 'mmm': #OK
1498	OrOps = {'g':{0:[1,3],1:[2,3]},'a':{1:[2,1],2:[3,1]},'b':{0:[1,2],2:[3,2]}}
1499	a = findMod(SSGData['modSymb'])
1500	if a == 'g':
1501	SSkl = [1,1,1]
1502	elif a == 'a':
1503	SSkl = [-1,1,-1]
1504	else:
1505	SSkl = [1,-1,-1]
1506	OrFrac = OrOps[a]
1507	for j in iFrac:
1508	for i in OrFrac[j]:
1509	SSGOps[i][0][3,j] = -2.eval(iFrac[j])SSkl[i-1]
1510	for i in [0,1,2]:
1511	SSGOps[i+1][0][3,3] = SSkl[i]
1512	SSGOps[i+1][1][3] = genQ[i]
1513	E,SSGOps = extendSSGOps(SSGOps)
1514	if not E:
1515	return E,SSGOps
1516	# tetragonal - all done & checked
1517	elif SGData['SGPtGrp'] == '4': #OK
1518	SSGOps[1][0][3,3] = SSGKl[0]
1519	SSGOps[1][1][3] = genQ[0]
1520	if '1/2' in SSGData['modSymb']:
1521	SSGOps[1][0][3,1] = -1
1522	elif SGData['SGPtGrp'] == '-4': #OK
1523	SSGOps[1][0][3,3] = SSGKl[0]
1524	if '1/2' in SSGData['modSymb']:
1525	SSGOps[1][0][3,1] = 1
1526	elif SGData['SGPtGrp'] in ['4/m',]: #OK
1527	if '1/2' in SSGData['modSymb']:
1528	SSGOps[1][0][3,1] = -SSGKl[0]
1529	for i,j in enumerate([1,3]):
1530	SSGOps[j][0][3,3] = 1
1531	if genQ[i]:
1532	SSGOps[j][1][3] = genQ[i]
1533	E,SSGOps = extendSSGOps(SSGOps)
1534	if not E:
1535	return E,SSGOps
1536	elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]: #OK
1537	iGens = [1,4,5]
1538	if SGData['SGPtGrp'] in ['4mm','-4m2',]:
1539	iGens = [1,6,7]
1540	for i,j in enumerate(iGens):
1541	if '1/2' in SSGData['modSymb'] and i < 2:
1542	SSGOps[j][0][3,1] = SSGKl[i]
1543	SSGOps[j][0][3,3] = SSGKl[i]
1544	if genQ[i]:
1545	if 's' in gensym and j == 6:
1546	SSGOps[j][1][3] = -genQ[i]
1547	else:
1548	SSGOps[j][1][3] = genQ[i]
1549	E,SSGOps = extendSSGOps(SSGOps)
1550	if not E:
1551	return E,SSGOps
1552	elif SGData['SGPtGrp'] in ['4/mmm',]:#OK
1553	if '1/2' in SSGData['modSymb']:
1554	SSGOps[1][0][3,1] = -SSGKl[0]
1555	SSGOps[6][0][3,1] = SSGKl[1]
1556	if modsym:
1557	SSGOps[1][1][3] = -genQ[3]
1558	for i,j in enumerate([1,2,6,7]):
1559	SSGOps[j][0][3,3] = 1
1560	SSGOps[j][1][3] = genQ[i]
1561	E,Result = extendSSGOps(SSGOps)
1562	if not E:
1563	return E,Result
1564	else:
1565	SSGOps = Result
1566
1567	# trigonal - all done & checked
1568	elif SGData['SGPtGrp'] == '3': #OK
1569	SSGOps[1][0][3,3] = SSGKl[0]
1570	if '1/3' in SSGData['modSymb']:
1571	SSGOps[1][0][3,1] = -1
1572	SSGOps[1][1][3] = genQ[0]
1573	elif SGData['SGPtGrp'] == '-3': #OK
1574	SSGOps[1][0][3,3] = -SSGKl[0]
1575	if '1/3' in SSGData['modSymb']:
1576	SSGOps[1][0][3,1] = -1
1577	SSGOps[1][1][3] = genQ[0]
1578	elif SGData['SGPtGrp'] in ['312','3m','-3m','-3m1','3m1']: #OK
1579	if '1/3' in SSGData['modSymb']:
1580	SSGOps[1][0][3,1] = -1
1581	for i,j in enumerate([1,5]):
1582	if SGData['SGPtGrp'] in ['3m','-3m']:
1583	SSGOps[j][0][3,3] = 1
1584	else:
1585	SSGOps[j][0][3,3] = SSGKl[i+1]
1586	if genQ[i]:
1587	SSGOps[j][1][3] = genQ[i]
1588	elif SGData['SGPtGrp'] in ['321','32']: #OK
1589	for i,j in enumerate([1,4]):
1590	SSGOps[j][0][3,3] = SSGKl[i]
1591	if genQ[i]:
1592	SSGOps[j][1][3] = genQ[i]
1593	elif SGData['SGPtGrp'] in ['31m','-31m']: #OK
1594	ids = [1,3]
1595	if SGData['SGPtGrp'] == '-31m':
1596	ids = [1,3]
1597	if '1/3' in SSGData['modSymb']:
1598	SSGOps[ids[0]][0][3,1] = -SSGKl[0]
1599	for i,j in enumerate(ids):
1600	SSGOps[j][0][3,3] = 1
1601	if genQ[i+1]:
1602	SSGOps[j][1][3] = genQ[i+1]
1603
1604	# hexagonal all done & checked
1605	elif SGData['SGPtGrp'] == '6': #OK
1606	SSGOps[1][0][3,3] = SSGKl[0]
1607	SSGOps[1][1][3] = genQ[0]
1608	elif SGData['SGPtGrp'] == '-6': #OK
1609	SSGOps[1][0][3,3] = SSGKl[0]
1610	elif SGData['SGPtGrp'] in ['6/m',]: #OK
1611	SSGOps[1][0][3,3] = -SSGKl[1]
1612	SSGOps[1][1][3] = genQ[0]
1613	SSGOps[2][1][3] = genQ[1]
1614	elif SGData['SGPtGrp'] in ['622',]: #OK
1615	for i,j in enumerate([1,9,8]):
1616	SSGOps[j][0][3,3] = SSGKl[i]
1617	if genQ[i]:
1618	SSGOps[j][1][3] = -genQ[i]
1619	E,SSGOps = extendSSGOps(SSGOps)
1620
1621	elif SGData['SGPtGrp'] in ['6mm','-62m','-6m2',]: #OK
1622	for i,j in enumerate([1,6,7]):
1623	SSGOps[j][0][3,3] = SSGKl[i]
1624	if genQ[i]:
1625	SSGOps[j][1][3] = genQ[i]
1626	E,SSGOps = extendSSGOps(SSGOps)
1627	elif SGData['SGPtGrp'] in ['6/mmm',]: # OK
1628	for i,j in enumerate([1,2,10,11]):
1629	SSGOps[j][0][3,3] = 1
1630	if genQ[i]:
1631	SSGOps[j][1][3] = genQ[i]
1632	E,SSGOps = extendSSGOps(SSGOps)
1633	elif SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
1634	return True,SSGOps
1635	E,SSGOps = extendSSGOps(SSGOps)
1636	return E,SSGOps
1637
1638	def specialGen(gensym,modsym):
1639	sym = ''.join(gensym)
1640	if SGData['SGPtGrp'] in ['2/m',] and 'n' in SGData['SpGrp']:
1641	if 's' in sym:
1642	gensym = 'ss'
1643	if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
1644	gensym = '0ss'
1645	elif SGData['SGPtGrp'] in ['222',]:
1646	if sym == '00s':
1647	gensym = '0ss'
1648	elif sym == '0s0':
1649	gensym = 'ss0'
1650	elif sym == 's00':
1651	gensym = 's0s'
1652	elif SGData['SGPtGrp'] in ['mmm',]:
1653	if 'g' in modsym:
1654	if sym == 's00':
1655	gensym = 's0s'
1656	elif sym == '0s0':
1657	gensym = '0ss'
1658	elif 'a' in modsym:
1659	if sym == '0s0':
1660	gensym = 'ss0'
1661	elif sym == '00s':
1662	gensym = 's0s'
1663	elif 'b' in modsym:
1664	if sym == '00s':
1665	gensym = '0ss'
1666	elif sym == 's00':
1667	gensym = 'ss0'
1668	return gensym
1669
1670	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.}
1671	if SGData['SGLaue'] in ['m3','m3m']:
1672	return '(3+1) superlattices not defined for cubic space groups',None
1673	elif SGData['SGLaue'] in ['3R','3mR']:
1674	return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
1675	try:
1676	modsym,gensym = splitSSsym(SSymbol)
1677	except ValueError:
1678	return 'Error in superspace symbol '+SSymbol,None
1679	modQ = [Fracs[mod] for mod in modsym]
1680	SSGKl = SGData['SSGKl'][:]
1681	if SGData['SGLaue'] in ['2/m','mmm']:
1682	SSGKl = fixMonoOrtho()
1683	Ngen = len(gensym)
1684	if SGData.get('SGGray',False):
1685	Ngen -= 1
1686	if len(gensym) and Ngen != len(SSGKl):
1687	return 'Wrong number of items in generator symbol '+''.join(gensym),None
1688	gensym = specialGen(gensym[:Ngen],modsym)
1689	genQ = [Fracs[mod] for mod in gensym[:Ngen]]
1690	if not genQ:
1691	genQ = [0,0,0,0]
1692	SSgSpc = SGData['SpGrp']+SSymbol
1693	if SGData['SGGray']:
1694	SSgSpc = SSgSpc.replace('('," 1'(")
1695	SSGData = {'SSpGrp':SSgSpc,'modQ':modQ,'modSymb':modsym,'SSGKl':SSGKl}
1696	SSCen = np.zeros((len(SGData['SGCen']),4))
1697	for icen,cen in enumerate(SGData['SGCen']):
1698	SSCen[icen,0:3] = cen
1699	if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1700	Ncen = len(SGData['SGCen'])
1701	for icen in range(Ncen//2,Ncen):
1702	SSCen[icen,3] = 0.5
1703	SSGData['SSGCen'] = SSCen%1.
1704	SSGData['SSGOps'] = []
1705	for iop,op in enumerate(SGData['SGOps']):
1706	T = np.zeros(4)
1707	ssop = np.zeros((4,4))
1708	ssop[:3,:3] = op[0]
1709	T[:3] = op[1]
1710	SSGData['SSGOps'].append([ssop,T])
1711	E,Result = genSSGOps()
1712	if E:
1713	SSGData['SSGOps'] = Result
1714	if DEBUG:
1715	print ('Super spacegroup operators for '+SSGData['SSpGrp'])
1716	for Op in Result:
1717	print (SSMT2text(Op).replace(' ',''))
1718	if SGData['SGInv']:
1719	for Op in Result:
1720	Op = [-Op[0],-Op[1]%1.]
1721	print (SSMT2text(Op).replace(' ',''))
1722	return None,SSGData
1723	else:
1724	return Result+'\nOperator conflict - incorrect superspace symbol',None
1725
1726	def SSChoice(SGData):
1727	'''
1728	Gets the unique set of possible super space groups for a given space group
1729	'''
1730	ptgpSS = {'1':['(abg)',],'-1':['(abg)',],
1731
1732	'2':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1733	'm':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1734	'2/m':['(a0g)','(a1/2g)','(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)'],
1735
1736	'222':['(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	'mm2':['(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	'm2m':['(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	'2mm':['(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	'mmm':['(00g)','(1/20g)','(01/2g)','(1/21/2g)','(10g)','(01g)',
1749	'(a00)','(a1/20)','(a01/2)','(a1/21/2)','(a10)','(a01)',
1750	'(0b0)','(1/2b0)','(0b1/2)','(1/2b1/2)','(1b0)','(0b1)',],
1751
1752	'4':['(00g)','(1/21/2g)'],'4mm':['(00g)','(1/21/2g)'],
1753	'4/m':['(00g)','(1/21/2g)'],
1754	'422':['(00g)','(1/21/2g)'],'-4m2':['(00g)','(1/21/2g)'],'-42m':['(00g)','(1/21/2g)'],
1755	'4/mmm':['(00g)','(1/21/2g)'],
1756
1757	'3':['(00g)','(1/31/3g)'],'-3':['(00g)','(1/31/3g)'],
1758	'32':['(00g)'],'3m':['(00g)'],'-3m':['(00g)'],
1759	'321':['(00g)'],'3m1':['(00g)'],'-3m1':['(00g)'],
1760	'312':['(00g)','(1/31/3g)'],'31m':['(00g)','(1/31/3g)'],'-31m':['(00g)','(1/31/3g)'],
1761
1762	'6':['(00g)',],'6/m':['(00g)',],'-62m':['(00g)',],'-6m2':['(00g)',],
1763	'622':['(00g)',],'6/mmm':['(00g)',],'6mm':['(00g)',],
1764
1765	'23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1766
1767	ptgpTS = {'1':['0',],'-1':['0',],
1768
1769	'2':['0','s'],'m':['0','s'],
1770	'2/m':['00','0s','ss','s0'],
1771
1772	'222':['000','s00','0s0','00s',],
1773	'mm2':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1774	'm2m':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1775	'2mm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1776	'mmm':['000','s00','0s0','00s','ss0','s0s','0ss','q00','0q0','00q','0qq','q0q','qq0'],
1777
1778	'4':['0','q','s'],'4mm':['000','q00','s00','s0s','ss0','0ss','qq0','qqs'],
1779	'4/m':['00','s0'],'-4m2':['000','0s0','0q0'],'-42m':['000','00s'],
1780	'422':['000','q00','s00','s0s','ss0','0ss','qq0','qqs','0q0'],
1781	'4/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1782
1783	'3':['0','t'],'-3':['0','t'],
1784	'32':['00','t0'],'3m':['00','0s'],'-3m':['00','0s'],
1785	'321':['000','t00'],'3m1':['000','0s0'],'-3m1':['000','0s0'],
1786	'312':['000','t00'],'31m':['000','00s'],'-31m':['000','00s'],
1787
1788	'6':['0','h','t','s'],
1789	'6/m':['00','s0'],'-62m':['000','00s'],'-6m2':['000','0s0'],
1790	'622':['000','h00','t00','s00',],'6mm':['000','ss0','s0s','0ss',],
1791	'6/mmm':['0000','s0s0','00ss','s00s','ss00','0ss0','0s0s'],
1792
1793	'23':['',],'m3':['',],'432':['',],'-43m':['',],'m3m':['',]}
1794
1795	ptgp = SGData['SGPtGrp']
1796	SSChoice = []
1797	for ax in ptgpSS[ptgp]:
1798	for sx in ptgpTS[ptgp]:
1799	SSChoice.append(ax+sx)
1800	if SGData['SGGray']: SSChoice[-1] += 's'
1801	ssChoice = []
1802	ssHash = []
1803	for item in SSChoice:
1804	E,SSG = SSpcGroup(SGData,item)
1805	if SSG:
1806	sshash = hash(str(SSGPrint(SGData,SSG)[1]))
1807	if sshash not in ssHash:
1808	ssHash.append(sshash)
1809	ssChoice.append(item)
1810	return ssChoice
1811
1812	def splitSSsym(SSymbol):
1813	'''
1814	Splits supersymmetry symbol into two lists of strings
1815	'''
1816	mssym = SSymbol.replace(' ','').split(')')
1817	if len(mssym) > 1:
1818	modsym,gensym = mssym
1819	else:
1820	modsym = mssym[0]
1821	gensym = ''
1822	modsym = modsym.replace(',','')
1823	if "1'" in modsym:
1824	gensym = gensym[:-1]
1825	modsym = modsym.replace("1'",'')
1826	if gensym in ['0','00','000','0000']: #get rid of extraneous symbols
1827	gensym = ''
1828	nfrac = modsym.count('/')
1829	modsym = modsym.lstrip('(')
1830	if nfrac == 0:
1831	modsym = list(modsym)
1832	elif nfrac == 1:
1833	pos = modsym.find('/')
1834	if pos == 1:
1835	modsym = [modsym[:3],modsym[3],modsym[4]]
1836	elif pos == 2:
1837	modsym = [modsym[0],modsym[1:4],modsym[4]]
1838	else:
1839	modsym = [modsym[0],modsym[1],modsym[2:]]
1840	else:
1841	lpos = modsym.find('/')
1842	rpos = modsym.rfind('/')
1843	if lpos == 1 and rpos == 4:
1844	modsym = [modsym[:3],modsym[3:6],modsym[6]]
1845	elif lpos == 1 and rpos == 5:
1846	modsym = [modsym[:3],modsym[3],modsym[4:]]
1847	else:
1848	modsym = [modsym[0],modsym[1:4],modsym[4:]]
1849	gensym = list(gensym)
1850	return modsym,gensym
1851
1852	def SSGPrint(SGData,SSGData,AddInv=False):
1853	'''
1854	Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
1855
1856	:param SGData: space group data structure as defined in SpcGroup above.
1857	:param SSGData: from :func:`SSpcGroup`
1858	:returns:
1859	SSGText - list of strings with the superspace group details
1860	SGTable - list of strings for each of the operations
1861	'''
1862	nCen = len(SSGData['SSGCen'])
1863	Mult = nCenlen(SSGData['SSGOps'])(int(SGData['SGInv'])+1)
1864	if SGData.get('SGFixed',False):
1865	Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])
1866	SSsymb = SSGData['SSpGrp']
1867	if 'BNSlattsym' in SGData and '_' in SGData['BNSlattsym'][0]:
1868	SSsymb = SGData['BNSlattsym'][0]+SSsymb[1:]
1869	SSGCen = list(SSGData['SSGCen'])
1870	if SGData.get('SGGray',False):
1871	if SGData.get('SGFixed',False):
1872	Mult //= 2
1873	else:
1874	SSGCen += list(SSGData['SSGCen']+[0,0,0,0.5])
1875	SSGCen = np.array(SSGCen)%1.
1876	else:
1877	if "1'" in SSsymb: #leftover in nonmag phase in mcif file
1878	nCen //= 2
1879	Mult //= 2
1880	SSsymb = SSsymb.replace("1'",'')[:-1]
1881	SSGText = []
1882	SSGText.append(' Superspace Group: '+SSsymb)
1883	CentStr = 'centrosymmetric'
1884	if not SGData['SGInv']:
1885	CentStr = 'non'+CentStr
1886	if SGData['SGLatt'] in 'ABCIFR':
1887	SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
1888	else:
1889	SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())
1890	SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
1891	SGptGp = SGData['SGPtGrp']
1892	if SGData['SGGray']:
1893	SGptGp += "1'"
1894	SSGText.append(' The superlattice point group is '+SGptGp+', '+''.join([str(i) for i in SSGData['SSGKl']]))
1895	SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
1896	SSGText.append(' Multiplicity of a general site is '+str(Mult))
1897	if SGData['SGUniq'] in ['a','b','c']:
1898	SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
1899	if SGData['SGInv']:
1900	SSGText.append(' The inversion center is located at 0,0,0')
1901	if SGData['SGPolax']:
1902	SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
1903	SSGText.append(' ')
1904	if len(SSGCen) > 1:
1905	SSGText.append(' The equivalent positions are:')
1906	SSGText.append(' ('+SSLatt2text(SSGCen)+')+\n')
1907	else:
1908	SSGText.append(' The equivalent positions are:\n')
1909	SSGTable = []
1910	for i,Opr in enumerate(SSGData['SSGOps']):
1911	SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
1912	if AddInv and SGData['SGInv']:
1913	for i,Opr in enumerate(SSGData['SSGOps']):
1914	IOpr = [-Opr[0],-Opr[1]]
1915	SSGTable.append('(%2d) %s'%(i+1+len(SSGData['SSGOps']),SSMT2text(IOpr)))
1916	return SSGText,SSGTable
1917
1918	def SSGModCheck(Vec,modSymb,newMod=True):
1919	''' Checks modulation vector compatibility with supersymmetry space group symbol.
1920	if newMod: Superspace group symbol takes precidence & the vector will be modified accordingly
1921	'''
1922	Fracs = {'1/2':0.5,'1/3':1./3,'1':1.0,'0':0.,'a':0.,'b':0.,'g':0.}
1923	modQ = [Fracs[mod] for mod in modSymb]
1924	if newMod:
1925	newVec = Vec
1926	if not np.any(Vec):
1927	newVec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
1928	return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(newVec,modSymb,modQ)], \
1929	[True if mod in ['a','b','g'] else False for mod in modSymb]
1930	else:
1931	return Vec,[True if mod in ['a','b','g'] else False for mod in modSymb]
1932
1933	def SSMT2text(Opr):
1934	"From superspace group matrix/translation operator returns text version"
1935	XYZS = ('x','y','z','t') #Stokes, Campbell & van Smaalen notation
1936	TRA = (' ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
1937	Fld = ''
1938	M,T = Opr
1939	for j in range(4):
1940	IJ = ''
1941	for k in range(4):
1942	txt = str(int(round(M[j][k])))
1943	txt = txt.replace('1',XYZS[k]).replace('0','')
1944	if '2' in txt:
1945	txt += XYZS[k]
1946	if IJ and M[j][k] > 0:
1947	IJ += '+'+txt
1948	else:
1949	IJ += txt
1950	IK = int(round(T[j]*12))%12
1951	if IK:
1952	if not IJ:
1953	break
1954	if IJ[0] == '-':
1955	Fld += (TRA[IK]+IJ).rjust(8)
1956	else:
1957	Fld += (TRA[IK]+'+'+IJ).rjust(8)
1958	else:
1959	Fld += IJ.rjust(8)
1960	if j != 3: Fld += ', '
1961	return Fld
1962
1963	def SSLatt2text(SSGCen):
1964	"Lattice centering vectors to text"
1965	lattTxt = ''
1966	lattDir = {4:'1/3',6:'1/2',8:'2/3',0:'0'}
1967	for vec in SSGCen:
1968	lattTxt += ' '
1969	for item in vec:
1970	lattTxt += '%s,'%(lattDir[int(item*12)])
1971	lattTxt = lattTxt.rstrip(',')
1972	lattTxt += ';'
1973	lattTxt = lattTxt.rstrip(';').lstrip(' ')
1974	return lattTxt
1975
1976	def SSpaceGroup(SGSymbol,SSymbol):
1977	'''
1978	Print the output of SSpcGroup in a nicely formatted way.
1979
1980	:param SGSymbol: space group symbol with spaces between axial fields.
1981	:param SSymbol: superspace group symbol extension (string).
1982	:returns: nothing
1983	'''
1984
1985	E,A = SpcGroup(SGSymbol)
1986	if E > 0:
1987	print (SGErrors(E))
1988	return
1989	E,B = SSpcGroup(A,SSymbol)
1990	if E > 0:
1991	print (E)
1992	return
1993	for l in SSGPrint(B):
1994	print (l)
1995
1996	def SGProd(OpA,OpB):
1997	'''
1998	Form space group operator product. OpA & OpB are [M,V] pairs;
1999	both must be of same dimension (3 or 4). Returns [M,V] pair
2000	'''
2001	A,U = OpA
2002	B,V = OpB
2003	M = np.inner(B,A.T)
2004	W = np.inner(B,U)+V
2005	return M,W
2006
2007	def GetLittleGrpOps(SGData,vec):
2008	''' Find rotation part of operators that leave vec unchanged
2009
2010	:param SGData: space group data structure as defined in SpcGroup above.
2011	:param vec: a numpy array of fractional vector coordinates
2012	:returns: Little - list of operators [M,T] that form the little gropu
2013	'''
2014	Little = []
2015	Ops = SGData['SGOps'][:]
2016	if SGData['SGInv']:
2017	Ops += [[-M,-T] for [M,T] in Ops]
2018	for [M,T] in Ops:
2019	tvec = np.inner(M,vec)%1.
2020	if np.allclose(tvec,vec%1.):
2021	Little.append([M,T])
2022	return Little
2023
2024	def MoveToUnitCell(xyz):
2025	'''
2026	Translates a set of coordinates so that all values are >=0 and < 1
2027
2028	:param xyz: a list or numpy array of fractional coordinates
2029	:returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
2030	'''
2031	XYZ = np.array(xyz)%1.
2032	cell = np.asarray(np.rint(XYZ-xyz),dtype=np.int32)
2033	return XYZ,cell
2034
2035	def Opposite(XYZ,toler=0.0002):
2036	'''
2037	Gives opposite corner, edge or face of unit cell for position within tolerance.
2038	Result may be just outside the cell within tolerance
2039
2040	:param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
2041	:param toler: unit cell fraction tolerance making opposite
2042	:returns:
2043	XYZ: dict of opposite positions; key=unit cell & always contains XYZ
2044	'''
2045	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]]
2046	TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
2047	perm = TB*perm3
2048	cperm = ['%d,%d,%d'%(i,j,k) for i,j,k in perm]
2049	D = dict(zip(cperm,perm))
2050	new = {}
2051	for key in D:
2052	new[key] = np.array(D[key])+np.array(XYZ)
2053	return new
2054
2055	def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
2056	'''
2057	Generates the equivalent positions for a specified coordinate and space group
2058
2059	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
2060	:param SGData: from :func:`SpcGroup`
2061	:param All: True return all equivalent positions including duplicates;
2062	False return only unique positions
2063	:param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
2064	:param Move: True move generated atom positions to be inside cell
2065	False do not move atoms
2066	:return: [[XYZEquiv],Idup,[UijEquiv],spnflp]
2067
2068	* [XYZEquiv] is list of equivalent positions (XYZ is first entry)
2069	* Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
2070	* is centering operator number (1-4) and - is for inversion
2071	Cell = unit cell translations needed to put new positions inside cell
2072	[UijEquiv] - equivalent Uij; absent if no Uij given
2073	* +1/-1 for spin inversion of operator - empty if not magnetic
2074
2075	'''
2076	XYZEquiv = []
2077	UijEquiv = []
2078	Idup = []
2079	Cell = []
2080	inv = int(SGData['SGInv']+1)
2081	icen = SGData['SGCen']
2082	if SGData.get('SGFixed',False):
2083	inv = 1
2084	SpnFlp = SGData.get('SpnFlp',[])
2085	spnflp = []
2086	X = np.array(XYZ)
2087	mj = 0
2088	for ic,cen in enumerate(icen):
2089	C = np.array(cen)
2090	for invers in range(inv):
2091	for io,[M,T] in enumerate(SGData['SGOps']):
2092	idup = ((io+1)+100ic)(1-2*invers)
2093	XT = np.inner(M,X)+T
2094	if len(Uij):
2095	U = Uij2U(Uij)
2096	U = np.inner(M,np.inner(U,M).T)
2097	newUij = U2Uij(U)
2098	if invers:
2099	XT = -XT
2100	XT += C
2101	cell = np.zeros(3,dtype=np.int32)
2102	if Move:
2103	newX,cell = MoveToUnitCell(XT)
2104	else:
2105	newX = XT
2106	if All:
2107	if np.allclose(newX,X,atol=0.0002): #do we want %1. here?
2108	idup = False
2109	else:
2110	if True in [np.allclose(newX%1.,oldX%1.,atol=0.0002) for oldX in XYZEquiv]:
2111	idup = False
2112	if All or idup:
2113	XYZEquiv.append(newX)
2114	Idup.append(idup)
2115	Cell.append(cell)
2116	if len(Uij):
2117	UijEquiv.append(newUij)
2118	if len(SpnFlp):
2119	spnflp.append(SpnFlp[mj])
2120	else:
2121	spnflp.append(1)
2122	mj += 1
2123	if len(Uij):
2124	return zip(XYZEquiv,UijEquiv,Idup,Cell,spnflp)
2125	else:
2126	return zip(XYZEquiv,Idup,Cell,spnflp)
2127
2128	def GenHKL(HKL,SGData):
2129	''' Generates all equivlent reflections including Friedel pairs
2130	:param HKL: [h,k,l] must be integral values
2131	:param SGData: space group data obtained from SpcGroup
2132	:returns: array Uniq: equivalent reflections
2133	'''
2134
2135	Ops = SGData['SGOps']
2136	OpM = np.array([op[0] for op in Ops])
2137	Uniq = np.inner(OpM,HKL)
2138	Uniq = list(Uniq)+list(-1*Uniq)
2139	return np.array(Uniq)
2140
2141	def GenHKLf(HKL,SGData):
2142	'''
2143	Uses old GSAS Fortran routine genhkl.for
2144
2145	:param HKL: [h,k,l] must be integral values for genhkl.for to work
2146	:param SGData: space group data obtained from SpcGroup
2147	:returns: iabsnt,mulp,Uniq,phi
2148
2149	* iabsnt = True if reflection is forbidden by symmetry
2150	* mulp = reflection multiplicity including Friedel pairs
2151	* Uniq = numpy array of equivalent hkl in descending order of h,k,l
2152	* phi = phase offset for each equivalent h,k,l
2153
2154	'''
2155	hklf = list(HKL)+[0,] #could be numpy array!
2156	Ops = SGData['SGOps']
2157	OpM = np.array([op[0] for op in Ops],order='F')
2158	OpT = np.array([op[1] for op in Ops])
2159	Cen = np.array([cen for cen in SGData['SGCen']],order='F')
2160
2161	import pyspg
2162	Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
2163	h,k,l,f = Uniq
2164	Uniq=np.array(list(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])))
2165	phi = f[:Nuniq]
2166	return iabsnt,mulp,Uniq,phi
2167
2168	def checkSSLaue(HKL,SGData,SSGData):
2169	#Laue check here - Toss HKL if outside unique Laue part
2170	h,k,l,m = HKL
2171	if SGData['SGLaue'] == '2/m':
2172	if SGData['SGUniq'] == 'a':
2173	if 'a' in SSGData['modSymb'] and h == 0 and m < 0:
2174	return False
2175	elif 'b' in SSGData['modSymb'] and k == 0 and l ==0 and m < 0:
2176	return False
2177	else:
2178	return True
2179	elif SGData['SGUniq'] == 'b':
2180	if 'b' in SSGData['modSymb'] and k == 0 and m < 0:
2181	return False
2182	elif 'a' in SSGData['modSymb'] and h == 0 and l ==0 and m < 0:
2183	return False
2184	else:
2185	return True
2186	elif SGData['SGUniq'] == 'c':
2187	if 'g' in SSGData['modSymb'] and l == 0 and m < 0:
2188	return False
2189	elif 'a' in SSGData['modSymb'] and h == 0 and k ==0 and m < 0:
2190	return False
2191	else:
2192	return True
2193	elif SGData['SGLaue'] == 'mmm':
2194	if 'a' in SSGData['modSymb']:
2195	if h == 0 and m < 0:
2196	return False
2197	else:
2198	return True
2199	elif 'b' in SSGData['modSymb']:
2200	if k == 0 and m < 0:
2201	return False
2202	else:
2203	return True
2204	elif 'g' in SSGData['modSymb']:
2205	if l == 0 and m < 0:
2206	return False
2207	else:
2208	return True
2209	else: #tetragonal, trigonal, hexagonal (& triclinic?)
2210	if l == 0 and m < 0:
2211	return False
2212	else:
2213	return True
2214
2215	def checkHKLextc(HKL,SGData):
2216	'''
2217	Checks if reflection extinct - does not check centering
2218
2219	:param HKL: [h,k,l]
2220	:param SGData: space group data obtained from SpcGroup
2221	:returns: True if extinct; False if allowed
2222
2223	'''
2224	Ops = SGData['SGOps']
2225	OpM = np.array([op[0] for op in Ops])
2226	OpT = np.array([op[1] for op in Ops])
2227	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2228	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,3))
2229	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2230	for dhkl,phkl in zip(DHKL,PHKL)[1:]: #skip identity
2231	if dhkl.any():
2232	continue
2233	else:
2234	if phkl%1.:
2235	return True
2236	return False
2237
2238	def checkMagextc(HKL,SGData):
2239	'''
2240	Checks if reflection magnetically extinct; does fullcheck (centering, too)
2241	uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
2242
2243	:param HKL: [h,k,l]
2244	:param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
2245	:returns: True if magnetically extinct; False if allowed (to match GenHKLf)
2246
2247	'''
2248	Ops = SGData['SGOps']
2249	Ncen = len(SGData['SGCen'])
2250	OpM = np.array([op[0] for op in Ops])
2251	OpT = np.array([op[1] for op in Ops])
2252	if SGData['SGInv'] and not SGData['SGFixed']:
2253	OpM = np.vstack((OpM,-OpM))
2254	OpT = np.vstack((OpT,-OpT))%1.
2255	OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
2256	OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
2257	Spn = SGData['SpnFlp'][:len(OpM)]
2258	Mag = np.array([nl.det(opm) for opm in OpM])*Spn
2259	DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
2260	PHKL = np.reshape(np.cos(2.0np.pinp.inner(HKL,OpT))Mag,(-1,))[:,nxs,nxs]OpM #compute T(R,theta) eq(7)
2261	Ftest = np.random.rand(3) #random magnetic moment
2262	Psum = np.zeros(3)
2263	nsum = 0.
2264	nA = 0
2265	for dhkl,phkl in zip(DHKL,PHKL):
2266	if not np.allclose(dhkl,HKL): #test for eq(5)
2267	continue
2268	else:
2269	nA += 1
2270	nsum += np.trace(phkl) #eq(8)
2271	pterm = np.inner(Ftest,phkl) #eq(9)
2272	Psum += pterm
2273	if nsum/nA > 1.: #only need to look at nA=1 frok eq(8)
2274	return False
2275	if np.allclose(Psum,np.zeros(3)):
2276	return True
2277	else:
2278	if np.inner(HKL,Psum):
2279	return True
2280	return False
2281
2282	def checkSSextc(HKL,SSGData):
2283	Ops = SSGData['SSGOps']
2284	OpM = np.array([op[0] for op in Ops])
2285	OpT = np.array([op[1] for op in Ops])
2286	HKLS = np.array([HKL,-HKL]) #Freidel's Law
2287	DHKL = np.reshape(np.inner(HKLS,OpM)-HKL,(-1,4))
2288	PHKL = np.reshape(np.inner(HKLS,OpT),(-1,))
2289	for dhkl,phkl in list(zip(DHKL,PHKL))[1:]: #skip identity
2290	if dhkl.any():
2291	continue
2292	else:
2293	if phkl%1.:
2294	return False
2295	return True
2296
2297	################################################################################
2298	#### Site symmetry tables
2299	################################################################################
2300
2301	OprName = {
2302	'-6643': ['-1',1],'6479' : ['2(z)',2],'-6479': ['m(z)',3],
2303	'6481' : ['m(y)',4],'-6481': ['2(y)',5],'6641' : ['m(x)',6],
2304	'-6641': ['2(x)',7],'6591' : ['m(+-0)',8],'-6591': ['2(+-0)',9],
2305	'6531' : ['m(110)',10],'-6531': ['2(110)',11],'6537' : ['4(z)',12],
2306	'-6537': ['-4(z)',13],'975' : ['3(111)',14],'6456' : ['3',15],
2307	'-489' : ['3(+--)',16],'483' : ['3(-+-)',17],'-969' : ['3(--+)',18],
2308	'819' : ['m(+0-)',19],'-819' : ['2(+0-)',20],'2431' : ['m(0+-)',21],
2309	'-2431': ['2(0+-)',22],'-657' : ['m(xz)',23],'657' : ['2(xz)',24],
2310	'1943' : ['-4(x)',25],'-1943': ['4(x)',26],'-2429': ['m(yz)',27],
2311	'2429' : ['2(yz)',28],'639' : ['-4(y)',29],'-639' : ['4(y)',30],
2312	'-6484': ['2(010)',4],'6484' : ['m(010)',5],'-6668': ['2(100)',6],
2313	'6668' : ['m(100)',7],'-6454': ['2(120)',18],'6454' : ['m(120)',19],
2314	'-6638': ['2(210)',20],'6638' : ['m(210)',21], #search in SytSym ends at m(210)
2315	'2223' : ['3(+++)2',39],
2316	'6538' : ['6(z)1',40],'-2169':['3(--+)2',41],'2151' : ['3(+--)2',42],
2317	'2205' :['-3(-+-)2',43],'-2205':[' (-+-)2',44],'489' :['-3(+--)1',45],
2318	'801' : ['4(y)1',46],'1945' : ['4(x)3',47],'-6585': ['-4(z)3 ',48],
2319	'6585' : ['4(z)3',49],'6584' : ['3(z)2',50],'6666' : ['6(z)5 ',51],
2320	'6643' : ['1',52],'-801' : ['-4(y)1',53],'-1945': ['-4(x)3 ',54],
2321	'-6666': ['-6(z)5',55],'-6538': ['-6(z)1',56],'-2223':['-3(+++)2',57],
2322	'-975' :['-3(+++)1',58],'-6456': ['-3(z)1',59],'-483' :['-3(-+-)1',60],
2323	'969' :['-3(--+)1',61],'-6584': ['-3(z)2',62],'2169' :['-3(--+)2',63],
2324	'-2151':['-3(+--)2',64], }
2325
2326	KNsym = {
2327	'0' :' 1 ','1' :' -1 ','64' :' 2(x)','32' :' m(x)',
2328	'97' :' 2/m(x)','16' :' 2(y)','8' :' m(y)','25' :' 2/m(y)',
2329	'2' :' 2(z)','4' :' m(z)','7' :' 2/m(z)','134217728' :' 2(yz)',
2330	'67108864' :' m(yz)','201326593' :' 2/m(yz)','2097152' :' 2(0+-)','1048576' :' m(0+-)',
2331	'3145729' :'2/m(0+-)','8388608' :' 2(xz)','4194304' :' m(xz)','12582913' :' 2/m(xz)',
2332	'524288' :' 2(+0-)','262144' :' m(+0-)','796433' :'2/m(+0-)','1024' :' 2(xy)',
2333	'512' :' m(xy)','1537' :' 2/m(xy)','256' :' 2(+-0)','128' :' m(+-0)',
2334	'385' :'2/m(+-0)','76' :' mm2(x)','52' :' mm2(y)','42' :' mm2(z)',
2335	'135266336' :' mm2(yz)','69206048' :'mm2(0+-)','8650760' :' mm2(xz)','4718600' :'mm2(+0-)',
2336	'1156' :' mm2(xy)','772' :'mm2(+-0)','82' :' 222 ','136314944' :' 222(x)',
2337	'8912912' :' 222(y)','1282' :' 222(z)','127' :' mmm ','204472417' :' mmm(x)',
2338	'13369369' :' mmm(y)','1927' :' mmm(z)','33554496' :' 4(x)','16777280' :' -4(x)',
2339	'50331745' :'4/m(x)' ,'169869394' :'422(x)','84934738' :'-42m(x)','101711948' :'4mm(x)',
2340	'254804095' :'4/mmm(x)','536870928 ':' 4(y)','268435472' :' -4(y)','805306393' :'4/m(y)',
2341	'545783890' :'422(y)','272891986' :'-42m(y)','541327412' :'4mm(y)','818675839' :'4/mmm(y)',
2342	'2050' :' 4(z)','4098' :' -4(z)','6151' :'4/m(z)','3410' :'422(z)',
2343	'4818' :'-42m(z)','2730' :'4mm(z)','8191' :'4/mmm(z)','8192' :' 3(111)',
2344	'8193' :' -3(111)','2629888' :' 32(111)','1319040' :' 3m(111)','3940737' :'-3m(111)',
2345	'32768' :' 3(+--)','32769' :' -3(+--)','10519552' :' 32(+--)','5276160' :' 3m(+--)',
2346	'15762945' :'-3m(+--)','65536' :' 3(-+-)','65537' :' -3(-+-)','134808576' :' 32(-+-)',
2347	'67437056' :' 3m(-+-)','202180097' :'-3m(-+-)','131072' :' 3(--+)','131073' :' -3(--+)',
2348	'142737664' :' 32(--+)','71434368' :' 3m(--+)','214040961' :'-3m(--+)','237650' :' 23 ',
2349	'237695' :' m3 ','715894098' :' 432 ','358068946' :' -43m ','1073725439':' m3m ',
2350	'68157504' :' mm2(d100)','4456464' :' mm2(d010)','642' :' mm2(d001)','153092172' :'-4m2(x)',
2351	'277348404' :'-4m2(y)','5418' :'-4m2(z)','1075726335':' 6/mmm ','1074414420':'-6m2(100)',
2352	'1075070124':'-6m2(120)','1075069650':' 6mm ','1074414890':' 622 ','1073758215':' 6/m ',
2353	'1073758212':' -6 ','1073758210':' 6 ','1073759865':'-3m(100)','1075724673':'-3m(120)',
2354	'1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
2355	'1073758209':' -3 ','1073758208':' 3 ','1074135143':'mmm(100)','1075314719':'mmm(010)',
2356	'1073743751':'mmm(110)','1074004034':' mm2(z100)','1074790418':' mm2(z010)','1073742466':' mm2(z110)',
2357	'1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
2358	'1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
2359	'1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
2360	'1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
2361	'1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
2362	'1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
2363	'1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
2364	'1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':' -1 ',
2365	}
2366
2367	NXUPQsym = {
2368	'1' :(28,29,28,28),'-1' :( 1,29,28, 0),'2(x)' :(12,18,12,25),'m(x)' :(25,18,12,25),
2369	'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),
2370	'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),
2371	'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),
2372	'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),
2373	'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),
2374	'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),
2375	'mm2(x)' :(12,10, 0,-1),'mm2(y)' :(13,10, 0,-1),'mm2(z)' :(14,10, 0,-1),
2376	'mm2(yz)' :(10,13, 0,-1),'mm2(0+-)' :(11,13, 0,-1),'mm2(xz)' :( 8,12, 0,-1),'mm2(+0-)' :( 9,12, 0,-1),
2377	'mm2(xy)' :( 6,11, 0,-1),'222' :( 1,10, 0,-1),'222(x)' :( 1,13, 0,-1),
2378	'222(y)' :( 1,12, 0,-1),'222(z)' :( 1,11, 0,-1),'mmm' :( 1,10, 0,-1),'mmm(x)' :( 1,13, 0,-1),
2379	'mmm(y)' :( 1,12, 0,-1),'mmm(z)' :( 1,11, 0,-1),'4(x)' :(12, 4,12, 0),'-4(x)' :( 1, 4,12, 0),
2380	'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),
2381	'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),
2382	'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,),
2383	'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),
2384	'-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),
2385	'-3(111)' :( 1, 5, 2, 0),'32(111)' :( 1, 5, 0, 2),'3m(111)' :( 2, 5, 0, 2),'-3m(111)' :( 1, 5, 0,-1),
2386	'3(+--)' :( 5, 8, 5, 0),'-3(+--)' :( 1, 8, 5, 0),'32(+--)' :( 1, 8, 0, 5),'3m(+--)' :( 5, 8, 0, 5),
2387	'-3m(+--)' :( 1, 8, 0,-1),'3(-+-)' :( 4, 7, 4, 0),'-3(-+-)' :( 1, 7, 4, 0),'32(-+-)' :( 1, 7, 0, 4),
2388	'3m(-+-)' :( 4, 7, 0, 4),'-3m(-+-)' :( 1, 7, 0,-1),'3(--+)' :( 3, 6, 3, 0),'-3(--+)' :( 1, 6, 3, 0),
2389	'32(--+)' :( 1, 6, 0, 3),'3m(--+)' :( 3, 6, 0, 3),'-3m(--+)' :( 1, 6, 0,-1),'23' :( 1, 1, 0, 0),
2390	'm3' :( 1, 1, 0, 0),'432' :( 1, 1, 0, 0),'-43m' :( 1, 1, 0, 0),'m3m' :( 1, 1, 0, 0),
2391	'mm2(d100)':(12,13, 0,-1),'mm2(d010)':(13,12, 0,-1),'mm2(d001)':(14,11, 0,-1),'-4m2(x)' :( 1, 4, 0,-1),
2392	'-4m2(y)' :( 1, 3, 0,-1),'-4m2(z)' :( 1, 2, 0,-1),'6/mmm' :( 1, 9, 0,-1),'-6m2(100)':( 1, 9, 0,-1),
2393	'-6m2(120)':( 1, 9, 0,-1),'6mm' :(14, 9, 0,-1),'622' :( 1, 9, 0,-1),'6/m' :( 1, 9,14,-1),
2394	'-6' :( 1, 9,14, 0),'6' :(14, 9,14, 0),'-3m(100)' :( 1, 9, 0,-1),'-3m(120)' :( 1, 9, 0,-1),
2395	'3m(100)' :(14, 9, 0,14),'3m(120)' :(14, 9, 0,14),'32(100)' :( 1, 9, 0,14),'32(120)' :( 1, 9, 0,14),
2396	'-3' :( 1, 9,14, 0),'3' :(14, 9,14, 0),'mmm(100)' :( 1,14, 0,-1),'mmm(010)' :( 1,15, 0,-1),
2397	'mmm(110)' :( 1,11, 0,-1),'mm2(z100)':(14,14, 0,-1),'mm2(z010)':(14,15, 0,-1),'mm2(z110)':(14,11, 0,-1),
2398	'mm2(100)' :(12,14, 0,-1),'mm2(010)' :(13,15, 0,-1),'mm2(110)' :( 6,11, 0,-1),'mm2(120)' :(15,14, 0,-1),
2399	'mm2(210)' :(16,15, 0,-1),'mm2(+-0)' :( 7,11, 0,-1),'222(100)' :( 1,14, 0,-1),'222(010)' :( 1,15, 0,-1),
2400	'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),
2401	'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),
2402	'm(001)' :(23,16,14,23),'m(100)' :(26,25,12,26),'m(010)' :(27,28,13,27),'m(110)' :(18,19, 6,18),
2403	'm(120)' :(24,27,15,24),'m(210)' :(25,26,16,25),'m(+-0)' :(17,20, 7,17),'2(001)' :(14,16,14,23),
2404	'2(100)' :(12,25,12,26),'2(010)' :(13,28,13,27),'2(110)' :( 6,19, 6,18),'2(120)' :(15,27,15,24),
2405	'2(210)' :(16,26,16,25),'2(+-0)' :( 7,20, 7,17),'-1' :( 1,29,28, 0)
2406	}
2407
2408	CSxinel = [[], # 0th empty - indices are Fortran style
2409	[[0,0,0],[ 0.0, 0.0, 0.0]], #1 0 0 0
2410	[[1,1,1],[ 1.0, 1.0, 1.0]], #2 X X X
2411	[[1,1,1],[ 1.0, 1.0,-1.0]], #3 X X -X
2412	[[1,1,1],[ 1.0,-1.0, 1.0]], #4 X -X X
2413	[[1,1,1],[ 1.0,-1.0,-1.0]], #5 -X X X
2414	[[1,1,0],[ 1.0, 1.0, 0.0]], #6 X X 0
2415	[[1,1,0],[ 1.0,-1.0, 0.0]], #7 X -X 0
2416	[[1,0,1],[ 1.0, 0.0, 1.0]], #8 X 0 X
2417	[[1,0,1],[ 1.0, 0.0,-1.0]], #9 X 0 -X
2418	[[0,1,1],[ 0.0, 1.0, 1.0]], #10 0 Y Y
2419	[[0,1,1],[ 0.0, 1.0,-1.0]], #11 0 Y -Y
2420	[[1,0,0],[ 1.0, 0.0, 0.0]], #12 X 0 0
2421	[[0,1,0],[ 0.0, 1.0, 0.0]], #13 0 Y 0
2422	[[0,0,1],[ 0.0, 0.0, 1.0]], #14 0 0 Z
2423	[[1,1,0],[ 1.0, 2.0, 0.0]], #15 X 2X 0
2424	[[1,1,0],[ 2.0, 1.0, 0.0]], #16 2X X 0
2425	[[1,1,2],[ 1.0, 1.0, 1.0]], #17 X X Z
2426	[[1,1,2],[ 1.0,-1.0, 1.0]], #18 X -X Z
2427	[[1,2,1],[ 1.0, 1.0, 1.0]], #19 X Y X
2428	[[1,2,1],[ 1.0, 1.0,-1.0]], #20 X Y -X
2429	[[1,2,2],[ 1.0, 1.0, 1.0]], #21 X Y Y
2430	[[1,2,2],[ 1.0, 1.0,-1.0]], #22 X Y -Y
2431	[[1,2,0],[ 1.0, 1.0, 0.0]], #23 X Y 0
2432	[[1,0,2],[ 1.0, 0.0, 1.0]], #24 X 0 Z
2433	[[0,1,2],[ 0.0, 1.0, 1.0]], #25 0 Y Z
2434	[[1,1,2],[ 1.0, 2.0, 1.0]], #26 X 2X Z
2435	[[1,1,2],[ 2.0, 1.0, 1.0]], #27 2X X Z
2436	[[1,2,3],[ 1.0, 1.0, 1.0]], #28 X Y Z
2437	]
2438
2439	CSuinel = [[], # 0th empty - indices are Fortran style
2440	[[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
2441	[[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
2442	[[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
2443	[[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
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]], #5 A A A D D D
2445	[[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
2446	[[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
2447	[[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
2448	[[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
2449	[[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
2450	[[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
2451	[[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
2452	[[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
2453	[[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
2454	[[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
2455	[[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
2456	[[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
2457	[[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
2458	[[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
2459	[[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
2460	[[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
2461	[[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
2462	[[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
2463	[[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
2464	[[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
2465	[[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
2466	[[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
2467	[[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
2468	[[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
2469	]
2470
2471	################################################################################
2472	#### Site symmetry routines
2473	################################################################################
2474
2475	def GetOprPtrName(key):
2476	'Needs a doc string'
2477	try:
2478	oprName = OprName[key][0]
2479	except KeyError:
2480	return key
2481	return oprName.replace('(','').replace(')','')
2482
2483	def GetOprPtrNumber(key):
2484	'Needs a doc string'
2485	try:
2486	return OprName[key][1]
2487	except KeyError:
2488	return key
2489
2490	def GetOprName(key):
2491	'Needs a doc string'
2492	return OprName[key][0]
2493
2494	def GetKNsym(key):
2495	'Needs a doc string'
2496	try:
2497	return KNsym[key].strip()
2498	except KeyError:
2499	return 'sp'
2500
2501	def GetNXUPQsym(siteSym):
2502	'''
2503	The codes XUPQ are for lookup of symmetry constraints for position(X), thermal parm(U) & magnetic moments (P & Q)
2504	'''
2505	return NXUPQsym[siteSym]
2506
2507	def GetCSxinel(siteSym):
2508	"returns Xyz terms, multipliers, GUI flags"
2509	indx = GetNXUPQsym(siteSym.strip())
2510	return CSxinel[indx[0]]
2511
2512	def GetCSuinel(siteSym):
2513	"returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
2514	indx = GetNXUPQsym(siteSym.strip())
2515	return CSuinel[indx[1]]
2516
2517	def GetCSpqinel(SpnFlp,dupDir):
2518	"returns Mxyz terms, multipliers, GUI flags"
2519	CSI = [[1,2,3],[1.0,1.0,1.0]]
2520	for sopr in dupDir:
2521	# print (sopr,dupDir[sopr])
2522	opr = sopr.replace("'",'')
2523	indx = GetNXUPQsym(opr)
2524	if SpnFlp[dupDir[sopr]] > 0:
2525	csi = CSxinel[indx[2]] #P
2526	else:
2527	csi = CSxinel[indx[3]] #Q
2528	# print(opr,indx,csi,CSI)
2529	if not len(csi):
2530	return [[0,0,0],[0.,0.,0.]]
2531	for kcs in [0,1,2]:
2532	if csi[0][kcs] == 0 and CSI[0][kcs] != 0:
2533	jcs = CSI[0][kcs]
2534	for ics in [0,1,2]:
2535	if CSI[0][ics] == jcs:
2536	CSI[0][ics] = 0
2537	CSI[1][ics] = 0.
2538	elif CSI[0][ics] > jcs:
2539	CSI[0][ics] = CSI[0][ics]-1
2540	elif (CSI[0][kcs] == csi[0][kcs]) and (CSI[1][kcs] != csi[1][kcs]):
2541	CSI[1][kcs] = csi[1][kcs]
2542	elif CSI[0][kcs] >= csi[0][kcs]:
2543	CSI[0][kcs] = min(CSI[0][kcs],csi[0][kcs])
2544	if CSI[1][kcs] != csi[1][kcs]:
2545	if CSI[1][kcs] == 1.:
2546	CSI[1][kcs] = csi[1][kcs]
2547	# print(CSI)
2548	return CSI
2549
2550	def getTauT(tau,sop,ssop,XYZ,wave=np.zeros(3)):
2551	phase = np.sum(XYZ*wave)
2552	ssopinv = nl.inv(ssop[0])
2553	mst = ssopinv[3][:3]
2554	epsinv = ssopinv[3][3]
2555	sdet = nl.det(sop[0])
2556	ssdet = nl.det(ssop[0])
2557	dtau = mst(XYZ-sop[1])-epsinvssop[1][3]
2558	dT = 1.0
2559	if np.any(dtau%.5):
2560	sumdtau = np.sum(dtau%.5)
2561	dT = 0.
2562	if np.abs(sumdtau-.5) > 1.e-4:
2563	dT = np.tan(np.pi*sumdtau)
2564	tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3]+phase)
2565	return sdet,ssdet,dtau,dT,tauT
2566
2567	def OpsfromStringOps(A,SGData,SSGData):
2568	SGOps = SGData['SGOps']
2569	SSGOps = SSGData['SSGOps']
2570	Ax = A.split('+')
2571	Ax[0] = int(Ax[0])
2572	iC = 1
2573	if Ax[0] < 0:
2574	iC = -1
2575	iAx = abs(Ax[0])
2576	nA = iAx%100-1
2577	nC = iAx//100
2578	unit = [0,0,0]
2579	if len(Ax) > 1:
2580	unit = eval('['+Ax[1]+']')
2581	return SGOps[nA],SSGOps[nA],iC,SGData['SGCen'][nC],unit
2582
2583	def GetSSfxuinel(waveType,Stype,nH,XYZ,SGData,SSGData,debug=False):
2584
2585	def orderParms(CSI):
2586	parms = [0,]
2587	for csi in CSI:
2588	for i in [0,1,2]:
2589	if csi[i] not in parms:
2590	parms.append(csi[i])
2591	for csi in CSI:
2592	for i in [0,1,2]:
2593	csi[i] = parms.index(csi[i])
2594	return CSI
2595
2596	def fracCrenel(tau,Toff,Twid):
2597	Tau = (tau-Toff[:,nxs])%1.
2598	A = np.where(Tau<Twid[:,nxs],1.,0.)
2599	return A
2600
2601	def fracFourier(tau,nH,fsin,fcos):
2602	SA = np.sin(2.nHnp.pi*tau)
2603	CB = np.cos(2.nHnp.pi*tau)
2604	A = SA[nxs,nxs,:]*fsin[:,:,nxs]
2605	B = CB[nxs,nxs,:]*fcos[:,:,nxs]
2606	return A+B
2607
2608	def posFourier(tau,nH,psin,pcos):
2609	SA = np.sin(2nHnp.pi*tau)
2610	CB = np.cos(2nHnp.pi*tau)
2611	A = SA[nxs,nxs,:]*psin[:,:,nxs]
2612	B = CB[nxs,nxs,:]*pcos[:,:,nxs]
2613	return A+B
2614
2615	def posZigZag(tau,Tmm,XYZmax):
2616	DT = Tmm[1]-Tmm[0]
2617	slopeUp = 2.*XYZmax/DT
2618	slopeDn = 2.*XYZmax/(1.-DT)
2619	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])
2620	return A
2621
2622	def posBlock(tau,Tmm,XYZmax):
2623	A = np.array([np.where(Tmm[0] < t <= Tmm[1],XYZmax,-XYZmax) for t in tau])
2624	return A
2625
2626	def DoFrac():
2627	delt2 = np.eye(2)*0.001
2628	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2629	dFTP = []
2630	if siteSym == '1':
2631	CSI = [[1,0],[2,0]],2*[[1.,0.],]
2632	elif siteSym == '-1':
2633	CSI = [[1,0],[0,0]],2*[[1.,0.],]
2634	else:
2635	FSC = np.ones(2,dtype='i')
2636	CSI = [np.zeros((2),dtype='i'),np.zeros(2)]
2637	if 'Crenel' in waveType:
2638	dF = np.zeros_like(tau)
2639	else:
2640	dF = fracFourier(tau,nH,delt2[:1],delt2[1:]).squeeze()
2641	dFT = np.zeros_like(dF)
2642	dFTP = []
2643	for i in SdIndx:
2644	sop = Sop[i]
2645	ssop = SSop[i]
2646	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2647	fsc = np.ones(2,dtype='i')
2648	if 'Crenel' in waveType:
2649	dFT = np.zeros_like(tau)
2650	fsc = [1,1]
2651	else: #Fourier
2652	dFT = fracFourier(tauT,nH,delt2[:1],delt2[1:]).squeeze()
2653	dFT = nl.det(sop[0])*dFT
2654	dFT = dFT[:,np.argsort(tauT)]
2655	dFT[0] *= ssdet
2656	dFT[1] *= sdet
2657	dFTP.append(dFT)
2658
2659	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2660	fsc = [1,1]
2661	if dT:
2662	CSI = [[[1,0],[1,0]],[[1.,0.],[1/dT,0.]]]
2663	else:
2664	CSI = [[[1,0],[0,0]],[[1.,0.],[0.,0.]]]
2665	FSC = np.zeros(2,dtype='i')
2666	return CSI,dF,dFTP
2667	else:
2668	for i in range(2):
2669	if np.allclose(dF[i,:],dFT[i,:],atol=1.e-6):
2670	fsc[i] = 1
2671	else:
2672	fsc[i] = 0
2673	FSC &= fsc
2674	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,fsc)
2675	n = -1
2676	for i,F in enumerate(FSC):
2677	if F:
2678	n += 1
2679	CSI[0][i] = n+1
2680	CSI[1][i] = 1.0
2681
2682	return CSI,dF,dFTP
2683
2684	def DoXYZ():
2685	delt5 = np.ones(5)*0.001
2686	delt6 = np.eye(6)*0.001
2687	if 'Fourier' in waveType:
2688	dX = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2689	#3x6x12 modulated position array (X,Spos,tau)& force positive
2690	elif waveType in ['ZigZag','Block']:
2691	if waveType == 'ZigZag':
2692	dX = posZigZag(tau,delt5[:2],delt5[2:])
2693	else:
2694	dX = posBlock(tau,delt5[:2],delt5[2:])
2695	dXTP = []
2696	if siteSym == '1':
2697	CSI = [[1,0,0],[2,0,0],[3,0,0], [4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2698	elif siteSym == '-1':
2699	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.],]
2700	else:
2701	if 'Fourier' in waveType:
2702	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2703	elif waveType in ['ZigZag','Block']:
2704	CSI = [np.array([[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0]]),
2705	np.array([[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0],[1.0,.0,.0]])]
2706	XSC = np.ones(6,dtype='i')
2707	dXTP = []
2708	for i in SdIndx:
2709	sop = Sop[i]
2710	ssop = SSop[i]
2711	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2712	xsc = np.ones(6,dtype='i')
2713	if 'Fourier' in waveType:
2714	dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2715	elif waveType == 'ZigZag':
2716	dXT = posZigZag(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2717	elif waveType == 'Block':
2718	dXT = posBlock(tauT,delt5[:2],delt5[2:])+np.array(XYZ)[:,nxs,nxs]
2719	dXT = np.inner(sop[0],dXT.T) # X modulations array(3x6x49) -> array(3x49x6)
2720	dXT = np.swapaxes(dXT,1,2) # back to array(3x6x49)
2721	dXT[:,:3,:] = (ssdetsdet) # modify the sin component
2722	dXTP.append(dXT)
2723	if waveType == 'Fourier':
2724	for i in range(3):
2725	if not np.allclose(dX[i,i,:],dXT[i,i,:]):
2726	xsc[i] = 0
2727	if not np.allclose(dX[i,i+3,:],dXT[i,i+3,:]):
2728	xsc[i+3] = 0
2729	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2730	xsc[3:6] = 0
2731	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2732	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2733	if dT:
2734	if '(x)' in siteSym:
2735	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2736	if 'm' in siteSym and len(SdIndx) == 1:
2737	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.]
2738	elif '(y)' in siteSym:
2739	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2740	if 'm' in siteSym and len(SdIndx) == 1:
2741	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2742	elif '(z)' in siteSym:
2743	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2744	if 'm' in siteSym and len(SdIndx) == 1:
2745	CSI[1][3:] = [1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2746	else:
2747	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2748	if '4/mmm' in laue:
2749	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2750	if '(xy)' in siteSym:
2751	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2752	if dT:
2753	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2754	else:
2755	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2756	if '(xy)' in siteSym or '(+-0)' in siteSym:
2757	mul = 1
2758	if '(+-0)' in siteSym:
2759	mul = -1
2760	if np.allclose(dX[0,0,:],dXT[1,0,:]):
2761	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2762	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2763	xsc[3:5] = 0
2764	if np.allclose(dX[0,3,:],dXT[0,4,:]):
2765	CSI[0][:2] = [[12,0,0],[12,0,0]]
2766	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2767	xsc[:2] = 0
2768	else:
2769	for i in range(3):
2770	if not np.allclose(dX[:,i],dXT[i,:,i]):
2771	xsc[i] = 0
2772	XSC &= xsc
2773	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,xsc)
2774	if waveType == 'Fourier':
2775	n = -1
2776	if debug: print (XSC)
2777	for i,X in enumerate(XSC):
2778	if X:
2779	n += 1
2780	CSI[0][i][0] = n+1
2781	CSI[1][i][0] = 1.0
2782
2783	return list(CSI),dX,dXTP
2784
2785	def DoUij():
2786	delt12 = np.eye(12)*0.0001
2787	dU = posFourier(tau,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x12 array
2788	dUTP = []
2789	if siteSym == '1':
2790	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2791	[7,0,0],[8,0,0],[9,0,0],[10,0,0],[11,0,0],[12,0,0]],12*[[1.,0.,0.],]
2792	elif siteSym == '-1':
2793	CSI = 6*[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]], \
2794	6*[[0.,0.,0.],]+[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]
2795	else:
2796	CSI = [np.zeros((12,3),dtype='i'),np.zeros((12,3))]
2797	USC = np.ones(12,dtype='i')
2798	dUTP = []
2799	dtau = 0.
2800	for i in SdIndx:
2801	sop = Sop[i]
2802	ssop = SSop[i]
2803	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2804	usc = np.ones(12,dtype='i')
2805	dUT = posFourier(tauT,nH,delt12[:6],delt12[6:]) #Uij modulations - 6x12x49 array
2806	dUijT = np.rollaxis(np.rollaxis(np.array(Uij2U(dUT)),3),3) #convert dUT to 12x49x3x3
2807	dUijT = np.rollaxis(np.inner(np.inner(sop[0],dUijT),sop[0].T),3) #transform by sop - 3x3x12x49
2808	dUT = np.array(U2Uij(dUijT)) #convert to 6x12x49
2809	dUT = dUT[:,:,np.argsort(tauT)]
2810	dUT[:,:6,:] =(ssdetsdet)
2811	dUTP.append(dUT)
2812	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2813	if dT:
2814	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2815	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2816	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2817	[1./dT,0.,0.],[1./dT,0.,0.],[1./dT,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2818	else:
2819	CSI = [[[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0],
2820	[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],
2821	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.],
2822	[0.,0.,0.],[0.,0.,0.],[0.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2823	if 'mm2(x)' in siteSym and dT:
2824	CSI[1][9:] = [0.,0.,0.],[-dT,0.,0.],[0.,0.,0.]
2825	USC = [1,1,1,0,1,0,1,1,1,0,1,0]
2826	elif '(xy)' in siteSym and dT:
2827	CSI[0] = [[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0],
2828	[1,0,0],[1,0,0],[2,0,0],[3,0,0],[4,0,0],[4,0,0]]
2829	CSI[1][9:] = [[1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]]
2830	USC = [1,1,1,1,1,1,1,1,1,1,1,1]
2831	elif '(x)' in siteSym and dT:
2832	CSI[1][9:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2833	elif '(y)' in siteSym and dT:
2834	CSI[1][9:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2835	elif '(z)' in siteSym and dT:
2836	CSI[1][9:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2837	for i in range(6):
2838	if not USC[i]:
2839	CSI[0][i] = [0,0,0]
2840	CSI[1][i] = [0.,0.,0.]
2841	CSI[0][i+6] = [0,0,0]
2842	CSI[1][i+6] = [0.,0.,0.]
2843	else:
2844	for i in range(6):
2845	if not np.allclose(dU[i,i,:],dUT[i,i,:]): #sin part
2846	usc[i] = 0
2847	if not np.allclose(dU[i,i+6,:],dUT[i,i+6,:]): #cos part
2848	usc[i+6] = 0
2849	if np.any(dUT[1,0,:]):
2850	if '4/m' in siteSym:
2851	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2852	if ssop[1][3]:
2853	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2854	usc[9] = 1
2855	else:
2856	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2857	usc[9] = 0
2858	elif '4' in siteSym:
2859	CSI[0][6:8] = [[12,0,0],[12,0,0]]
2860	CSI[0][:2] = [[11,0,0],[11,0,0]]
2861	if ssop[1][3]:
2862	CSI[1][:2] = [[1.,0.,0.],[-1.,0.,0.]]
2863	CSI[1][6:8] = [[1.,0.,0.],[-1.,0.,0.]]
2864	usc[2] = 0
2865	usc[8] = 0
2866	usc[3] = 1
2867	usc[9] = 1
2868	else:
2869	CSI[1][:2] = [[1.,0.,0.],[1.,0.,0.]]
2870	CSI[1][6:8] = [[1.,0.,0.],[1.,0.,0.]]
2871	usc[2] = 1
2872	usc[8] = 1
2873	usc[3] = 0
2874	usc[9] = 0
2875	elif 'xy' in siteSym or '+-0' in siteSym:
2876	if np.allclose(dU[0,0,:],dUT[0,1,:]*sdet):
2877	CSI[0][4:6] = [[12,0,0],[12,0,0]]
2878	CSI[0][6:8] = [[11,0,0],[11,0,0]]
2879	CSI[1][4:6] = [[1.,0.,0.],[sdet,0.,0.]]
2880	CSI[1][6:8] = [[1.,0.,0.],[sdet,0.,0.]]
2881	usc[4:6] = 0
2882	usc[6:8] = 0
2883
2884	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,usc)
2885	USC &= usc
2886	if debug: print (USC)
2887	if not np.any(dtau%.5):
2888	n = -1
2889	for i,U in enumerate(USC):
2890	if U:
2891	n += 1
2892	CSI[0][i][0] = n+1
2893	CSI[1][i][0] = 1.0
2894
2895	return list(CSI),dU,dUTP
2896
2897	def DoMag():
2898	delt6 = np.eye(6)*0.001
2899	dM = posFourier(tau,nH,delt6[:3],delt6[3:]) #+np.array(Mxyz)[:,nxs,nxs]
2900	dMTP = []
2901	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2902	if siteSym == '1':
2903	CSI = [[1,0,0],[2,0,0],[3,0,0],[4,0,0],[5,0,0],[6,0,0]],6*[[1.,0.,0.],]
2904	elif siteSym in ['-1','mmm',]:
2905	CSI = 3[[0,0,0],]+[[1,0,0],[2,0,0],[3,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2906	elif siteSym in ['4(z)','422(z)']:
2907	CSI[0][0][0] = CSI[0][4][1] = 1
2908	CSI[1][0][0] = 1.0
2909	CSI[1][4][1] = -1.0
2910	elif siteSym in ['-4m2(z)','422(z)',]:
2911	CSI[0][5][0] = 1
2912	CSI[1][5][0] = 1.0
2913	elif siteSym in ['-32(100)','-3',]:
2914	CSI[0][2][0] = 1
2915	CSI[1][2][0] = 1.0
2916	elif siteSym in ['3',]:
2917	CSI[0][0][0] = CSI[0][3][0] = CSI[0][4][0] = 1
2918	CSI[1][0][0] = -np.sqrt(3.0)
2919	CSI[1][3][0] = 2.0
2920	CSI[1][4][0] = 1.0
2921	elif siteSym in ['622','2(100)','32(100)',]:
2922	CSI[0][0][0] = CSI[0][1][0] = CSI[0][3][0] = 1
2923	CSI[1][0][0] = 1.0
2924	CSI[1][1][0] = 2.0
2925	CSI[1][3][0] = np.sqrt(3.0)
2926	else:
2927	#3x6x12 modulated moment array (M,Spos,tau)& force positive
2928	CSI = [np.zeros((6,3),dtype='i'),np.zeros((6,3))]
2929	MSC = np.ones(6,dtype='i')
2930	dMTP = []
2931	for i in SdIndx:
2932	sop = Sop[i]
2933	ssop = SSop[i]
2934	sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
2935	msc = np.ones(6,dtype='i')
2936	dMT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:]) #+np.array(XYZ)[:,nxs,nxs]
2937	dMT = np.inner(sop[0],dMT.T) # X modulations array(3x6x49) -> array(3x49x6)
2938	dMT = np.swapaxes(dMT,1,2) # back to array(3x6x49)
2939	dMT[:,:3,:] = (ssdetsdet) # modify the sin component
2940	dMTP.append(dMT)
2941	for i in range(3):
2942	if not np.allclose(dM[i,i,:],sdet*dMT[i,i,:]):
2943	msc[i] = 0
2944	if not np.allclose(dM[i,i+3,:],sdet*dMT[i,i+3,:]):
2945	msc[i+3] = 0
2946	if np.any(dtau%.5) and ('1/2' in SSGData['modSymb'] or '1' in SSGData['modSymb']):
2947	msc[3:6] = 0
2948	CSI = [[[1,0,0],[2,0,0],[3,0,0], [1,0,0],[2,0,0],[3,0,0]],
2949	[[1.,0.,0.],[1.,0.,0.],[1.,0.,0.], [1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]]
2950	if dT:
2951	if '(x)' in siteSym:
2952	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2953	if 'm' in siteSym and len(SdIndx) == 1:
2954	CSI[1][3:] = [1./dT,0.,0.],[-dT,0.,0.],[-dT,0.,0.]
2955	elif '(y)' in siteSym:
2956	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2957	if 'm' in siteSym and len(SdIndx) == 1:
2958	CSI[1][3:] = [-dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]
2959	elif '(z)' in siteSym:
2960	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2961	if 'm' in siteSym and len(SdIndx) == 1:
2962	CSI[1][3:] = [-dT,0.,0.],[-dT,0.,0.],[1./dT,0.,0.]
2963	else:
2964	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2965	if '4/mmm' in laue:
2966	if siteSym in ['4/mmm(z)',]:
2967	CSI = 3[[0,0,0],]+[[0,0,0],[0,0,0],[1,0,0]],3[[0.,0.,0.],]+3*[[1.,0.,0.],]
2968	if np.any(dtau%.5) and '1/2' in SSGData['modSymb']:
2969	if '(xy)' in siteSym:
2970	CSI[0] = [[1,0,0],[1,0,0],[2,0,0], [1,0,0],[1,0,0],[2,0,0]]
2971	if dT:
2972	CSI[1][3:] = [[1./dT,0.,0.],[1./dT,0.,0.],[-dT,0.,0.]]
2973	else:
2974	CSI[1][3:] = [0.,0.,0.],[0.,0.,0.],[0.,0.,0.]
2975	if '(xy)' in siteSym or '(+-0)' in siteSym:
2976	mul = 1
2977	if '(+-0)' in siteSym:
2978	mul = -1
2979	if np.allclose(dM[0,0,:],dMT[1,0,:]):
2980	CSI[0][3:5] = [[11,0,0],[11,0,0]]
2981	CSI[1][3:5] = [[1.,0,0],[mul,0,0]]
2982	msc[3:5] = 0
2983	if np.allclose(dM[0,3,:],dMT[0,4,:]):
2984	CSI[0][:2] = [[12,0,0],[12,0,0]]
2985	CSI[1][:2] = [[1.,0,0],[mul,0,0]]
2986	msc[:2] = 0
2987	MSC &= msc
2988	if debug: print (SSMT2text(ssop).replace(' ',''),sdet,ssdet,epsinv,msc)
2989	n = -1
2990	if debug: print (MSC)
2991	for i,M in enumerate(MSC):
2992	if M:
2993	n += 1
2994	CSI[0][i][0] = n+1
2995	CSI[1][i][0] = 1.0
2996
2997	return list(CSI),dM,dMTP
2998
2999	if debug: print ('super space group: '+SSGData['SSpGrp'])
3000	xyz = np.array(XYZ)%1.
3001	SGOps = copy.deepcopy(SGData['SGOps'])
3002	laue = SGData['SGLaue']
3003	siteSym = SytSym(XYZ,SGData)[0].strip()
3004	if debug: print ('siteSym: '+siteSym)
3005	SSGOps = copy.deepcopy(SSGData['SSGOps'])
3006	#expand ops to include inversions if any
3007	if SGData['SGInv'] and not SGData['SGFixed']:
3008	for op,sop in zip(SGData['SGOps'],SSGData['SSGOps']):
3009	SGOps.append([-op[0],-op[1]%1.])
3010	SSGOps.append([-sop[0],-sop[1]%1.])
3011	#build set of sym ops around special position
3012	SSop = []
3013	Sop = []
3014	Sdtau = []
3015	for iop,Op in enumerate(SGOps):
3016	nxyz = (np.inner(Op[0],xyz)+Op[1])%1.
3017	if np.allclose(xyz,nxyz,1.e-4) and iop and MT2text(Op).replace(' ','') != '-X,-Y,-Z':
3018	SSop.append(SSGOps[iop])
3019	Sop.append(SGOps[iop])
3020	ssopinv = nl.inv(SSGOps[iop][0])
3021	mst = ssopinv[3][:3]
3022	epsinv = ssopinv[3][3]
3023	Sdtau.append(np.sum(mst(XYZ-SGOps[iop][1])-epsinvSSGOps[iop][1][3]))
3024	SdIndx = np.argsort(np.array(Sdtau)) # just to do in sensible order
3025	if debug: print ('special pos super operators: ',[SSMT2text(ss).replace(' ','') for ss in SSop])
3026	#setup displacement arrays
3027	tau = np.linspace(-1,1,49,True)
3028	#make modulation arrays - one parameter at a time
3029	if Stype == 'Sfrac':
3030	CSI,dF,dFTP = DoFrac()
3031	elif Stype == 'Spos':
3032	CSI,dF,dFTP = DoXYZ()
3033	CSI[0] = orderParms(CSI[0])
3034	elif Stype == 'Sadp':
3035	CSI,dF,dFTP = DoUij()
3036	CSI[0] = orderParms(CSI[0])
3037	elif Stype == 'Smag':
3038	CSI,dF,dFTP = DoMag()
3039
3040	if debug:
3041	return CSI,dF,dFTP
3042	else:
3043	return CSI,[],[]
3044
3045	def MustrainNames(SGData):
3046	'Needs a doc string'
3047	laue = SGData['SGLaue']
3048	uniq = SGData['SGUniq']
3049	if laue in ['m3','m3m']:
3050	return ['S400','S220']
3051	elif laue in ['6/m','6/mmm','3m1','3']:
3052	return ['S400','S004','S202']
3053	elif laue in ['31m',]:
3054	return ['S400','S004','S202','S301']
3055	elif laue in ['3R','3mR']:
3056	return ['S400','S220','S310','S211']
3057	elif laue in ['4/m','4/mmm']:
3058	return ['S400','S004','S220','S022']
3059	elif laue in ['mmm']:
3060	return ['S400','S040','S004','S220','S202','S022']
3061	elif laue in ['2/m']:
3062	SHKL = ['S400','S040','S004','S220','S202','S022']
3063	if uniq == 'a':
3064	SHKL += ['S013','S031','S211']
3065	elif uniq == 'b':
3066	SHKL += ['S301','S103','S121']
3067	elif uniq == 'c':
3068	SHKL += ['S130','S310','S112']
3069	return SHKL
3070	else:
3071	SHKL = ['S400','S040','S004','S220','S202','S022']
3072	SHKL += ['S310','S103','S031','S130','S301','S013']
3073	SHKL += ['S211','S121','S112']
3074	return SHKL
3075
3076	def HStrainVals(HSvals,SGData):
3077	laue = SGData['SGLaue']
3078	uniq = SGData['SGUniq']
3079	DIJ = np.zeros(6)
3080	if laue in ['m3','m3m']:
3081	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[0]]
3082	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3083	DIJ[:4] = [HSvals[0],HSvals[0],HSvals[1],HSvals[0]]
3084	elif laue in ['3R','3mR']:
3085	DIJ = [HSvals[0],HSvals[0],HSvals[0],HSvals[1],HSvals[1],HSvals[1]]
3086	elif laue in ['4/m','4/mmm']:
3087	DIJ[:3] = [HSvals[0],HSvals[0],HSvals[1]]
3088	elif laue in ['mmm']:
3089	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3090	elif laue in ['2/m']:
3091	DIJ[:3] = [HSvals[0],HSvals[1],HSvals[2]]
3092	if uniq == 'a':
3093	DIJ[5] = HSvals[3]
3094	elif uniq == 'b':
3095	DIJ[4] = HSvals[3]
3096	elif uniq == 'c':
3097	DIJ[3] = HSvals[3]
3098	else:
3099	DIJ = [HSvals[0],HSvals[1],HSvals[2],HSvals[3],HSvals[4],HSvals[5]]
3100	return DIJ
3101
3102	def HStrainNames(SGData):
3103	'Needs a doc string'
3104	laue = SGData['SGLaue']
3105	uniq = SGData['SGUniq']
3106	if laue in ['m3','m3m']:
3107	return ['D11','eA'] #add cubic strain term
3108	elif laue in ['6/m','6/mmm','3m1','31m','3']:
3109	return ['D11','D33']
3110	elif laue in ['3R','3mR']:
3111	return ['D11','D12']
3112	elif laue in ['4/m','4/mmm']:
3113	return ['D11','D33']
3114	elif laue in ['mmm']:
3115	return ['D11','D22','D33']
3116	elif laue in ['2/m']:
3117	Dij = ['D11','D22','D33']
3118	if uniq == 'a':
3119	Dij += ['D23']
3120	elif uniq == 'b':
3121	Dij += ['D13']
3122	elif uniq == 'c':
3123	Dij += ['D12']
3124	return Dij
3125	else:
3126	Dij = ['D11','D22','D33','D12','D13','D23']
3127	return Dij
3128
3129	def MustrainCoeff(HKL,SGData):
3130	'Needs a doc string'
3131	#NB: order of terms is the same as returned by MustrainNames
3132	laue = SGData['SGLaue']
3133	uniq = SGData['SGUniq']
3134	h,k,l = HKL
3135	Strm = []
3136	if laue in ['m3','m3m']:
3137	Strm.append(h4+k4+l**4)
3138	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3139	elif laue in ['6/m','6/mmm','3m1','3']:
3140	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3141	Strm.append(l**4)
3142	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3143	elif laue in ['31m',]:
3144	Strm.append(h4+k4+2.0kh*3+2.0hk3+3.0(hk)*2)
3145	Strm.append(l**4)
3146	Strm.append(3.0((hl)*2+(kl)*2+hkl*2))
3147	Strm.append(4.0lh**3)
3148	elif laue in ['3R','3mR']:
3149	Strm.append(h4+k4+l**4)
3150	Strm.append(3.0((hk)*2+(hl)*2+(kl)**2))
3151	Strm.append(2.0(hl*3+lk*3+kh*3)+2.0(lh3+kl*3+lk**3))
3152	Strm.append(4.0(klh2+hlk2+hkl*2))
3153	elif laue in ['4/m','4/mmm']:
3154	Strm.append(h4+k4)
3155	Strm.append(l**4)
3156	Strm.append(3.0(hk)**2)
3157	Strm.append(3.0((hl)*2+(kl)**2))
3158	elif laue in ['mmm']:
3159	Strm.append(h**4)
3160	Strm.append(k**4)
3161	Strm.append(l**4)
3162	Strm.append(3.0(hk)**2)
3163	Strm.append(3.0(hl)**2)
3164	Strm.append(3.0(kl)**2)
3165	elif laue in ['2/m']:
3166	Strm.append(h**4)
3167	Strm.append(k**4)
3168	Strm.append(l**4)
3169	Strm.append(3.0(hk)**2)
3170	Strm.append(3.0(hl)**2)
3171	Strm.append(3.0(kl)**2)
3172	if uniq == 'a':
3173	Strm.append(2.0kl**3)
3174	Strm.append(2.0lk**3)
3175	Strm.append(4.0klh*2)
3176	elif uniq == 'b':
3177	Strm.append(2.0lh**3)
3178	Strm.append(2.0hl**3)
3179	Strm.append(4.0hlk*2)
3180	elif uniq == 'c':
3181	Strm.append(2.0hk**3)
3182	Strm.append(2.0kh**3)
3183	Strm.append(4.0hkl*2)
3184	else:
3185	Strm.append(h**4)
3186	Strm.append(k**4)
3187	Strm.append(l**4)
3188	Strm.append(3.0(hk)**2)
3189	Strm.append(3.0(hl)**2)
3190	Strm.append(3.0(kl)**2)
3191	Strm.append(2.0kh**3)
3192	Strm.append(2.0hl**3)
3193	Strm.append(2.0lk**3)
3194	Strm.append(2.0hk**3)
3195	Strm.append(2.0lh**3)
3196	Strm.append(2.0kl**3)
3197	Strm.append(4.0klh*2)
3198	Strm.append(4.0hlk*2)
3199	Strm.append(4.0khl*2)
3200	return Strm
3201
3202	def MuShklMean(SGData,Amat,Shkl):
3203
3204	def genMustrain(xyz,Shkl):
3205	uvw = np.inner(Amat.T,xyz)
3206	Strm = np.array(MustrainCoeff(uvw,SGData))
3207	Sum = np.sum(np.multiply(Shkl,Strm))
3208	Sum = np.where(Sum > 0.01,Sum,0.01)
3209	Sum = np.sqrt(Sum)
3210	return Sum*xyz
3211
3212	PHI = np.linspace(0.,360.,30,True)
3213	PSI = np.linspace(0.,180.,30,True)
3214	X = np.outer(npcosd(PHI),npsind(PSI))
3215	Y = np.outer(npsind(PHI),npsind(PSI))
3216	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3217	XYZ = np.dstack((X,Y,Z))
3218	XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,Shkl))
3219	return np.sqrt(np.sum(XYZ**2)/900.)
3220
3221	def Muiso2Shkl(muiso,SGData,cell):
3222	"this is to convert isotropic mustrain to generalized Shkls"
3223	import GSASIIlattice as G2lat
3224	A = G2lat.cell2AB(cell)[0]
3225
3226	def minMus(Shkl,muiso,H,SGData,A):
3227	U = np.inner(A.T,H)
3228	S = np.array(MustrainCoeff(U,SGData))
3229	nS = S.shape[0]
3230	Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:nS,nxs]),axis=0))
3231	rad = np.sqrt(np.sum((Sum[:,nxs]H)*2,axis=1))
3232	return (muiso-rad)**2
3233
3234	laue = SGData['SGLaue']
3235	PHI = np.linspace(0.,360.,60,True)
3236	PSI = np.linspace(0.,180.,60,True)
3237	X = np.outer(npsind(PHI),npsind(PSI))
3238	Y = np.outer(npcosd(PHI),npsind(PSI))
3239	Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
3240	HKL = np.dstack((X,Y,Z))
3241	if laue in ['m3','m3m']:
3242	S0 = [1000.,1000.]
3243	elif laue in ['6/m','6/mmm']:
3244	S0 = [1000.,1000.,1000.]
3245	elif laue in ['31m','3','3m1']:
3246	S0 = [1000.,1000.,1000.]
3247	elif laue in ['3R','3mR']:
3248	S0 = [1000.,1000.,1000.,1000.]
3249	elif laue in ['4/m','4/mmm']:
3250	S0 = [1000.,1000.,1000.,1000.]
3251	elif laue in ['mmm']:
3252	S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
3253	elif laue in ['2/m']:
3254	S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
3255	else:
3256	S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000.,
3257	1000.,1000.,0.,0.,0.]
3258	S0 = np.array(S0)
3259	HKL = np.reshape(HKL,(-1,3))
3260	result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
3261	return result[0]
3262
3263	def PackRot(SGOps):
3264	IRT = []
3265	for ops in SGOps:
3266	M = ops[0]
3267	irt = 0
3268	for j in range(2,-1,-1):
3269	for k in range(2,-1,-1):
3270	irt *= 3
3271	irt += M[k][j]
3272	IRT.append(int(irt))
3273	return IRT
3274
3275	def SytSym(XYZ,SGData):
3276	'''
3277	Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
3278
3279	:param XYZ: an array, tuple or list containing 3 elements: x, y & z
3280	:param SGData: from SpcGroup
3281	:Returns: a four element tuple:
3282
3283	* The 1st element is a code for the site symmetry (see GetKNsym)
3284	* The 2nd element is the site multiplicity
3285	* Ndup number of overlapping operators
3286	* dupDir Dict - dictionary of overlapping operators
3287
3288	'''
3289	if SGData['SpGrp'] == 'P 1':
3290	return '1',1,1,{}
3291	Mult = 1
3292	Isym = 0
3293	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3294	Isym = 1073741824
3295	Jdup = 0
3296	Ndup = 0
3297	dupDir = {}
3298	inv = SGData['SGInv']+1
3299	icen = SGData['SGCen']
3300	Ncen = len(icen)
3301	if SGData['SGFixed']: #already in list of operators
3302	inv = 1
3303	if SGData['SGGray'] and Ncen > 1: Ncen //= 2
3304	Xeqv = list(GenAtom(np.array(XYZ)%1.,SGData,True))
3305	# for xeqv in Xeqv: print(xeqv)
3306	IRT = PackRot(SGData['SGOps'])
3307	L = -1
3308	for ic,cen in enumerate(icen[:Ncen]):
3309	for invers in range(int(inv)):
3310	for io,ops in enumerate(SGData['SGOps']):
3311	irtx = (1-2invers)IRT[io]
3312	L += 1
3313	if not Xeqv[L][1]:
3314	Ndup = io
3315	Jdup += 1
3316	jx = GetOprPtrNumber(str(irtx)) #[KN table no,op name,KNsym ptr]
3317	if jx < 39:
3318	px = GetOprName(str(irtx))
3319	if Xeqv[L][-1] < 0:
3320	if '(' in px:
3321	px = px.split('(')
3322	px[0] += "'"
3323	px = '('.join(px)
3324	else:
3325	px += "'"
3326	dupDir[px] = L
3327	Isym += 2**(jx-1)
3328	if Isym == 1073741824: Isym = 0
3329	try:
3330	Mult = len(SGData['SGOps'])Nceninv//Jdup
3331	except: # patch because Jdup is not getting incremented for most atoms!
3332	Mult = 0
3333
3334	return GetKNsym(str(Isym)),Mult,Ndup,dupDir
3335
3336	def MagSytSym(SytSym,dupDir,SGData):
3337	'''
3338	site sym operations: 1,-1,2,3,-3,4,-4,6,-6,m need to be marked if spin inversion
3339	'''
3340	SGData['GenSym'],SGData['GenFlg'] = GetGenSym(SGData)[:2]
3341	# print('SGPtGrp',SGData['SGPtGrp'],'SytSym',SytSym,'MagSpGrp',SGData['MagSpGrp'])
3342	# print('dupDir',dupDir)
3343	SplitSytSym = SytSym.split('(')
3344	if SGData['SGGray']:
3345	return SytSym+"1'"
3346	if SytSym in ['1','-1','2','3','6','m']: #axial position
3347	return SytSym
3348	if SplitSytSym[0] == SGData['SGPtGrp']: #simple cases
3349	try:
3350	MagSytSym = SGData['MagSpGrp'].split()[1]
3351	except IndexError:
3352	MagSytSym = SGData['MagSpGrp'][1:].strip("1'")
3353	if len(SplitSytSym) > 1:
3354	MagSytSym += '('+SplitSytSym[1]
3355	return MagSytSym
3356	if len(dupDir) == 1:
3357	return list(dupDir.keys())[0]
3358
3359
3360	if '2/m' in SytSym: #done I think; last 2wo might be not needed
3361	ops = {'(x)':['2(x)','m(x)'],'(y)':['2(y)','m(y)'],'(z)':['2(z)','m(z)'],
3362	'(100)':['2(100)','m(100)'],'(010)':['2(010)','m(010)'],'(001)':['2(001)','m(001)'],
3363	'(120)':['2(120)','m(120)'],'(210)':['2(210)','m(210)'],'(+-0)':['2(+-0)','m(+-0)'],
3364	'(110)':['2(110)','m(110)']}
3365
3366	elif '4/mmm' in SytSym:
3367	ops = {'(x)':['4(x)','m(x)','m(y)','m(0+-)'], #m(0+-) for cubic m3m?
3368	'(y)':['4(y)','m(y)','m(z)','m(+0-)'], #m(+0-)
3369	'(z)':['4(z)','m(z)','m(x)','m(+-0)']} #m(+-0)
3370	elif '4mm' in SytSym:
3371	ops = {'(x)':['4(x)','m(y)','m(yz)'],'(y)':['4(y)','m(z)','m(xz)'],'(z)':['4(z)','m(x)','m(110)']}
3372	elif '422' in SytSym:
3373	ops = {'(x)':['4(x)','2(y)','2(yz)'],'(y)':['4(y)','2(z)','2(xz)'],'(z)':['4(z)','2(x)','2(110)']}
3374	elif '-4m2' in SytSym:
3375	ops = {'(x)':['-4(x)','m(x)','2(yz)'],'(y)':['-4(y)','m(y)','2(xz)'],'(z)':['-4(z)','m(z)','2(110)']}
3376	elif '-42m' in SytSym:
3377	ops = {'(x)':['-4(x)','2(y)','m(yz)'],'(y)':['-4(y)','2(z)','m(xz)'],'(z)':['-4(z)','2(x)','m(110)']}
3378	elif '-4' in SytSym:
3379	ops = {'(x)':['-4(x)',],'(y)':['-4(y)',],'(z)':['-4(z)',],}
3380	elif '4' in SytSym:
3381	ops = {'(x)':['4(x)',],'(y)':['4(y)',],'(z)':['4(z)',],}
3382
3383	elif '222' in SytSym:
3384	ops = {'':['2(x)','2(y)','2(z)'],
3385	'(x)':['2(y)','2(z)','2(x)'],'(y)':['2(x)','2(z)','2(y)'],'(z)':['2(x)','2(y)','2(z)'],
3386	'(100)':['2(z)','2(100)','2(120)',],'(010)':['2(z)','2(010)','2(210)',],
3387	'(110)':['2(z)','2(110)','2(+-0)',],}
3388	elif 'mm2' in SytSym:
3389	ops = {'(x)':['m(y)','m(z)','2(x)'],'(y)':['m(x)','m(z)','2(y)'],'(z)':['m(x)','m(y)','2(z)'],
3390	'(xy)':['m(+-0)','m(z)','2(110)'],'(yz)':['m(0+-)','m(xz)','2(yz)'], #not 2(xy)!
3391	'(xz)':['m(+0-)','m(y)','2(xz)'],'(z100)':['m(100)','m(120)','2(z)'],
3392	'(z010)':['m(010)','m(210)','2(z)'],'(z110)':['m(110)','m(+-0)','2(z)'],
3393	'(+-0)':[ 'm(110)','m(z)','2(+-0)'],'(d100)':['m(yz)','m(0+-)','2(xz)'],
3394	'(d010)':['m(xz)','m(+0-)','2(y)'],'(d001)':['m(110)','m(+-0)','2(z)'],
3395	'(210)':['m(z)','m(010)','2(210)'],'(120)':['m(z)','m(100)','2(120)'],
3396	'(100)':['m(z)','m(120)','2(100)',],'(010)':['m(z)','m(210)','2(010)',],
3397	'(110)':['m(z)','m(+-0)','2(110)',],}
3398	elif 'mmm' in SytSym:
3399	ops = {'':['m(x)','m(y)','m(z)'],
3400	'(100)':['m(z)','m(100)','m(120)',],'(010)':['m(z)','m(010)','m(210)',],
3401	'(110)':['m(z)','m(110)','m(+-0)',],
3402	'(x)':['m(x)','m(y)','m(z)'],'(y)':['m(x)','m(y)','m(z)'],'(z)':['m(x)','m(y)','m(z)'],}
3403
3404	elif '32' in SytSym:
3405	ops = {'(120)':['3','2(120)',],'(100)':['3','2(100)'],'(111)':['3(111)','2(x)']}
3406	elif '23' in SytSym:
3407	ops = {'':['2(x)','3(111)']}
3408	elif 'm3' in SytSym:
3409	ops = {'(100)':['(+-0)',],'(+--)':[],'(-+-)':[],'(--+)':[]}
3410	elif '3m' in SytSym:
3411	ops = {'(111)':['3(111)','m(+-0)',],'(+--)':['3(+--)','m(0+-)',],
3412	'(-+-)':['3(-+-)','m(+0-)',],'(--+)':['3(--+)','m(+-0)',],
3413	'(100)':['3','m(100)'],'(120)':['3','m(210)',]}
3414
3415	if SytSym.split('(')[0] in ['6/m','6mm','-6m2','622','-6','-3','-3m','-43m',]: #not simple cases
3416	MagSytSym = SytSym
3417	if "-1'" in dupDir:
3418	if '-6' in SytSym:
3419	MagSytSym = MagSytSym.replace('-6',"-6'")
3420	elif '-3m' in SytSym:
3421	MagSytSym = MagSytSym.replace('-3m',"-3'm'")
3422	elif '-3' in SytSym:
3423	MagSytSym = MagSytSym.replace('-3',"-3'")
3424	elif '-6m2' in SytSym:
3425	if "m'(110)" in dupDir:
3426	MagSytSym = "-6m'2'("+SytSym.split('(')[1]
3427	elif '6/m' in SytSym:
3428	if "m'(z)" in dupDir:
3429	MagSytSym = "6'/m'"
3430	elif '6mm' in SytSym:
3431	if "m'(110)" in dupDir:
3432	MagSytSym = "6'm'm"
3433	elif '-43m' in SytSym:
3434	if "m'(110)" in dupDir:
3435	MagSytSym = "-43m'"
3436	return MagSytSym
3437	try:
3438	axis = '('+SytSym.split('(')[1]
3439	except IndexError:
3440	axis = ''
3441	MagSytSym = ''
3442	for m in ops[axis]:
3443	if m in dupDir:
3444	MagSytSym += m.split('(')[0]
3445	else:
3446	MagSytSym += m.split('(')[0]+"'"
3447	if '2/m' in SytSym and '2' in m:
3448	MagSytSym += '/'
3449	if '-3/m' in SytSym:
3450	MagSytSym = '-'+MagSytSym
3451
3452	MagSytSym += axis
3453	# some exceptions & special rules
3454	if MagSytSym == "4'/m'm'm'": MagSytSym = "4/m'm'm'"
3455	return MagSytSym
3456
3457	# if len(GenSym) == 3:
3458	# if SGSpin[1] < 0:
3459	# if 'mm2' in SytSym:
3460	# MagSytSym = "m'm'2"+'('+SplitSytSym[1]
3461	# else: #bad rule for I41/a
3462	# MagSytSym = SplitSytSym[0]+"'"
3463	# if len(SplitSytSym) > 1:
3464	# MagSytSym += '('+SplitSytSym[1]
3465	# else:
3466	# MagSytSym = SytSym
3467	# if len(SplitSytSym) >1:
3468	# if "-4'"+'('+SplitSytSym[1] in dupDir:
3469	# MagSytSym = MagSytSym.replace('-4',"-4'")
3470	# if "-6'"+'('+SplitSytSym[1] in dupDir:
3471	# MagSytSym = MagSytSym.replace('-6',"-6'")
3472	# return MagSytSym
3473	#
3474	return SytSym
3475
3476	def UpdateSytSym(Phase):
3477	''' Update site symmetry/site multiplicity after space group/BNS lattice change
3478	'''
3479	generalData = Phase['General']
3480	SGData = generalData['SGData']
3481	Atoms = Phase['Atoms']
3482	cx,ct,cs,cia = generalData['AtomPtrs']
3483	for atom in Atoms:
3484	XYZ = atom[cx:cx+3]
3485	sytsym,Mult = SytSym(XYZ,SGData)[:2]
3486	sytSym,Mul,Nop,dupDir = SytSym(XYZ,SGData)
3487	atom[cs] = sytsym
3488	if generalData['Type'] == 'magnetic':
3489	magSytSym = MagSytSym(sytSym,dupDir,SGData)
3490	atom[cs] = magSytSym
3491	atom[cs+1] = Mult
3492	return
3493
3494	def ElemPosition(SGData):
3495	''' Under development.
3496	Object here is to return a list of symmetry element types and locations suitable
3497	for say drawing them.
3498	So far I have the element type... getting all possible locations without lookup may be impossible!
3499	'''
3500	Inv = SGData['SGInv']
3501	eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
3502	# get operators & expand if centrosymmetric
3503	SymElements = []
3504	Ops = SGData['SGOps']
3505	opM = np.array([op[0].T for op in Ops])
3506	opT = np.array([op[1] for op in Ops])
3507	if Inv:
3508	opM = np.concatenate((opM,-opM))
3509	opT = np.concatenate((opT,-opT))
3510	opMT = list(zip(opM,opT))
3511	for M,T in opMT[1:]: #skip I
3512	Dt = int(nl.det(M))
3513	Tr = int(np.trace(M))
3514	Dt = -(Dt-1)//2
3515	Es = eleSym[Tr][Dt]
3516	if Dt: #rotation-inversion
3517	I = np.eye(3)
3518	if Tr == 1: #mirrors/glides
3519	if np.any(T): #glide
3520	M2 = np.inner(M,M)
3521	MT = np.inner(M,T)+T
3522	print ('glide',Es,MT)
3523	print (M2)
3524	else: #mirror
3525	print ('mirror',Es,T)
3526	print (I-M)
3527	X = [-1,-1,-1]
3528	elif Tr == -3: # pure inversion
3529	X = np.inner(nl.inv(I-M),T)
3530	print ('inversion',Es,X)
3531	else: #other rotation-inversion
3532	M2 = np.inner(M,M)
3533	MT = np.inner(M,T)+T
3534	print ('rot-inv',Es,MT)
3535	print (M2)
3536	X = [-1,-1,-1]
3537	else: #rotations
3538	print ('rotation',Es)
3539	X = [-1,-1,-1]
3540	SymElements.append([Es,X])
3541
3542	return SymElements
3543
3544	def ApplyStringOps(A,SGData,X,Uij=[]):
3545	'Needs a doc string'
3546	SGOps = SGData['SGOps']
3547	SGCen = SGData['SGCen']
3548	Ax = A.split('+')
3549	Ax[0] = int(Ax[0])
3550	iC = 1
3551	if Ax[0] < 0:
3552	iC = -1
3553	Ax[0] = abs(Ax[0])
3554	nA = Ax[0]%100-1
3555	cA = Ax[0]//100
3556	Cen = SGCen[cA]
3557	M,T = SGOps[nA]
3558	if len(Ax)>1:
3559	cellA = Ax[1].split(',')
3560	cellA = np.array([int(a) for a in cellA])
3561	else:
3562	cellA = np.zeros(3)
3563	newX = Cen+iC*(np.inner(M,X).T+T)+cellA
3564	if len(Uij):
3565	U = Uij2U(Uij)
3566	U = np.inner(M,np.inner(U,M).T)
3567	newUij = U2Uij(U)
3568	return [newX,newUij]
3569	else:
3570	return newX
3571
3572	def ApplyStringOpsMom(A,SGData,SSGData,Mom):
3573	'''Applies string operations to modulated magnetic moment components used in drawing
3574	Drawing matches Bilbao MVISUALIZE
3575	'''
3576	SGOps = SGData['SGOps']
3577	Ax = A.split('+')
3578	Ax[0] = int(Ax[0])
3579	iAx = abs(Ax[0])
3580	nA = iAx%100-1
3581	if SGData['SGInv'] and not SGData['SGFixed']:
3582	nC = 2len(SGOps)(iAx//100)
3583	else:
3584	nC = len(SGOps)*(iAx//100)
3585	NA = nA
3586	if Ax[0] < 0:
3587	NA += len(SGOps)
3588	M,T = SGOps[nA]
3589	newMom = np.inner(Mom,M).T*SGData['SpnFlp'][NA+nC]
3590	if SSGData is not None:
3591	if SSGData['SSGCen'][iAx//100][3]: #flip spin for BNS centered atoms
3592	newMom *= -1.
3593	return newMom
3594
3595	def StringOpsProd(A,B,SGData):
3596	"""
3597	Find AB where A & B are in strings '-' + '100c+n' + '+ijk'
3598	where '-' indicates inversion, c(>0) is the cell centering operator,
3599	n is operator number from SgOps and ijk are unit cell translations (each may be <0).
3600	Should return resultant string - C. SGData - dictionary using entries:
3601
3602	* 'SGCen': cell centering vectors [0,0,0] at least
3603	* 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
3604
3605	"""
3606	SGOps = SGData['SGOps']
3607	SGCen = SGData['SGCen']
3608	#1st split out the cell translation part & work on the operator parts
3609	Ax = A.split('+'); Bx = B.split('+')
3610	Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
3611	iC = 0
3612	if Ax[0]*Bx[0] < 0:
3613	iC = 1
3614	Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
3615	nA = Ax[0]%100-1; nB = Bx[0]%100-1
3616	cA = Ax[0]//100; cB = Bx[0]//100
3617	Cen = (SGCen[cA]+SGCen[cB])%1.0
3618	cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
3619	Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
3620	Mc = np.inner(Ma,Mb.T)
3621	# print Ma,Mb,Mc
3622	Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
3623	# print Ta,Tb,Tc
3624	# print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
3625	nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
3626	#now the cell translation part
3627	if len(Ax)>1:
3628	cellA = Ax[1].split(',')
3629	cellA = [int(a) for a in cellA]
3630	else:
3631	cellA = [0,0,0]
3632	if len(Bx)>1:
3633	cellB = Bx[1].split(',')
3634	cellB = [int(b) for b in cellB]
3635	else:
3636	cellB = [0,0,0]
3637	cellC = np.add(cellA,cellB)
3638	C = str(((nC+1)+(100cC))(1-2*iC))+'+'+ \
3639	str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
3640	return C
3641
3642	def U2Uij(U):
3643	#returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
3644	return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
3645
3646	def Uij2U(Uij):
3647	#returns the thermal motion tensor U from Uij as numpy array
3648	return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
3649
3650	def StandardizeSpcName(spcgroup):
3651	'''Accept a spacegroup name where spaces may have not been used
3652	in the names according to the GSAS convention (spaces between symmetry
3653	for each axis) and return the space group name as used in GSAS
3654	'''
3655	rspc = spcgroup.replace(' ','').upper()
3656	# deal with rhombohedral and hexagonal setting designations
3657	rhomb = ''
3658	if rspc[-1:] == 'R':
3659	rspc = rspc[:-1]
3660	rhomb = ' R'
3661	gray = ''
3662	if "1'" in rspc:
3663	gray = " 1'"
3664	rspc = rspc.replace("1'",'')
3665	rspc = rspc.replace("'",'')
3666	if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
3667	rspc = rspc[:-1]
3668	if rspc[1:3] in ['M3','N3','A3','D3']: #fix cubic old style
3669	rspc.replace('3','-3')
3670	bns = -1
3671	try:
3672	bns = rspc.index('_')
3673	rspc = rspc.replace(rspc[bns:bns+2],'')
3674	except ValueError:
3675	pass
3676	# look for a match in the spacegroup lists
3677	for i in spglist.values():
3678	for spc in i:
3679	if rspc == spc.replace(' ','').upper():
3680	return spc+gray+rhomb
3681	# how about the post-2002 orthorhombic names?
3682	if rspc in sgequiv_2002_orthorhombic:
3683	return sgequiv_2002_orthorhombic[rspc]+gray
3684	else:
3685	# not found
3686	return ''
3687
3688	def SpaceGroupNumber(spcgroup):
3689	SGNo = -1
3690	SpcGp = StandardizeSpcName(spcgroup)
3691	if not SpcGp:
3692	return SGNo
3693	try:
3694	SGNo = spgbyNum.index(SpcGp)
3695	except ValueError:
3696	pass
3697	return SGNo
3698
3699
3700	spgbyNum = []
3701	'''Space groups indexed by number'''
3702	spgbyNum = [None,
3703	'P 1','P -1', #1-2
3704	'P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
3705	'C 2/m','P 2/c','P 21/c','C 2/c', #3-15
3706	'P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21',
3707	'C 2 2 21','C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
3708	'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
3709	'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
3710	'C m m 2','C m c 21','C c c 2',
3711	'A m m 2','A b m 2','A m a 2','A b a 2',
3712	'F m m 2','F d d 2','I m m 2','I b a 2','I m a 2',
3713	'P m m m','P n n n','P c c m','P b a n',
3714	'P m m a','P n n a','P m n a','P c c a','P b a m','P c c n',
3715	'P b c m','P n n m','P m m n','P b c n','P b c a','P n m a',
3716	'C m c m','C m c a','C m m m','C c c m','C m m a','C c c a',
3717	'F m m m', 'F d d d',
3718	'I m m m','I b a m','I b c a','I m m a', #16-74
3719	'P 4','P 41','P 42','P 43',
3720	'I 4','I 41',
3721	'P -4','I -4','P 4/m','P 42/m','P 4/n','P 42/n',
3722	'I 4/m','I 41/a',
3723	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3724	'P 42 21 2','P 43 2 2','P 43 21 2',
3725	'I 4 2 2','I 41 2 2',
3726	'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c','P 4 n c',
3727	'P 42 m c','P 42 b c',
3728	'I 4 m m','I 4 c m','I 41 m d','I 41 c d',
3729	'P -4 2 m','P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2',
3730	'P -4 c 2','P -4 b 2','P -4 n 2',
3731	'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d',
3732	'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',
3733	'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',
3734	'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',
3735	'P 42/n c m',
3736	'I 4/m m m','I 4/m c m','I 41/a m d','I 41/a c d',
3737	'P 3','P 31','P 32','R 3','P -3','R -3',
3738	'P 3 1 2','P 3 2 1','P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1',
3739	'R 3 2',
3740	'P 3 m 1','P 3 1 m','P 3 c 1','P 3 1 c',
3741	'R 3 m','R 3 c',
3742	'P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1',
3743	'R -3 m','R -3 c', #75-167
3744	'P 6','P 61',
3745	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3746	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3747	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3748	'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
3749	'P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m 3','P n 3',
3750	'F m -3','F d -3','I m -3',
3751	'P a -3','I a -3','P 4 3 2','P 42 3 2','F 4 3 2','F 41 3 2',
3752	'I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2','P -4 3 m',
3753	'F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
3754	'P m -3 m','P n -3 n','P m -3 n','P n -3 m',
3755	'F m -3 m','F m -3 c','F d -3 m','F d -3 c',
3756	'I m -3 m','I a -3 d',] #195-230
3757	altSettingOrtho = {}
3758	''' A dictionary of alternate settings for orthorhombic unit cells
3759	'''
3760	altSettingOrtho = {
3761	'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'},
3762	'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'},
3763	'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'},
3764	'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'},
3765	'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'},
3766	'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'},
3767	'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'},
3768	'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'},
3769	'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'},
3770	'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'},
3771	'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'},
3772	'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'},
3773	'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'},
3774	'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'},
3775	'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'},
3776	'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'},
3777	'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'},
3778	'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'},
3779	'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'},
3780	'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'},
3781	'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'},
3782	'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'},
3783	'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'},
3784	'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'},
3785	'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'},
3786	'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'},
3787	'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'},
3788	'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'},
3789	'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'},
3790	'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'},
3791	'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'},
3792	'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'},
3793	'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'},
3794	'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'},
3795	'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'},
3796	'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'},
3797	'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'},
3798	'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'},
3799	'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'},
3800	'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'},
3801	'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'},
3802	'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'},
3803	'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'},
3804	'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'},
3805	'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'},
3806	'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'},
3807	'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'},
3808	'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'},
3809	'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'},
3810	'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'},
3811	'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'},
3812	'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'},
3813	'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'},
3814	'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'},
3815	'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'},
3816	'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'},
3817	'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'},
3818	'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'},
3819	'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'},
3820	}
3821	spg2origins = {}
3822	''' A dictionary of all spacegroups that have 2nd settings; the value is the
3823	1st --> 2nd setting transformation vector as X(2nd) = X(1st)-V, nonstandard ones are included.
3824	'''
3825	spg2origins = {
3826	'P n n n':[-.25,-.25,-.25],
3827	'P b a n':[-.25,-.25,0],'P n c b':[0,-.25,-.25],'P c n a':[-.25,0,-.25],
3828	'P m m n':[-.25,-.25,0],'P n m m':[0,-.25,-.25],'P m n m':[-.25,0,-.25],
3829	'C c c a':[0,-.25,-.25],'C c c b':[-.25,0,-.25],'A b a a':[-.25,0,-.25],
3830	'A c a a':[-.25,-.25,0],'B b c b':[-.25,-.25,0],'B b a b':[0,-.25,-.25],
3831	'F d d d':[-.125,-.125,-.125],
3832	'P 4/n':[-.25,-.25,0],'P 42/n':[-.25,-.25,-.25],'I 41/a':[0,-.25,-.125],
3833	'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],
3834	'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],
3835	'I 41/a m d':[0,.25,-.125],'I 41/a c d':[0,.25,-.125],
3836	'p n -3':[-.25,-.25,-.25],'F d -3':[-.125,-.125,-.125],'P n -3 n':[-.25,-.25,-.25],
3837	'P n -3 m':[-.25,-.25,-.25],'F d -3 m':[-.125,-.125,-.125],'F d -3 c':[-.375,-.375,-.375],
3838	'p n 3':[-.25,-.25,-.25],'F d 3':[-.125,-.125,-.125],'P n 3 n':[-.25,-.25,-.25],
3839	'P n 3 m':[-.25,-.25,-.25],'F d 3 m':[-.125,-.125,-.125],'F d - c':[-.375,-.375,-.375]}
3840	spglist = {}
3841	'''A dictionary of space groups as ordered and named in the pre-2002 International
3842	Tables Volume A, except that spaces are used following the GSAS convention to
3843	separate the different crystallographic directions.
3844	Note that the symmetry codes here will recognize many non-standard space group
3845	symbols with different settings. They are ordered by Laue group
3846	'''
3847	spglist = {
3848	'P1' : ('P 1','P -1',), # 1-2
3849	'C1' : ('C 1','C -1',),
3850	'P2/m': ('P 2','P 21','P m','P a','P c','P n',
3851	'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
3852	'C2/m':('C 2','C m','C c','C n',
3853	'C 2/m','C 2/c','C 2/n',),
3854	'A2/m':('A 2','A m','A a','A n',
3855	'A 2/m','A 2/a','A 2/n',),
3856	'I2/m':('I 2','I m','I a','I n','I c',
3857	'I 2/m','I 2/a','I 2/c','I 2/n',),
3858	'Pmmm':('P 2 2 2',
3859	'P 2 2 21','P 21 2 2','P 2 21 2',
3860	'P 21 21 2','P 2 21 21','P 21 2 21',
3861	'P 21 21 21',
3862	'P m m 2','P 2 m m','P m 2 m',
3863	'P m c 21','P 21 m a','P b 21 m','P m 21 b','P c m 21','P 21 a m',
3864	'P c c 2','P 2 a a','P b 2 b',
3865	'P m a 2','P 2 m b','P c 2 m','P m 2 a','P b m 2','P 2 c m',
3866	'P c a 21','P 21 a b','P c 21 b','P b 21 a','P b c 21','P 21 c a',
3867	'P n c 2','P 2 n a','P b 2 n','P n 2 b','P c n 2','P 2 a n',
3868	'P m n 21','P 21 m n','P n 21 m','P m 21 n','P n m 21','P 21 n m',
3869	'P b a 2','P 2 c b','P c 2 a',
3870	'P n a 21','P 21 n b','P c 21 n','P n 21 a','P b n 21','P 21 c n',
3871	'P n n 2','P 2 n n','P n 2 n',
3872	'P m m m','P n n n',
3873	'P c c m','P m a a','P b m b',
3874	'P b a n','P n c b','P c n a',
3875	'P m m a','P b m m','P m c m','P m a m','P m m b','P c m m',
3876	'P n n a','P b n n','P n c n','P n a n','P n n b','P c n n',
3877	'P m n a','P b m n','P n c m','P m a n','P n m b','P c n m',
3878	'P c c a','P b a a','P b c b','P b a b','P c c b','P c a a',
3879	'P b a m','P m c b','P c m a',
3880	'P c c n','P n a a','P b n b',
3881	'P b c m','P m c a','P b m a','P c m b','P c a m','P m a b',
3882	'P n n m','P m n n','P n m n',
3883	'P m m n','P n m m','P m n m',
3884	'P b c n','P n c a','P b n a','P c n b','P c a n','P n a b',
3885	'P b c a','P c a b',
3886	'P n m a','P b n m','P m c n','P n a m','P m n b','P c m n',
3887	),
3888	'Cmmm':('C 2 2 21','C 2 2 2','C m m 2',
3889	'C m c 21','C c m 21','C c c 2','C m 2 m','C 2 m m',
3890	'C m 2 a','C 2 m b','C c 2 m','C 2 c m','C c 2 a','C 2 c b',
3891	'C m c m','C c m m','C m c a','C c m b',
3892	'C m m m','C c c m','C m m a','C m m b','C c c a','C c c b',),
3893	'Ammm':('A 21 2 2','A 2 2 2','A 2 m m',
3894	'A 21 m a','A 21 a m','A 2 a a','A m 2 m','A m m 2',
3895	'A b m 2','A c 2 m','A m a 2','A m 2 a','A b a 2','A c 2 a',
3896	'A m m a','A m a m','A b m a','A c a m',
3897	'A m m m','A m a a','A b m m','A c m m','A c a a','A b a a',),
3898	'Bmmm':('B 2 21 2','B 2 2 2','B m 2 m',
3899	'B m 21 b','B b 21 m','B b 2 b','B m m 2','B 2 m m',
3900	'B 2 c m','B m a 2','B 2 m b','B b m 2','B 2 c b','B b a 2',
3901	'B b m m','B m m b','B b c m','B m a b',
3902	'B m m m','B b m b','B m a m','B m c m','B b a b','B b c b',),
3903	'Immm':('I 2 2 2','I 21 21 21',
3904	'I m m 2','I m 2 m','I 2 m m',
3905	'I b a 2','I 2 c b','I c 2 a',
3906	'I m a 2','I 2 m b','I c 2 m','I m 2 a','I b m 2','I 2 c m',
3907	'I m m m','I b a m','I m c b','I c m a',
3908	'I b c a','I c a b',
3909	'I m m a','I b m m ','I m c m','I m a m','I m m b','I c m m',),
3910	'Fmmm':('F 2 2 2','F m m m', 'F d d d',
3911	'F m m 2','F m 2 m','F 2 m m',
3912	'F d d 2','F d 2 d','F 2 d d',),
3913	'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
3914	'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
3915	'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
3916	'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
3917	'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
3918	'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',
3919	'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',
3920	'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',
3921	'P 42/n c m',),
3922	'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',
3923	'I 4 c m','I 41 m d','I 41 c d',
3924	'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',
3925	'I 41/a m d','I 41/a c d'),
3926	'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
3927	'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
3928	'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
3929	'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
3930	'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
3931	'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
3932	'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
3933	'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
3934	'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',
3935	'P 4 3 2','P 42 3 2','P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n',
3936	'P m 3 m','P m -3 m','P n 3 n','P n -3 n',
3937	'P m 3 n','P m -3 n','P n 3 m','P n -3 m',),
3938	'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',
3939	'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'),
3940	'Fm3m':('F 2 3','F m 3','F m -3','F d 3','F d -3',
3941	'F 4 3 2','F 41 3 2','F -4 3 m','F -4 3 c',
3942	'F m 3 m','F m -3 m','F m 3 c','F m -3 c',
3943	'F d 3 m','F d -3 m','F d 3 c','F d -3 c',),
3944	}
3945	sgequiv_2002_orthorhombic = {}
3946	''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
3947	along with the pre-2002 name. The e designates a double glide-plane
3948	'''
3949	sgequiv_2002_orthorhombic = {
3950	'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
3951	'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
3952	'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
3953	'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
3954	'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
3955	'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
3956	'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
3957	'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
3958
3959	#'A few non-standard space groups for test use'
3960	nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h',
3961	'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
3962	'R 3 c r','R -3 c r','R -3 m r',),
3963
3964	#Use the space groups types in this order to list the symbols in the
3965	#order they are listed in the International Tables, vol. A'''
3966	symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal',
3967	'trigonal', 'hexagonal', 'cubic')
3968
3969	# self-test materials follow. Requires files in directory testinp
3970	selftestlist = []
3971	'''Defines a list of self-tests'''
3972	selftestquiet = True
3973	def _ReportTest():
3974	'Report name and doc string of current routine when ``selftestquiet`` is False'
3975	if not selftestquiet:
3976	import inspect
3977	caller = inspect.stack()[1][3]
3978	doc = eval(caller).__doc__
3979	if doc is not None:
3980	print('testing '+__file__+' with '+caller+' ('+doc+')')
3981	else:
3982	print('testing '+__file__()+" with "+caller)
3983	def test0():
3984	'''self-test #0: exercise MoveToUnitCell'''
3985	_ReportTest()
3986	msg = "MoveToUnitCell failed"
3987	assert (MoveToUnitCell([1,2,3])[0] == [0,0,0]).all, msg
3988	assert (MoveToUnitCell([2,-1,-2])[0] == [0,0,0]).all, msg
3989	assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9)[0] < 1e-6, msg
3990	assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1)[0] < 1e-6, msg
3991	selftestlist.append(test0)
3992
3993	def test1():
3994	'''self-test #1: SpcGroup against previous results'''
3995	#'''self-test #1: SpcGroup and SGPrint against previous results'''
3996	_ReportTest()
3997	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
3998	if ospath.exists(testdir):
3999	if testdir not in sys.path: sys.path.insert(0,testdir)
4000	import spctestinp
4001	def CompareSpcGroup(spc, referr, refdict, reflist):
4002	'Compare output from GSASIIspc.SpcGroup with results from a previous run'
4003	# if an error is reported, the dictionary can be ignored
4004	msg0 = "CompareSpcGroup failed on space group %s" % spc
4005	result = SpcGroup(spc)
4006	if result[0] == referr and referr > 0: return True
4007	# #print result[1]['SpGrp']
4008	#msg = msg0 + " in list lengths"
4009	#assert len(keys) == len(refdict.keys()), msg
4010	for key in refdict.keys():
4011	if key == 'SGOps' or key == 'SGCen':
4012	msg = msg0 + (" in key %s length" % key)
4013	assert len(refdict[key]) == len(result[1][key]), msg
4014	for i in range(len(refdict[key])):
4015	msg = msg0 + (" in key %s level 0" % key)
4016	assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
4017	msg = msg0 + (" in key %s level 1" % key)
4018	assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
4019	else:
4020	msg = msg0 + (" in key %s" % key)
4021	assert result[1][key] == refdict[key], msg
4022	msg = msg0 + (" in key %s reflist" % key)
4023	#for (l1,l2) in zip(reflist, SGPrint(result[1])):
4024	# assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
4025	# for now disable SGPrint testing, output has changed
4026	#assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
4027	for spc in spctestinp.SGdat:
4028	CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
4029	selftestlist.append(test1)
4030
4031	def test2():
4032	'''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
4033	_ReportTest()
4034	testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
4035	if ospath.exists(testdir):
4036	if testdir not in sys.path: sys.path.insert(0,testdir)
4037	import sgtbxtestinp
4038	def CompareWcctbx(spcname, cctbx_in, debug=0):
4039	'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
4040	cctbx = cctbx_in[:] # make copy so we don't delete from the original
4041	spc = (SpcGroup(spcname))[1]
4042	if debug: print (spc['SpGrp'])
4043	if debug: print (spc['SGCen'])
4044	latticetype = spcname.strip().upper()[0]
4045	# lattice type of R implies Hexagonal centering", fix the rhombohedral settings
4046	if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
4047	assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
4048	onebar = [1]
4049	if spc['SGInv']: onebar.append(-1)
4050	for (op,off) in spc['SGOps']:
4051	for inv in onebar:
4052	for cen in spc['SGCen']:
4053	noff = off + cen
4054	noff = MoveToUnitCell(noff)[0]
4055	mult = tuple((op*inv).ravel().tolist())
4056	if debug: print ("\n%s: %s + %s" % (spcname,mult,noff))
4057	for refop in cctbx:
4058	if debug: print (refop)
4059	# check the transform
4060	if refop[:9] != mult: continue
4061	if debug: print ("mult match")
4062	# check the translation
4063	reftrans = list(refop[-3:])
4064	reftrans = MoveToUnitCell(reftrans)[0]
4065	if all(abs(noff - reftrans) < 1.e-5):
4066	cctbx.remove(refop)
4067	break
4068	else:
4069	assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
4070	for key in sgtbxtestinp.sgtbx:
4071	CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
4072	selftestlist.append(test2)
4073
4074	def test3():
4075	'''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
4076	for selected space groups against info in IT Volume A '''
4077	_ReportTest()
4078	def ExerciseSiteSym (spc, crdlist):
4079	'compare site symmetries and multiplicities for a specified space group'
4080	msg = "failed on site sym test for %s" % spc
4081	(E,S) = SpcGroup(spc)
4082	assert not E, msg
4083	for t in crdlist:
4084	symb, m, n, od = SytSym(t[0],S)
4085	if symb.strip() != t[2].strip() or m != t[1]:
4086	print (spc,t[0],m,n,symb,t[2],od)
4087	assert m == t[1]
4088	#assert symb.strip() == t[2].strip()
4089
4090	ExerciseSiteSym('p 1',[
4091	((0.13,0.22,0.31),1,'1'),
4092	((0,0,0),1,'1'),
4093	])
4094	ExerciseSiteSym('p -1',[
4095	((0.13,0.22,0.31),2,'1'),
4096	((0,0.5,0),1,'-1'),
4097	])
4098	ExerciseSiteSym('C 2/c',[
4099	((0.13,0.22,0.31),8,'1'),
4100	((0.0,.31,0.25),4,'2(y)'),
4101	((0.25,.25,0.5),4,'-1'),
4102	((0,0.5,0),4,'-1'),
4103	])
4104	ExerciseSiteSym('p 2 2 2',[
4105	((0.13,0.22,0.31),4,'1'),
4106	((0,0.5,.31),2,'2(z)'),
4107	((0.5,.31,0.5),2,'2(y)'),
4108	((.11,0,0),2,'2(x)'),
4109	((0,0.5,0),1,'222'),
4110	])
4111	ExerciseSiteSym('p 4/n',[
4112	((0.13,0.22,0.31),8,'1'),
4113	((0.25,0.75,.31),4,'2(z)'),
4114	((0.5,0.5,0.5),4,'-1'),
4115	((0,0.5,0),4,'-1'),
4116	((0.25,0.25,.31),2,'4(001)'),
4117	((0.25,.75,0.5),2,'-4(001)'),
4118	((0.25,.75,0.0),2,'-4(001)'),
4119	])
4120	ExerciseSiteSym('p 31 2 1',[
4121	((0.13,0.22,0.31),6,'1'),
4122	((0.13,0.0,0.833333333),3,'2(100)'),
4123	((0.13,0.13,0.),3,'2(110)'),
4124	])
4125	ExerciseSiteSym('R 3 c',[
4126	((0.13,0.22,0.31),18,'1'),
4127	((0.0,0.0,0.31),6,'3'),
4128	])
4129	ExerciseSiteSym('R 3 c R',[
4130	((0.13,0.22,0.31),6,'1'),
4131	((0.31,0.31,0.31),2,'3(111)'),
4132	])
4133	ExerciseSiteSym('P 63 m c',[
4134	((0.13,0.22,0.31),12,'1'),
4135	((0.11,0.22,0.31),6,'m(100)'),
4136	((0.333333,0.6666667,0.31),2,'3m(100)'),
4137	((0,0,0.31),2,'3m(100)'),
4138	])
4139	ExerciseSiteSym('I a -3',[
4140	((0.13,0.22,0.31),48,'1'),
4141	((0.11,0,0.25),24,'2(x)'),
4142	((0.11,0.11,0.11),16,'3(111)'),
4143	((0,0,0),8,'-3(111)'),
4144	])
4145	selftestlist.append(test3)
4146
4147	if __name__ == '__main__':
4148	# run self-tests
4149	selftestquiet = False
4150	for test in selftestlist:
4151	test()
4152	print ("OK")

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

Download in other formats: