source: trunk/GSASIIspc.py @ 3425

Last change on this file since 3425 was 3425, checked in by vondreele, 5 years ago

use nxs = np.newaxis in G2spc.
Add 2 functions for checking symm extinct reflections (not in use yet)
fix wx.Colour object showing up in gpx file problem - convert to list[:3] (takes off alpha)

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