source: trunk/GSASIIspc.py @ 3194

Last change on this file since 3194 was 3194, checked in by vondreele, 4 years ago

fix to CifFile?
fixes to cbf & tif importers for new Pilatus 2M detector
fixes to mag structure display

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