source: trunk/GSASIIspc.py @ 1547

Last change on this file since 1547 was 1547, checked in by vondreele, 9 years ago

change "Crystal size" to "Domain size"
use lookups for allowed super symmetries for GSAS-II standard space groups as used in the indexing routine.
supersymmetry for trigonal/rhomahedral & hexagonal all complete & checked
supersylmmetry for monoclinic & tetragonal all complete but get errors
orthorhombic not done yet
min Nc/No? set to 2 (at 1 it failed) for indexing
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 97.0 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: 2014-10-29 11:56:39 +0000 (Wed, 29 Oct 2014) $
12# $Author: vondreele $
13# $Revision: 1547 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 1547 2014-10-29 11:56:39Z vondreele $
16########### SVN repository information ###################
17import numpy as np
18import numpy.ma as ma
19import numpy.linalg as nl
20import scipy.optimize as so
21import math
22import sys
23import os.path as ospath
24
25import GSASIIpath
26GSASIIpath.SetVersionNumber("$Revision: 1547 $")
27import pyspg
28
29npsind = lambda x: np.sin(x*np.pi/180.)
30npcosd = lambda x: np.cos(x*np.pi/180.)
31   
32################################################################################
33#### Space group codes
34################################################################################
35
36def SpcGroup(SGSymbol):
37    """
38    Determines cell and symmetry information from a short H-M space group name
39
40    :param SGSymbol: space group symbol (string) with spaces between axial fields
41    :returns: (SGError,SGData)
42       * 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)
58                - 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    SGInfo = pyspg.sgforpy(SGSymbol)
69    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
70    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
71    SGData['SGInv'] = bool(SGInfo[1])
72    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
73    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
74    if SGData['SGLatt'] == 'P':
75        SGData['SGCen'] = np.array(([0,0,0],))
76    elif SGData['SGLatt'] == 'A':
77        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
78    elif SGData['SGLatt'] == 'B':
79        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
80    elif SGData['SGLatt'] == 'C':
81        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
82    elif SGData['SGLatt'] == 'I':
83        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
84    elif SGData['SGLatt'] == 'F':
85        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
86    elif SGData['SGLatt'] == 'R':
87        SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
88    SGData['SGOps'] = []
89    for i in range(SGInfo[5]):
90        Mat = np.array(SGInfo[6][i])
91        Trns = np.array(SGInfo[7][i])
92        SGData['SGOps'].append([Mat,Trns])
93    if SGData['SGLaue'] in '-1':
94        SGData['SGSys'] = SysSym[0]
95    elif SGData['SGLaue'] in '2/m':
96        SGData['SGSys'] = SysSym[1]
97    elif SGData['SGLaue'] in 'mmm':
98        SGData['SGSys'] = SysSym[2]
99    elif SGData['SGLaue'] in ['4/m','4/mmm']:
100        SGData['SGSys'] = SysSym[3]
101    elif SGData['SGLaue'] in ['3R','3mR']:
102        SGData['SGSys'] = SysSym[4]
103    elif SGData['SGLaue'] in ['3','3m1','31m']:
104        SGData['SGSys'] = SysSym[5]
105    elif SGData['SGLaue'] in ['6/m','6/mmm']:
106        SGData['SGSys'] = SysSym[6]
107    elif SGData['SGLaue'] in ['m3','m3m']:
108        SGData['SGSys'] = SysSym[7]
109    SGData['SGPolax'] = SGpolar(SGData)
110    SGData['SGPtGrp'],SGData['SSGKl'] = SGPtGroup(SGData)
111    return SGInfo[8],SGData
112
113def SGErrors(IErr):
114    '''
115    Interprets the error message code from SpcGroup. Used in SpaceGroup.
116   
117    :param IErr: see SGError in :func:`SpcGroup`
118    :returns:
119        ErrString - a string with the error message or "Unknown error"
120    '''
121
122    ErrString = [' ',
123        'Less than 2 operator fields were found',
124        'Illegal Lattice type, not P, A, B, C, I, F or R',
125        'Rhombohedral lattice requires a 3-axis',
126        'Minus sign does not preceed 1, 2, 3, 4 or 6',
127        'Either a 5-axis anywhere or a 3-axis in field not allowed',
128        ' ',
129        'I for COMPUTED GO TO out of range.',
130        'An a-glide mirror normal to A not allowed',
131        'A b-glide mirror normal to B not allowed',
132        'A c-glide mirror normal to C not allowed',
133        'D-glide in a primitive lattice not allowed',
134        'A 4-axis not allowed in the 2nd operator field',
135        'A 6-axis not allowed in the 2nd operator field',
136        'More than 24 matrices needed to define group',
137        ' ',
138        'Improper construction of a rotation operator',
139        'Mirror following a / not allowed',
140        'A translation conflict between operators',
141        'The 2bar operator is not allowed',
142        '3 fields are legal only in R & m3 cubic groups',
143        'Syntax error. Expected I -4 3 d at this point',
144        ' ',
145        'A or B centered tetragonal not allowed',
146        ' ','unknown error in sgroup',' ',' ',' ',
147        'Illegal character in the space group symbol',
148        ]
149    try:
150        return ErrString[IErr]
151    except:
152        return "Unknown error"
153
154def SGpolar(SGData):
155    '''
156    Determine identity of polar axes if any
157    '''
158    POL = ('','x','y','x y','z','x z','y z','xyz','111')
159    NP = [1,2,4]
160    NPZ = [0,1]
161    for M,T in SGData['SGOps']:
162        for i in range(3):
163            if M[i][i] <= 0.: NP[i] = 0
164        if M[0][2] > 0: NPZ[0] = 8
165        if M[1][2] > 0: NPZ[1] = 0
166    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
167    return POL[NPol]
168   
169def SGPtGroup(SGData):
170    '''
171    Determine point group of the space group - done after space group symbol has
172    been evaluated by SpcGroup. Only short symbols are allowed
173   
174    :param SGData: from :func SpcGroup
175    returns SSGPtGrp & SSGKl (only defaults for Mono & Ortho)
176    '''
177    Flds = SGData['SpGrp'].split(' ')
178    if SGData['SGLaue'] == '-1':    #triclinic
179        if '-' in Flds[1]:
180            return '-1',[-1,]
181        else:
182            return '1',[1,]
183    elif SGData['SGLaue'] == '2/m': #monoclinic - default for 2D modulation vector
184        if '/' in SGData['SpGrp']:
185            return '2/m',[-1,1]
186        elif '2' in SGData['SpGrp']:
187            return '2',[-1,]
188        else:
189            return 'm',[1,]
190    elif SGData['SGLaue'] == 'mmm': #orthorhombic
191        if SGData['SpGrp'].count('2') == 3:
192            return '222',[1,1,1]
193        elif SGData['SpGrp'].count('2') == 1:
194            if SGData['SGPolax'] == 'x':
195                return '2mm',[1,1,1]
196            elif SGData['SGPolax'] == 'y':
197                return 'm2m',[1,1,1]
198            elif SGData['SGPolax'] == 'z':
199                return 'mm2',[1,1,1]
200        else:
201            return 'mmm',[1,1,-1]
202    elif SGData['SGLaue'] == '4/m': #tetragonal
203        if '/' in SGData['SpGrp']:
204            return '4/m',[1,-1]
205        elif '-' in Flds[1]:
206            return '-4',[-1,]
207        else:
208            return '4',[1,]
209    elif SGData['SGLaue'] == '4/mmm':
210        if '/' in SGData['SpGrp']:
211            return '4/mmm',[1,-1,1,1]
212        elif '-' in Flds[1]:
213            if '2' in Flds[2]:
214                return '-42m',[-1,-1,1]
215            else:
216                return '-4m2',[-1,1,-1]             
217        elif '2' in Flds[2:]:
218            return '422',[1,-1,-1]
219        else:
220            return '4mm',[1,1,1]
221    elif SGData['SGLaue'] in ['3','3R']:  #trigonal/rhombohedral
222        if '-' in Flds[1]:
223            return '-3',[-1,]
224        else:
225            return '3',[1,]
226    elif SGData['SGLaue'] == '3mR' or 'R' in Flds[0]:
227        if '2' in Flds[2]:
228            return '32',[1,-1]
229        elif '-' in Flds[1]:
230            return '-3m',[-1,1]
231        else:
232            return '3m',[1,1]
233    elif SGData['SGLaue'] == '3m1':
234        if '2' in Flds[2]:
235            return '321',[1,-1,1]
236        elif '-' in Flds[1]:
237            return '-3m1',[-1,1,1]
238        else:
239            return '3m1',[1,1,1]
240    elif SGData['SGLaue'] == '31m':
241        if '2' in Flds[3]:
242            return '312',[1,1,-1]
243        elif '-' in Flds[1]:
244            return '-31m',[-1,1,1]
245        else:
246            return '31m',[1,1,1]
247    elif SGData['SGLaue'] == '6/m': #hexagonal
248        if '/' in SGData['SpGrp']:
249            return '6/m',[1,-1]
250        elif '-' in SGData['SpGrp']:
251            return '-6',[-1,]
252        else:
253            return '6',[1,]
254    elif SGData['SGLaue'] == '6/mmm':
255        if '/' in SGData['SpGrp']:
256            return '6/mmm',[1,-1,1,1]
257        elif '-' in Flds[1]:
258            if '2' in Flds[2]:
259                return '-62m',[-1,-1,1]
260            else:
261                return '-6m2',[-1,1,-1]                 
262        elif '2' in Flds[2:]:
263            return '622',[1,-1,-1]
264        else:
265            return '6mm',[1,1,1]   
266    elif SGData['SGLaue'] == 'm3':      #cubic - no (3+1) supersymmetry
267        if '2' in Flds[1]:
268            return '23',[]
269        else: 
270            return 'm3',[]
271    elif SGData['SGLaue'] == 'm3m':
272        if '4' in Flds[1]:
273            if '-' in Flds[1]:
274                return '-43m',[]
275            else:
276                return '432',[]
277        else:
278            return 'm-3m',[]
279   
280def SGPrint(SGData):
281    '''
282    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
283
284    :param SGData: from :func:`SpcGroup`
285    :returns:
286        SGText - list of strings with the space group details
287        SGTable - list of strings for each of the operations
288    '''
289    Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
290    SGText = []
291    SGText.append(' Space Group: '+SGData['SpGrp'])
292    CentStr = 'centrosymmetric'
293    if not SGData['SGInv']:
294        CentStr = 'non'+CentStr
295    if SGData['SGLatt'] in 'ABCIFR':
296        SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
297    else:
298        SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower()) 
299    SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
300    if 'SGPtGrp' in SGData:         #patch
301        SGText.append(' The lattice point group is '+SGData['SGPtGrp'])
302    SGText.append(' Multiplicity of a general site is '+str(Mult))
303    if SGData['SGUniq'] in ['a','b','c']:
304        SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
305    if SGData['SGInv']:
306        SGText.append(' The inversion center is located at 0,0,0')
307    if SGData['SGPolax']:
308        SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
309    SGText.append(' ')
310    if SGData['SGLatt'] == 'P':
311        SGText.append(' The equivalent positions are:\n')
312    else:   
313        SGText.append(' The equivalent positions are:')
314        SGText.append(' ('+Latt2text(SGData['SGLatt'])+')+\n')
315    SGTable = []
316    for i,Opr in enumerate(SGData['SGOps']):
317        SGTable.append('(%2d) %s'%(i+1,MT2text(Opr)))
318    return SGText,SGTable
319
320def AllOps(SGData):
321    '''
322    Returns a list of all operators for a space group, including those for
323    centering and a center of symmetry
324   
325    :param SGData: from :func:`SpcGroup`
326    :returns: (SGTextList,offsetList,symOpList,G2oprList) where
327
328      * SGTextList: a list of strings with formatted and normalized
329        symmetry operators.
330      * offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
331        symmetry operation to the operator in SGTextList and symOpList.
332        these dx (etc.) values are added to the GSAS-II generated
333        positions to provide the positions that are generated
334        by the normalized symmetry operators.       
335      * symOpList: a list of tuples with the normalized symmetry
336        operations as (M,T) values
337        (see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
338      * G2oprList: The GSAS-II operations for each symmetry operation as
339        a tuple with (center,mult,opnum), where center is (0,0,0), (0.5,0,0),
340        (0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
341        and opnum is the number for the symmetry operation, in ``SGOps``
342        (starting with 0).
343    '''
344    SGTextList = []
345    offsetList = []
346    symOpList = []
347    G2oprList = []
348    onebar = (1,)
349    if SGData['SGInv']:
350        onebar += (-1,)
351    for cen in SGData['SGCen']:
352        for mult in onebar:
353            for j,(M,T) in enumerate(SGData['SGOps']):
354                offset = [0,0,0]
355                Tprime = (mult*T)+cen
356                for i in range(3):
357                    while Tprime[i] < 0:
358                        Tprime[i] += 1
359                        offset[i] += 1
360                    while Tprime[i] >= 1:
361                        Tprime[i] += -1
362                        offset[i] += -1
363                Opr = [mult*M,Tprime]
364                OPtxt = MT2text(Opr)
365                SGTextList.append(OPtxt.replace(' ',''))
366                offsetList.append(tuple(offset))
367                symOpList.append((mult*M,Tprime))
368                G2oprList.append((cen,mult,j))
369    return SGTextList,offsetList,symOpList,G2oprList
370   
371def MT2text(Opr):
372    "From space group matrix/translation operator returns text version"
373    XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
374    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
375    Fld = ''
376    M,T = Opr
377    for j in range(3):
378        IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4))%12
379        IK = int(round(T[j]*12))%12
380        if IK:
381            if IJ < 3:
382                Fld += (TRA[IK]+XYZ[IJ]).rjust(5)
383            else:
384                Fld += (TRA[IK]+'+'+XYZ[IJ]).rjust(5)
385        else:
386            Fld += XYZ[IJ].rjust(5)
387        if j != 2: Fld += ', '
388    return Fld
389   
390def Latt2text(Latt):
391    "From lattice type ('P',A', etc.) returns ';' delimited cell centering vectors"
392    lattTxt = {'A':'0,0,0; 0,1/2,1/2','B':'0,0,0; 1/2,0,1/2',
393        'C':'0,0,0; 1/2,1/2,0','I':'0,0,0; 1/2,1/2,1/2',
394        'F':'0,0,0; 0,1/2,1/2; 1/2,0,1/2; 1/2,1/2,0',
395        'R':'0,0,0; 1/3,2/3,2/3; 2/3,1/3,1/3','P':'0,0,0'}
396    return lattTxt[Latt]   
397       
398def SpaceGroup(SGSymbol):
399    '''
400    Print the output of SpcGroup in a nicely formatted way.
401
402    :param SGSymbol: space group symbol (string) with spaces between axial fields
403    :returns: nothing
404    '''
405    E,A = SpcGroup(SGSymbol)
406    if E > 0:
407        print SGErrors(E)
408        return
409    for l in SGPrint(A):
410        print l
411       
412################################################################################
413#### Superspace group codes
414################################################################################
415       
416def SSpcGroup(SGData,SSymbol):
417    """
418    Determines supersymmetry information from superspace group name; currently only for (3+1) superlattices
419
420    :param SGData: space group data structure as defined in SpcGroup above.
421    :param SSymbol: superspace group symbol extension (string) defining modulation direction & generator info.
422    :returns: (SSGError,SSGData)
423       * SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
424       * SSGData - is a dict (see :ref:`Superspace Group object<SSGData_table>`) with entries:
425       
426             * 'SSpGrp': superspace group symbol extension to space group symbol, accidental spaces removed
427             * 'SSGCen': 4D cell centering vectors [0,0,0,0] at least
428             * 'SSGOps': 4D symmetry operations as [M,T] so that M*x+T = x'
429
430    """
431   
432    def splitSSsym(SSymbol):
433        '''
434        Splits supersymmetry symbol into two lists of strings
435        '''
436        modsym,gensym = SSymbol.replace(' ','').split(')')
437        nfrac = modsym.count('/')
438        modsym = modsym.lstrip('(')
439        if nfrac == 0:
440            modsym = list(modsym)
441        elif nfrac == 1:
442            pos = modsym.find('/')
443            if pos == 1:
444                modsym = [modsym[:3],modsym[3],modsym[4]]
445            elif pos == 2:
446                modsym = [modsym[0],modsym[1:4],modsym[4]]
447            else:
448                modsym = [modsym[0],modsym[1],modsym[2:]]
449        else:
450            lpos = modsym.find('/')
451            rpos = modsym.rfind('/')
452            if lpos == 1 and rpos == 4:
453                modsym = [modsym[:3],modsym[3:6],modsym[6]]
454            elif lpos == 1 and rpos == 5:
455                modsym = [modsym[:3],modsym[3],modsym[4:]]
456            else:
457                modsym = [modsym[0],modsym[1:4],modsym[4:]]
458        gensym = list(gensym)
459        return modsym,gensym
460       
461    def checkModSym():
462        '''
463        Checks to see if proposed modulation form is allowed for Laue group
464        '''
465        if LaueId in [0,] and LaueModId in [0,]:
466            return True
467        elif LaueId in [1,]:
468            try:
469                if modsym.index('1/2') != ['A','B','C'].index(SGData['SGLatt']):
470                    return False
471                if 'I'.index(SGData['SGLatt']) and modsym.count('1/2') not in [0,2]:
472                    return False
473            except ValueError:
474                pass
475            if SGData['SGUniq'] == 'a' and LaueModId in [5,6,7,8,9,]:
476                return True
477            elif SGData['SGUniq'] == 'b' and LaueModId in [3,4,13,14,15,]:
478                return True
479            elif SGData['SGUniq'] == 'c' and LaueModId in [1,2,19,20,21,]:
480                return True
481        elif LaueId in [2,] and LaueModId in [i+7 for i in range(18)]:
482            try:
483                if modsym.index('1/2') != ['A','B','C'].index(SGData['SGLatt']):
484                    return False
485                if SGData['SGLatt'] in ['I','F',] and modsym.index('1/2'):
486                    return False
487            except ValueError:
488                pass
489            return True
490        elif LaueId in [3,4,] and LaueModId in [19,22,]:
491            try:
492                if SGData['SGLatt'] == 'I' and modsym.count('1/2'):
493                    return False
494            except ValueError:
495                pass
496            return True
497        elif LaueId in [7,8,9,] and LaueModId in [19,25,]:
498            if (SGData['SGLatt'] == 'R' or SGData['SGPtGrp'] in ['3m1','-3m1']) and modsym.count('1/3'):
499                return False
500            return True
501        elif LaueId in [10,11,] and LaueModId in [19,]:
502            return True
503        return False
504       
505    def fixMonoOrtho():
506        mod = ''.join(modsym).replace('1/2','0').replace('1','0')
507        if SGData['SGPtGrp'] in ['2','m']:  #OK
508            if mod in ['a00','0b0','00g']:
509                result = [i*-1 for i in SGData['SSGKl']]
510            else:
511                result = SGData['SSGKl'][:]
512            if '/' in mod:
513                return [i*-1 for i in result]
514            else:
515                return result
516        elif SGData['SGPtGrp'] == '2/m':    #OK
517            if mod in ['a00','0b0','00g']:
518                result =  SGData['SSGKl'][:]
519            else:
520                result = [i*-1 for i in SGData['SSGKl']]
521            if '/' in mod:
522                return [i*-1 for i in result]
523            else:
524                return result
525        else:   #orthorhombic
526            if SGData['SGPtGrp'] == '222':
527                return [1 if i in ['a','b','g'] else -1 for i in mod]
528            elif SGData['SGPtGrp'] == 'mm2':
529                if 'g' in mod:
530                    return [1,1,1]
531                elif 'b' in mod:
532                    return [1,-1,-1]
533                else:
534                    return [-1,1,-1]
535            elif SGData['SGPtGrp'] == 'm2m':
536                if 'b' in mod:
537                    return [1,1,1]
538                elif 'g' in mod:
539                    return [1,-1,-1]
540                else:
541                    return [-1,-1,1]               
542            elif SGData['SGPtGrp'] == '2mm':
543                if 'a' in mod:
544                    return [1,1,1]
545                elif 'b' in mod:
546                    return [-1,-1,1]
547                else:
548                    return [-1,1,-1]
549            else:
550                return [-1 if i in ['a','b','g'] else 1 for i in mod]
551               
552    def extendSSGOps(SSGOps):
553        nOps = len(SSGOps)
554        for i in range(nOps):
555            if np.allclose(SSGOps[i][0][3],np.zeros(4)):
556                continue
557            for j in range(nOps):
558                if np.allclose(SSGOps[j][0][3],np.zeros(4)):
559                    continue
560                OpC = list(SGProd(SSGOps[j],SSGOps[i]))
561                OpC[1] %= 1.
562                for k in range(nOps):
563                    OpD = SSGOps[k]
564                    if SSMT2text(OpC) == SSMT2text(OpD):
565                        continue
566                    elif np.allclose(OpC[0][:3,:3],OpD[0][:3,:3]):
567                        if np.allclose(OpD[0][3],np.zeros(4)):
568                            SSGOps[k] = OpC
569                        elif np.any([np.allclose(OpC[0][3][:3],cen) for cen in SGData['SGCen']]):   #?
570                            continue
571                        else:
572                            OpCtxt = SSMT2text(OpC).replace(' ','')
573                            OpDtxt = SSMT2text(OpD).replace(' ','')
574                            print 'OpC',OpCtxt,'OpD',OpDtxt
575                            return False,OpCtxt+' conflict with '+OpDtxt
576        return True,SSGOps
577               
578    def genSSGOps():
579        SSGOps = SSGData['SSGOps'][:]
580        iFrac = {}
581        for i,frac in enumerate(SSGData['modSymb']):
582            if frac in ['1/2','1/3','1/4','1/6','1']:
583                iFrac[i] = frac
584        print SGData['SpGrp']+SSymbol
585        print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac
586# set identity & 1,-1; triclinic
587        SSGOps[0][0][3,3] = 1.
588# expand if centrosymmetric
589        if SGData['SGInv']:
590            SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
591# monoclinic - all done
592        if SGData['SGPtGrp'] in ['2','m']:  #OK
593            SSGOps[1][0][3,3] = SSGKl[0]
594            SSGOps[1][1][3] = genQ[0]
595            for i in iFrac:
596                SSGOps[1][0][3,i] = -SSGKl[0]
597        elif SGData['SGPtGrp'] == '2/m':    #OK
598            for i,j in enumerate([1,3]):
599                SSGOps[j][0][3,3] = -SSGKl[i]
600                if genQ[i]:
601                    SSGOps[j][1][3] = genQ[i]
602                for k in iFrac:
603                    SSGOps[j][0][3,k] = SSGKl[i]
604                E,SSGOps = extendSSGOps(SSGOps)
605           
606# orthorhombic
607        elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm','mmm']:
608            for i in [0,1,2]:
609                SSGOps[i+1][0][3,3] = SSGKl[i]
610                SSGOps[i+1][1][3] = genQ[i]
611            for i in iFrac:
612                SSGOps[1][0][3,i] = -1
613            print SSMT2text(SSGOps[1]).replace(' ',''),SSMT2text(SSGOps[2]).replace(' ',''), \
614                SSMT2text(SSGOps[3]).replace(' ','')
615               
616# tetragonal
617        elif SGData['SGPtGrp'] == '4':  #OK
618            SSGOps[1][0][3,3] = SSGKl[0]
619            SSGOps[1][1][3] = genQ[0]
620            if '1/2' in SSGData['modSymb']:
621                SSGOps[1][0][3,1] = -1
622        elif SGData['SGPtGrp'] == '-4': #OK
623            SSGOps[1][0][3,3] = SSGKl[0]
624            if '1/2' in SSGData['modSymb']:
625                SSGOps[1][0][3,1] = 1
626        elif SGData['SGPtGrp'] in ['4/m',]:
627            if '1/2' in SSGData['modSymb']:
628                SSGOps[1][0][3,1] = -1
629            for i,j in enumerate([1,6]):
630                SSGOps[j][0][3,3] = SSGKl[i]
631                if genQ[i]:
632                    SSGOps[j][1][3] = genQ[i]
633                E,SSGOps = extendSSGOps(SSGOps)
634        elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]:
635            if '1/2' in SSGData['modSymb']:
636                SSGOps[1][0][3,1] = -1
637            for i,j in enumerate([1,4,5]):
638                SSGOps[j][0][3,3] = SSGKl[i]
639                if genQ[i]:
640                    SSGOps[j][1][3] = genQ[i]
641                E,SSGOps = extendSSGOps(SSGOps)
642        elif SGData['SGPtGrp'] in ['4/mmm',]:
643            if '1/2' in SSGData['modSymb']:
644                SSGOps[1][0][3,1] = -1
645            for i,j in enumerate([1,10,6,7]):
646                SSGOps[j][0][3,3] = SSGKl[i]
647                if genQ[i]:
648                    SSGOps[j][1][3] = genQ[i]
649#                for k in iFrac:
650#                    SSGOps[j][0][3,k] = SSGKl[i]
651                E,SSGOps = extendSSGOps(SSGOps)
652               
653# trigonal - all done
654        elif SGData['SGPtGrp'] == '3':  #OK
655            SSGOps[1][0][3,3] = SSGKl[0]
656            if '1/3' in SSGData['modSymb']:
657                SSGOps[1][0][3,1] = -1
658            SSGOps[1][1][3] = genQ[0]
659        elif SGData['SGPtGrp'] == '-3': #OK
660            SSGOps[1][0][3,3] = -SSGKl[0]
661            if '1/3' in SSGData['modSymb']:
662                SSGOps[1][0][3,1] = -1
663            SSGOps[1][1][3] = genQ[0]
664        elif SGData['SGPtGrp'] in ['312','3m','-3m','-3m1','3m1']:   #OK
665            if '1/3' in SSGData['modSymb']:
666                SSGOps[1][0][3,1] = -1
667            for i,j in enumerate([1,5]):
668                if SGData['SGPtGrp'] in ['3m','-3m']:
669                    SSGOps[j][0][3,3] = SSGKl[i]
670                else:                   
671                    SSGOps[j][0][3,3] = SSGKl[i+1]
672                if genQ[i]:
673                    SSGOps[j][1][3] = genQ[i]
674        elif SGData['SGPtGrp'] in ['321','32']:   #OK
675            for i,j in enumerate([1,4]):
676                SSGOps[j][0][3,3] = SSGKl[i]
677                if genQ[i]:
678                    SSGOps[j][1][3] = genQ[i]
679        elif SGData['SGPtGrp'] in ['31m','-31m']:   #OK
680            ids = [1,3]
681            if SGData['SGPtGrp'] == '-31m':
682                ids = [7,3]
683            if '1/3' in SSGData['modSymb']:
684                SSGOps[ids[0]][0][3,1] = -1
685            for i,j in enumerate(ids):
686                SSGOps[j][0][3,3] = SSGKl[i]
687                if genQ[i+1]:
688                    SSGOps[j][1][3] = genQ[i+1]
689                     
690# hexagonal - all done
691        elif SGData['SGPtGrp'] == '6':  #OK
692            SSGOps[1][0][3,3] = SSGKl[0]
693            SSGOps[1][1][3] = genQ[0]
694        elif SGData['SGPtGrp'] == '-6': #OK
695            SSGOps[1][0][3,3] = SSGKl[0]
696        elif SGData['SGPtGrp'] in ['6/m',]: #OK
697            SSGOps[1][0][3,3] = -SSGKl[1]
698            SSGOps[1][1][3] = genQ[0]
699            SSGOps[2][1][3] = genQ[1]
700        elif SGData['SGPtGrp'] in ['622','6mm','-62m','-6m2',]: #OK
701            for i,j in enumerate([1,10,11]):
702                SSGOps[j][0][3,3] = SSGKl[i]
703                if genQ[i]:
704                    SSGOps[j][1][3] = genQ[i]
705                E,SSGOps = extendSSGOps(SSGOps)
706        elif SGData['SGPtGrp'] in ['6/mmm',]: #OK
707            for i,j in enumerate([1,15,10,11]):
708                SSGOps[j][0][3,3] = SSGKl[i]
709                if genQ[i]:
710                    SSGOps[j][1][3] = genQ[i]
711                E,SSGOps = extendSSGOps(SSGOps)
712        elif SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
713            return True,SSGOps
714        E,SSGOps = extendSSGOps(SSGOps)
715        return E,SSGOps
716       
717    def specialGen(gensym):
718        sym = ''.join(gensym)
719        if SGData['SGPtGrp'] in ['2/m',] and 'n' in SGData['SpGrp']:
720            if 's' in sym:
721                gensym = 'ss'
722        if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
723            gensym = '0ss'
724        elif SGData['SGPtGrp'] in ['222',]:
725            if sym == '00s':
726                gensym = '0ss'
727            elif sym == '0s0':
728                gensym = 'ss0'
729            elif sym == 's00':
730                gensym = 's0s'
731        return gensym
732                   
733    def checkGen(gensym):
734        sym = ''.join(gensym)
735        print str(SSGKl),sym
736# monoclinic - all done
737        if str(SSGKl) == '[-1]' and sym == 's':
738            return False
739        elif SGData['SGPtGrp'] in ['2/m',]:
740            if str(SSGKl) == '[-1, 1]' and sym == '0s':
741                return False
742            elif str(SSGKl) == '[1, -1]' and sym == 's0':
743                return False
744#orthorhombic - all
745        elif SGData['SGPtGrp'] in ['222',] and sym not in ['','s00','0s0','00s']:
746            return False 
747        elif SGData['SGPtGrp'] in ['2mm','m2m','mm2','mmm'] and sym not in GenSymList[4:15]:
748            return False 
749#tetragonal - all done
750        elif SGData['SGPtGrp'] in ['4',] and sym not in ['','s','q']:
751            return False 
752        elif SGData['SGPtGrp'] in ['-4',] and sym not in ['',]:
753            return False             
754        elif SGData['SGPtGrp'] in ['4/m',] and sym not in ['','s0','q0']:
755            return False
756        elif SGData['SGPtGrp'] in ['422',] and sym not in ['','q00','s00']:
757            return False         
758        elif SGData['SGPtGrp'] in ['4mm',] and sym not in ['','ss0','s0s','0ss','qq0','qqs']:
759            return False
760        elif SGData['SGPtGrp'] in ['-4m2',] and sym not in ['','00s','00q']:
761            return False
762        elif SGData['SGPtGrp'] in ['-42m',] and sym not in ['','0s0','0q0']:
763            return False
764        elif SGData['SGPtGrp'] in ['4/mmm',] and sym not in ['','s00s','s0s0','00ss','q0q0','q0qs']:
765            return False
766#trigonal/rhombohedral - all done
767        elif SGData['SGPtGrp'] in ['3',] and sym not in ['','t']:
768            return False 
769        elif SGData['SGPtGrp'] in ['-3',] and sym not in ['',]:
770            return False 
771        elif SGData['SGPtGrp'] in ['32',] and sym not in ['','t0']:
772            return False 
773        elif SGData['SGPtGrp'] in ['321','312'] and sym not in ['','t00']:
774            return False 
775        elif SGData['SGPtGrp'] in ['3m','-3m'] and sym not in ['','0s']:
776            return False 
777        elif SGData['SGPtGrp'] in ['3m1','-3m1'] and sym not in ['','0s0']:
778            return False 
779        elif SGData['SGPtGrp'] in ['31m','-31m'] and sym not in ['','00s']:
780            return False 
781#hexagonal - all done
782        elif SGData['SGPtGrp'] in ['6',] and sym not in ['','s','h','t']:
783            return False 
784        elif SGData['SGPtGrp'] in ['-6',] and sym not in ['',]:
785            return False
786        elif SGData['SGPtGrp'] in ['6/m',] and sym not in ['','s0']:
787            return False
788        elif SGData['SGPtGrp'] in ['622',] and sym not in ['','h00','t00','s00']:
789            return False         
790        elif SGData['SGPtGrp'] in ['6mm',] and sym not in ['','ss0','s0s','0ss']:
791            return False
792        elif SGData['SGPtGrp'] in ['-6m2',] and sym not in ['','0s0']:
793            return False
794        elif SGData['SGPtGrp'] in ['-62m',] and sym not in ['','0ss']:
795            return False
796        elif SGData['SGPtGrp'] in ['6/mmm',] and sym not in ['','s00s','s0s0','00ss']:
797            return False
798        return True
799       
800    LaueModList = ['abg', 'ab0', 'ab1/2', 'a0g', 'a1/2g','0bg', '1/2bg',
801               'a00', 'a01/2', 'a1/20', 'a1/21/2', 'a01', 'a10', 
802               '0b0', '0b1/2', '1/2b0', '1/2b1/2', '0b1', '1b0',
803               '00g', '01/2g', '1/20g', '1/21/2g', '01g', '10g','1/31/3g']
804    LaueList = ['-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m']
805    GenSymList = ['','s','0s','s0','00s','0s0','s00','s0s','ss0','0ss','q00','0q0','00q','qq0','q0q','0qq',
806        'q','q0','0q','qqs','s0s0','00ss','s00s','q0q0','q0qs','t','t00','t0','h','h00']
807    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.}
808    LaueId = LaueList.index(SGData['SGLaue'])
809    if SGData['SGLaue'] in ['m3','m3m']:
810        return '(3+1) superlattices not defined for cubic space groups',None
811    elif SGData['SGLaue'] in ['3R','3mR']:
812        return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
813    try:
814        modsym,gensym = splitSSsym(SSymbol)
815    except ValueError:
816        return 'Error in superspace symbol '+SSymbol,None
817    if ''.join(gensym) not in GenSymList:
818        return 'unknown generator symbol '+''.join(gensym),None
819    try:
820        print modsym,''.join(modsym)
821        LaueModId = LaueModList.index(''.join(modsym))
822    except ValueError:
823        return 'Unknown modulation symbol '+''.join(modsym),None
824    if not checkModSym():
825        return 'Modulation '+''.join(modsym)+' not consistent with space group '+SGData['SpGrp'],None
826    modQ = [Fracs[mod] for mod in modsym]
827    SSGKl = SGData['SSGKl'][:]
828    if SGData['SGLaue'] in ['2/m','mmm']:
829        SSGKl = fixMonoOrtho()
830    if len(gensym) and len(gensym) != len(SSGKl):
831        return 'Wrong number of items in generator symbol '+''.join(gensym),None
832    if not checkGen(gensym):
833        return 'Generator '+''.join(gensym)+' not consistent with space group '+SGData['SpGrp'],None
834    gensym = specialGen(gensym)
835    genQ = [Fracs[mod] for mod in gensym]
836    if not genQ:
837        genQ = [0,0,0,0]
838    SSGData = {'SSpGrp':SGData['SpGrp']+SSymbol,'modQ':modQ,'modSymb':modsym}
839    SSCen = np.zeros((len(SGData['SGCen']),4))
840    for icen,cen in enumerate(SGData['SGCen']):
841        SSCen[icen,0:3] = cen
842    SSCen[0] = np.zeros(4)
843    SSGData['SSGCen'] = SSCen
844    SSGData['SSGOps'] = []
845    for iop,op in enumerate(SGData['SGOps']):
846        T = np.zeros(4)
847        ssop = np.zeros((4,4))
848        ssop[:3,:3] = op[0]
849        T[:3] = op[1]
850        SSGData['SSGOps'].append([ssop,T])
851    E,Result = genSSGOps()
852    if E:
853        SSGData['SSGOps'] = Result                     
854        return None,SSGData
855    else:
856        return Result+'\nOperator conflict - incorrect superspace symbol',None
857
858def SSGPrint(SGData,SSGData):
859    '''
860    Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
861
862    :param SGData: space group data structure as defined in SpcGroup above.
863    :param SSGData: from :func:`SSpcGroup`
864    :returns:
865        SSGText - list of strings with the superspace group details
866        SGTable - list of strings for each of the operations
867    '''
868    Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])
869    SSGText = []
870    SSGText.append(' Superspace Group: '+SSGData['SSpGrp'])
871    CentStr = 'centrosymmetric'
872    if not SGData['SGInv']:
873        CentStr = 'non'+CentStr
874    if SGData['SGLatt'] in 'ABCIFR':
875        SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
876    else:
877        SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())       
878    SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
879    SSGText.append(' The superlattice point group is '+SGData['SGPtGrp']+','+''.join([str(i) for i in SGData['SSGKl']]))
880    SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
881    SSGText.append(' Multiplicity of a general site is '+str(Mult))
882    if SGData['SGUniq'] in ['a','b','c']:
883        SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
884    if SGData['SGInv']:
885        SSGText.append(' The inversion center is located at 0,0,0')
886    if SGData['SGPolax']:
887        SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
888    SSGText.append(' ')
889    if len(SSGData['SSGCen']) > 1:
890        SSGText.append(' The equivalent positions are:')
891        SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
892    else:
893        SSGText.append(' The equivalent positions are:\n')
894    SSGTable = []
895    for i,Opr in enumerate(SSGData['SSGOps']):
896        SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
897    return SSGText,SSGTable
898   
899def SSGModCheck(Vec,SSGData):
900    ''' Checks modulation vector compatibility with supersymmetry space group symbol.
901    Superspace group symbol takes precidence & the vector will be modified accordingly
902    '''
903    modQ = SSGData['modQ']
904    modSymb = SSGData['modSymb']
905    Vec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
906    return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(Vec,modSymb,modQ)]
907
908def SSMT2text(Opr):
909    "From superspace group matrix/translation operator returns text version"
910    XYZS = ('x','y','z','t')    #Stokes, Campbell & van Smaalen notation
911    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
912    Fld = ''
913    M,T = Opr
914    for j in range(4):
915        IJ = ''
916        for k in range(4):
917            txt = str(int(round(M[j][k])))
918            txt = txt.replace('1',XYZS[k]).replace('0','')
919            if '2' in txt:
920                txt += XYZS[k]
921            if IJ and M[j][k] > 0:
922                IJ += '+'+txt
923            else:
924                IJ += txt
925        IK = int(round(T[j]*12))%12
926        if IK:
927            if not IJ:
928                break
929            if IJ[0] == '-':
930                Fld += (TRA[IK]+IJ).rjust(8)
931            else:
932                Fld += (TRA[IK]+'+'+IJ).rjust(8)
933        else:
934            Fld += IJ.rjust(8)
935        if j != 3: Fld += ', '
936    return Fld
937   
938def SSLatt2text(SSGCen):
939    "Lattice centering vectors to text"
940    lattTxt = ''
941    for vec in SSGCen:
942        lattTxt += ' '
943        for item in vec:
944            if int(item*12.):
945                lattTxt += '1/%d,'%(12/int(item*12))
946            else:
947                lattTxt += '0,'
948        lattTxt = lattTxt.rstrip(',')
949        lattTxt += ';'
950    lattTxt = lattTxt.rstrip(';').lstrip(' ')
951    return lattTxt
952       
953def SSpaceGroup(SGSymbol,SSymbol):
954    '''
955    Print the output of SSpcGroup in a nicely formatted way.
956
957    :param SGSymbol: space group symbol with spaces between axial fields.
958    :param SSymbol: superspace group symbol extension (string).
959    :returns: nothing
960    '''
961
962    E,A = SpcGroup(SGSymbol)
963    if E > 0:
964        print SGErrors(E)
965        return
966    E,B = SSpcGroup(A,SSymbol)   
967    if E > 0:
968        print E
969        return
970    for l in SSGPrint(B):
971        print l
972       
973def SGProd(OpA,OpB):
974    '''
975    Form space group operator product. OpA & OpB are [M,V] pairs;
976        both must be of same dimension (3 or 4). Returns [M,V] pair
977    '''
978    A,U = OpA
979    B,V = OpB
980    M = np.inner(B.T,A)
981    W = np.inner(B,U)+V
982    return M.T,W
983       
984def MoveToUnitCell(xyz):
985    '''
986    Translates a set of coordinates so that all values are >=0 and < 1
987
988    :param xyz: a list or numpy array of fractional coordinates
989    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
990    '''
991    XYZ = np.zeros(3)
992    for i,x in enumerate(xyz):
993        XYZ[i] = (x-int(x))%1.0
994    return XYZ
995       
996def Opposite(XYZ,toler=0.0002):
997    '''
998    Gives opposite corner, edge or face of unit cell for position within tolerance.
999        Result may be just outside the cell within tolerance
1000
1001    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
1002    :param toler: unit cell fraction tolerance making opposite
1003    :returns:
1004        XYZ: array of opposite positions; always contains XYZ
1005    '''
1006    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]]
1007    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
1008    perm = TB*perm3
1009    cperm = ['%d%d%d'%(i,j,k) for i,j,k in perm]
1010    D = dict(zip(cperm,perm))
1011    new = []
1012    for key in D:
1013        new.append(np.array(D[key])+np.array(XYZ))
1014    return new
1015       
1016def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
1017    '''
1018    Generates the equivalent positions for a specified coordinate and space group
1019
1020    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1021    :param SGData: from :func:`SpcGroup`
1022    :param All: True return all equivalent positions including duplicates;
1023      False return only unique positions
1024    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
1025    :param Move: True move generated atom positions to be inside cell
1026      False do not move atoms       
1027    :return: [[XYZEquiv],Idup,[UijEquiv]]
1028
1029      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
1030      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
1031      * is centering operator number (1-4) and - is for inversion
1032        Cell = unit cell translations needed to put new positions inside cell
1033        [UijEquiv] - equivalent Uij; absent if no Uij given
1034       
1035    '''
1036    XYZEquiv = []
1037    UijEquiv = []
1038    Idup = []
1039    Cell = []
1040    X = np.array(XYZ)
1041    if Move:
1042        X = MoveToUnitCell(X)
1043    for ic,cen in enumerate(SGData['SGCen']):
1044        C = np.array(cen)
1045        for invers in range(int(SGData['SGInv']+1)):
1046            for io,[M,T] in enumerate(SGData['SGOps']):
1047                idup = ((io+1)+100*ic)*(1-2*invers)
1048                XT = np.inner(M,X)+T
1049                if len(Uij):
1050                    U = Uij2U(Uij)
1051                    U = np.inner(M,np.inner(U,M).T)
1052                    newUij = U2Uij(U)
1053                if invers:
1054                    XT = -XT
1055                XT += C
1056                if Move:
1057                    newX = MoveToUnitCell(XT)
1058                else:
1059                    newX = XT
1060                cell = np.asarray(np.rint(newX-XT),dtype=np.int32)
1061                if All:
1062                    if np.allclose(newX,X,atol=0.0002):
1063                        idup = False
1064                else:
1065                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1066                        idup = False
1067                if All or idup:
1068                    XYZEquiv.append(newX)
1069                    Idup.append(idup)
1070                    Cell.append(cell)
1071                    if len(Uij):
1072                        UijEquiv.append(newUij)                   
1073    if len(Uij):
1074        return zip(XYZEquiv,UijEquiv,Idup,Cell)
1075    else:
1076        return zip(XYZEquiv,Idup,Cell)
1077
1078def GenHKLf(HKL,SGData):
1079    '''
1080    Uses old GSAS Fortran routine genhkl.for
1081
1082    :param HKL:  [h,k,l]
1083    :param SGData: space group data obtained from SpcGroup
1084    :returns: iabsnt,mulp,Uniq,phi
1085
1086     *   iabsnt = True if reflection is forbidden by symmetry
1087     *   mulp = reflection multiplicity including Friedel pairs
1088     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
1089
1090    '''
1091    hklf = HKL+[0,]
1092    Ops = SGData['SGOps']
1093    OpM = np.array([op[0] for op in Ops])
1094    OpT = np.array([op[1] for op in Ops])
1095    Inv = SGData['SGInv']
1096    Cen = np.array([cen for cen in SGData['SGCen']])
1097   
1098    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1099    h,k,l,f = Uniq
1100    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
1101    phi = f[:Nuniq]
1102   
1103    return iabsnt,mulp,Uniq,phi
1104                                 
1105def GetOprPtrName(key):
1106    'Needs a doc string'
1107    OprPtrName = {
1108        '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
1109        '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
1110        '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
1111        '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
1112        '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
1113        '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
1114        '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
1115        '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
1116        '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
1117        '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
1118        '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
1119        '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
1120        '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
1121        '2223' :[  68,' 3+++2',39],
1122        '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
1123        '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
1124        '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
1125        '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
1126        '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
1127        '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
1128        '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
1129        '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
1130        '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
1131        }
1132    return OprPtrName[key]
1133
1134def GetKNsym(key):
1135    'Needs a doc string'
1136    KNsym = {
1137        '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
1138        '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
1139        '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
1140        '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
1141        '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
1142        '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
1143        '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
1144        '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
1145        '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
1146        '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
1147        '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
1148        '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(100)','16777280'  :' -4(100)',
1149        '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
1150        '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
1151        '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
1152        '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
1153        '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
1154        '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
1155        '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
1156        '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
1157        '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
1158        '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
1159        '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
1160        '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
1161        '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
1162        '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
1163        '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
1164        '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
1165        '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
1166        '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
1167        '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
1168        '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
1169        '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
1170        '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
1171        '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
1172        '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
1173        '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
1174        '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   '
1175        }
1176    return KNsym[key]       
1177
1178def GetNXUPQsym(siteSym):       
1179    'Needs a doc string'
1180    NXUPQsym = {
1181        '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'    2(x)':(12,18,12,25),'    m(x)':(25,18,12,25),
1182        '  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),
1183        '    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),
1184        '   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),
1185        '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),
1186        '  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),
1187        '   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),
1188        '2/m(+-0)':( 1,20, 0,-1),'  mm2(x)':(12,10, 0,-1),'  mm2(y)':(13,10, 0,-1),'  mm2(z)':(14,10, 0,-1),
1189        ' mm2(yz)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),' mm2(xz)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
1190        ' mm2(xy)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'  222(x)':( 1,13, 0,-1),
1191        '  222(y)':( 1,12, 0,-1),'  222(z)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'  mmm(x)':( 1,13, 0,-1),
1192        '  mmm(y)':( 1,12, 0,-1),'  mmm(z)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
1193        '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),
1194        '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),
1195        '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
1196        '  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),
1197        '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
1198        ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
1199        '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
1200        '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
1201        ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
1202        ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
1203        '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
1204        ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
1205        '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
1206        '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
1207        '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
1208        ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
1209        '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
1210        'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
1211        'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
1212        'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
1213        '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),
1214        '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),
1215        ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
1216        ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
1217        ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
1218        ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
1219        }
1220    return NXUPQsym[siteSym]
1221
1222def GetCSxinel(siteSym): 
1223    'Needs a doc string'
1224    CSxinel = [[],                         # 0th empty - indices are Fortran style
1225        [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
1226        [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
1227        [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
1228        [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
1229        [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
1230        [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
1231        [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
1232        [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
1233        [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
1234        [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
1235        [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
1236        [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
1237        [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
1238        [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
1239        [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
1240        [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
1241        [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
1242        [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
1243        [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
1244        [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
1245        [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
1246        [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
1247        [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
1248        [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
1249        [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
1250        [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
1251        [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
1252        [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
1253        ]
1254    indx = GetNXUPQsym(siteSym)
1255    return CSxinel[indx[0]]
1256   
1257def GetCSuinel(siteSym):
1258    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
1259    CSuinel = [[],                                             # 0th empty - indices are Fortran style
1260        [[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
1261        [[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
1262        [[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
1263        [[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
1264        [[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
1265        [[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
1266        [[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
1267        [[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
1268        [[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
1269        [[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
1270        [[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
1271        [[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
1272        [[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
1273        [[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
1274        [[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
1275        [[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
1276        [[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
1277        [[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
1278        [[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
1279        [[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
1280        [[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
1281        [[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
1282        [[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
1283        [[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
1284        [[1,2,3,2,4,4],[ 1.0, 1.0, 1.0, 0.5, 0.5, 1.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
1285        [[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
1286        [[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
1287        [[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
1288        [[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
1289        ]
1290    indx = GetNXUPQsym(siteSym)
1291    return CSuinel[indx[1]]
1292   
1293def MustrainNames(SGData):
1294    'Needs a doc string'
1295    laue = SGData['SGLaue']
1296    uniq = SGData['SGUniq']
1297    if laue in ['m3','m3m']:
1298        return ['S400','S220']
1299    elif laue in ['6/m','6/mmm','3m1']:
1300        return ['S400','S004','S202']
1301    elif laue in ['31m','3']:
1302        return ['S400','S004','S202','S211']
1303    elif laue in ['3R','3mR']:
1304        return ['S400','S220','S310','S211']
1305    elif laue in ['4/m','4/mmm']:
1306        return ['S400','S004','S220','S022']
1307    elif laue in ['mmm']:
1308        return ['S400','S040','S004','S220','S202','S022']
1309    elif laue in ['2/m']:
1310        SHKL = ['S400','S040','S004','S220','S202','S022']
1311        if uniq == 'a':
1312            SHKL += ['S013','S031','S211']
1313        elif uniq == 'b':
1314            SHKL += ['S301','S103','S121']
1315        elif uniq == 'c':
1316            SHKL += ['S130','S310','S112']
1317        return SHKL
1318    else:
1319        SHKL = ['S400','S040','S004','S220','S202','S022']
1320        SHKL += ['S310','S103','S031','S130','S301','S013']
1321        SHKL += ['S211','S121','S112']
1322        return SHKL
1323
1324def HStrainNames(SGData):
1325    'Needs a doc string'
1326    laue = SGData['SGLaue']
1327    uniq = SGData['SGUniq']
1328    if laue in ['m3','m3m']:
1329        return ['D11','eA']         #add cubic strain term
1330    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1331        return ['D11','D33']
1332    elif laue in ['3R','3mR']:
1333        return ['D11','D12']
1334    elif laue in ['4/m','4/mmm']:
1335        return ['D11','D33']
1336    elif laue in ['mmm']:
1337        return ['D11','D22','D33']
1338    elif laue in ['2/m']:
1339        Dij = ['D11','D22','D33']
1340        if uniq == 'a':
1341            Dij += ['D23']
1342        elif uniq == 'b':
1343            Dij += ['D13']
1344        elif uniq == 'c':
1345            Dij += ['D12']
1346        return Dij
1347    else:
1348        Dij = ['D11','D22','D33','D12','D13','D23']
1349        return Dij
1350   
1351def MustrainCoeff(HKL,SGData):
1352    'Needs a doc string'
1353    #NB: order of terms is the same as returned by MustrainNames
1354    laue = SGData['SGLaue']
1355    uniq = SGData['SGUniq']
1356    h,k,l = HKL
1357    Strm = []
1358    if laue in ['m3','m3m']:
1359        Strm.append(h**4+k**4+l**4)
1360        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
1361    elif laue in ['6/m','6/mmm','3m1']:
1362        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
1363        Strm.append(l**4)
1364        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
1365    elif laue in ['31m','3']:
1366        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
1367        Strm.append(l**4)
1368        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
1369        Strm.append(4.0*h*k*l*(h+k))
1370    elif laue in ['3R','3mR']:
1371        Strm.append(h**4+k**4+l**4)
1372        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
1373        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
1374        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
1375    elif laue in ['4/m','4/mmm']:
1376        Strm.append(h**4+k**4)
1377        Strm.append(l**4)
1378        Strm.append(3.0*(h*k)**2)
1379        Strm.append(3.0*((h*l)**2+(k*l)**2))
1380    elif laue in ['mmm']:
1381        Strm.append(h**4)
1382        Strm.append(k**4)
1383        Strm.append(l**4)
1384        Strm.append(3.0*(h*k)**2)
1385        Strm.append(3.0*(h*l)**2)
1386        Strm.append(3.0*(k*l)**2)
1387    elif laue in ['2/m']:
1388        Strm.append(h**4)
1389        Strm.append(k**4)
1390        Strm.append(l**4)
1391        Strm.append(3.0*(h*k)**2)
1392        Strm.append(3.0*(h*l)**2)
1393        Strm.append(3.0*(k*l)**2)
1394        if uniq == 'a':
1395            Strm.append(2.0*k*l**3)
1396            Strm.append(2.0*l*k**3)
1397            Strm.append(4.0*k*l*h**2)
1398        elif uniq == 'b':
1399            Strm.append(2.0*l*h**3)
1400            Strm.append(2.0*h*l**3)
1401            Strm.append(4.0*h*l*k**2)
1402        elif uniq == 'c':
1403            Strm.append(2.0*h*k**3)
1404            Strm.append(2.0*k*h**3)
1405            Strm.append(4.0*h*k*l**2)
1406    else:
1407        Strm.append(h**4)
1408        Strm.append(k**4)
1409        Strm.append(l**4)
1410        Strm.append(3.0*(h*k)**2)
1411        Strm.append(3.0*(h*l)**2)
1412        Strm.append(3.0*(k*l)**2)
1413        Strm.append(2.0*k*h**3)
1414        Strm.append(2.0*h*l**3)
1415        Strm.append(2.0*l*k**3)
1416        Strm.append(2.0*h*k**3)
1417        Strm.append(2.0*l*h**3)
1418        Strm.append(2.0*k*l**3)
1419        Strm.append(4.0*k*l*h**2)
1420        Strm.append(4.0*h*l*k**2)
1421        Strm.append(4.0*k*h*l**2)
1422    return Strm
1423   
1424def Muiso2Shkl(muiso,SGData,cell):
1425    "this is to convert isotropic mustrain to generalized Shkls"
1426    import GSASIIlattice as G2lat
1427    A = G2lat.cell2AB(cell)[0]
1428   
1429    def minMus(Shkl,muiso,H,SGData,A):
1430        U = np.inner(A.T,H)
1431        S = np.array(MustrainCoeff(U,SGData))
1432        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
1433        rad = np.sqrt(np.sum((Sum[:,np.newaxis]*H)**2,axis=1))
1434        return (muiso-rad)**2
1435       
1436    laue = SGData['SGLaue']
1437    PHI = np.linspace(0.,360.,60,True)
1438    PSI = np.linspace(0.,180.,60,True)
1439    X = np.outer(npsind(PHI),npsind(PSI))
1440    Y = np.outer(npcosd(PHI),npsind(PSI))
1441    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1442    HKL = np.dstack((X,Y,Z))
1443    if laue in ['m3','m3m']:
1444        S0 = [1000.,1000.]
1445    elif laue in ['6/m','6/mmm','3m1']:
1446        S0 = [1000.,1000.,1000.]
1447    elif laue in ['31m','3']:
1448        S0 = [1000.,1000.,1000.,1000.]
1449    elif laue in ['3R','3mR']:
1450        S0 = [1000.,1000.,1000.,1000.]
1451    elif laue in ['4/m','4/mmm']:
1452        S0 = [1000.,1000.,1000.,1000.]
1453    elif laue in ['mmm']:
1454        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
1455    elif laue in ['2/m']:
1456        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
1457    else:
1458        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
1459            1000.,1000.,0.,0.,0.]
1460    S0 = np.array(S0)
1461    HKL = np.reshape(HKL,(-1,3))
1462    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
1463    return result[0]
1464       
1465def SytSym(XYZ,SGData):
1466    '''
1467    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
1468
1469    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1470    :param SGData: from SpcGroup
1471    :Returns: a two element tuple:
1472
1473     * The 1st element is a code for the site symmetry (see GetKNsym)
1474     * The 2nd element is the site multiplicity
1475
1476    '''
1477    def PackRot(SGOps):
1478        IRT = []
1479        for ops in SGOps:
1480            M = ops[0]
1481            irt = 0
1482            for j in range(2,-1,-1):
1483                for k in range(2,-1,-1):
1484                    irt *= 3
1485                    irt += M[k][j]
1486            IRT.append(int(irt))
1487        return IRT
1488       
1489    SymName = ''
1490    Mult = 1
1491    Isym = 0
1492    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1493        Isym = 1073741824
1494    Jdup = 0
1495    Xeqv = GenAtom(XYZ,SGData,True)
1496    IRT = PackRot(SGData['SGOps'])
1497    L = -1
1498    for ic,cen in enumerate(SGData['SGCen']):
1499        for invers in range(int(SGData['SGInv']+1)):
1500            for io,ops in enumerate(SGData['SGOps']):
1501                irtx = (1-2*invers)*IRT[io]
1502                L += 1
1503                if not Xeqv[L][1]:
1504                    Jdup += 1
1505                    jx = GetOprPtrName(str(irtx))
1506                    if jx[2] < 39:
1507                        Isym += 2**(jx[2]-1)
1508    if Isym == 1073741824: Isym = 0
1509    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)/Jdup
1510         
1511    return GetKNsym(str(Isym)),Mult
1512   
1513def ElemPosition(SGData):
1514    ''' Under development.
1515    Object here is to return a list of symmetry element types and locations suitable
1516    for say drawing them.
1517    So far I have the element type... getting all possible locations without lookup may be impossible!
1518    '''
1519    SymElements = []
1520    Inv = SGData['SGInv']
1521    Cen = SGData['SGCen']
1522    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
1523    # get operators & expand if centrosymmetric
1524    Ops = SGData['SGOps']
1525    opM = np.array([op[0].T for op in Ops])
1526    opT = np.array([op[1] for op in Ops])
1527    if Inv:
1528        opM = np.concatenate((opM,-opM))
1529        opT = np.concatenate((opT,-opT))
1530    opMT = zip(opM,opT)
1531    for M,T in opMT[1:]:        #skip I
1532        Dt = int(nl.det(M))
1533        Tr = int(np.trace(M))
1534        Dt = -(Dt-1)/2
1535        Es = eleSym[Tr][Dt]
1536        if Dt:              #rotation-inversion
1537            I = np.eye(3)
1538            if Tr == 1:     #mirrors/glides
1539                if np.any(T):       #glide
1540                    M2 = np.inner(M,M)
1541                    MT = np.inner(M,T)+T
1542                    print 'glide',Es,MT
1543                    print M2
1544                else:               #mirror
1545                    print 'mirror',Es,T
1546                    print I-M
1547                X = [-1,-1,-1]
1548            elif Tr == -3:  # pure inversion
1549                X = np.inner(nl.inv(I-M),T)
1550                print 'inversion',Es,X
1551            else:           #other rotation-inversion
1552                M2 = np.inner(M,M)
1553                MT = np.inner(M,T)+T
1554                print 'rot-inv',Es,MT
1555                print M2
1556                X = [-1,-1,-1]
1557        else:               #rotations
1558            print 'rotation',Es
1559            X = [-1,-1,-1]
1560        #SymElements.append([Es,X])
1561       
1562    return #SymElements
1563   
1564def ApplyStringOps(A,SGData,X,Uij=[]):
1565    'Needs a doc string'
1566    SGOps = SGData['SGOps']
1567    SGCen = SGData['SGCen']
1568    Ax = A.split('+')
1569    Ax[0] = int(Ax[0])
1570    iC = 0
1571    if Ax[0] < 0:
1572        iC = 1
1573    Ax[0] = abs(Ax[0])
1574    nA = Ax[0]%100-1
1575    cA = Ax[0]/100
1576    Cen = SGCen[cA]
1577    M,T = SGOps[nA]
1578    if len(Ax)>1:
1579        cellA = Ax[1].split(',')
1580        cellA = np.array([int(a) for a in cellA])
1581    else:
1582        cellA = np.zeros(3)
1583    newX = (1-2*iC)*(Cen+np.inner(M,X)+T)+cellA
1584    if len(Uij):
1585        U = Uij2U(Uij)
1586        U = np.inner(M,np.inner(U,M).T)
1587        newUij = U2Uij(U)
1588        return [newX,newUij]
1589    else:
1590        return newX
1591       
1592def StringOpsProd(A,B,SGData):
1593    """
1594    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
1595    where '-' indicates inversion, c(>0) is the cell centering operator,
1596    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
1597    Should return resultant string - C. SGData - dictionary using entries:
1598
1599       *  'SGCen': cell centering vectors [0,0,0] at least
1600       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
1601
1602    """
1603    SGOps = SGData['SGOps']
1604    SGCen = SGData['SGCen']
1605    #1st split out the cell translation part & work on the operator parts
1606    Ax = A.split('+'); Bx = B.split('+')
1607    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
1608    iC = 0
1609    if Ax[0]*Bx[0] < 0:
1610        iC = 1
1611    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
1612    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
1613    cA = Ax[0]/100;  cB = Bx[0]/100
1614    Cen = (SGCen[cA]+SGCen[cB])%1.0
1615    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
1616    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
1617    Mc = np.inner(Ma,Mb.T)
1618#    print Ma,Mb,Mc
1619    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
1620#    print Ta,Tb,Tc
1621#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
1622    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
1623    #now the cell translation part
1624    if len(Ax)>1:
1625        cellA = Ax[1].split(',')
1626        cellA = [int(a) for a in cellA]
1627    else:
1628        cellA = [0,0,0]
1629    if len(Bx)>1:
1630        cellB = Bx[1].split(',')
1631        cellB = [int(b) for b in cellB]
1632    else:
1633        cellB = [0,0,0]
1634    cellC = np.add(cellA,cellB)
1635    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
1636        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
1637    return C
1638           
1639def U2Uij(U):
1640    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
1641    return [U[0][0],U[1][1],U[2][2],2.*U[0][1],2.*U[0][2],2.*U[1][2]]
1642   
1643def Uij2U(Uij):
1644    #returns the thermal motion tensor U from Uij as numpy array
1645    return np.array([[Uij[0],Uij[3]/2.,Uij[4]/2.],[Uij[3]/2.,Uij[1],Uij[5]/2.],[Uij[4]/2.,Uij[5]/2.,Uij[2]]])
1646
1647def StandardizeSpcName(spcgroup):
1648    '''Accept a spacegroup name where spaces may have not been used
1649    in the names according to the GSAS convention (spaces between symmetry
1650    for each axis) and return the space group name as used in GSAS
1651    '''
1652    rspc = spcgroup.replace(' ','').upper()
1653    # deal with rhombohedral and hexagonal setting designations
1654    rhomb = ''
1655    if rspc[-1:] == 'R':
1656        rspc = rspc[:-1]
1657        rhomb = ' R'
1658    if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
1659        rspc = rspc[:-1]
1660    # look for a match in the spacegroup lists
1661    for i in spglist.values():
1662        for spc in i:
1663            if rspc == spc.replace(' ','').upper():
1664                return spc + rhomb
1665    # how about the post-2002 orthorhombic names?
1666    for i,spc in sgequiv_2002_orthorhombic:
1667        if rspc == i.replace(' ','').upper():
1668            return spc
1669    # not found
1670    return ''
1671
1672   
1673spglist = {}
1674'''A dictionary of space groups as ordered and named in the pre-2002 International
1675Tables Volume A, except that spaces are used following the GSAS convention to
1676separate the different crystallographic directions.
1677Note that the symmetry codes here will recognize many non-standard space group
1678symbols with different settings. They are ordered by Laue group
1679'''
1680spglist = {
1681    'P1' : ('P 1','P -1',), # 1-2
1682    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
1683        '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
1684    'C2/m':('C 2','C m','C c','C n',
1685        'C 2/m','C 2/c','C 2/n',),
1686    'Pmmm':('P 2 2 2',
1687        'P 2 2 21','P 2 21 2','P 21 2 2',
1688        'P 21 21 2','P 21 2 21','P 2 21 21',
1689        'P 21 21 21',
1690        'P m m 2','P m 2 m','P 2 m m',
1691        'P m c 21','P c m 21','P 21 m a','P 21 a m','P b 21 m','P m 21 b',
1692        'P c c 2','P 2 a a','P b 2 b',
1693        'P m a 2','P b m 2','P 2 m b','P 2 c m','P c 2 m','P m 2 a',
1694        'P c a 21','P b c 21','P 21 a b','P 21 c a','P c 21 b','P b 21 a',
1695        'P n c 2','P c n 2','P 2 n a','P 2 a n','P b 2 n','P n 2 b',
1696        'P m n 21','P n m 21','P 21 m n','P 21 n m','P n 21 m','P m 21 n',
1697        'P b a 2','P 2 c b','P c 2 a',
1698        'P n a 21','P b n 21','P 21 n b','P 21 c n','P c 21 n','P n 21 a',
1699        'P n n 2','P 2 n n','P n 2 n',
1700        'P m m m','P n n n',
1701        'P c c m','P m a a','P b m b',
1702        'P b a n','P n c b','P c n a',
1703        'P m m a','P m m b','P b m m','P c m m','P m c m','P m a m',
1704        'P n n a','P n n b','P b n n','P c n n','P n c n','P n a n',
1705        'P m n a','P n m b','P b m n','P c n m','P n c m','P m a n',
1706        'P c c a','P c c b','P b a a','P c a a','P b c b','P b a b',
1707        'P b a m','P m c b','P c m a',
1708        'P c c n','P n a a','P b n b',
1709        'P b c m','P c a m','P m c a','P m a b','P b m a','P c m b',
1710        'P n n m','P m n n','P n m n',
1711        'P m m n','P n m m','P m n m',
1712        'P b c n','P c a n','P n c a','P n a b','P b n a','P c n b',
1713        'P b c a','P c a b',
1714        'P n m a','P m n b','P b n m','P c m n','P m c n','P n a m',
1715        ),
1716    'Cmmm':('C 2 2 21','C 2 2 2','C m m 2','C m c 21','C c c 2','C m 2 m','C 2 m m',
1717        'C m 2 a','C 2 m b','C 2 c m','C c 2 m','C 2 c m',
1718        'C m c a','C m m m','C c c m','C m m a','C c c a','C m c m',),
1719    'Immm':('I 2 2 2','I 21 21 21','I m m m',
1720        'I m m 2','I m 2 m','I 2 m m',
1721        'I b a 2','I 2 c b','I c 2 a',
1722        'I m a 2','I b m 2','I 2 m b','I 2 c m','I c 2 m','I m 2 a',
1723        'I b a m','I m c b','I c m a',
1724        'I b c a','I c a b',
1725        'I m m a','I m m b','I b m m ','I c m m','I m c m','I m a m',),
1726    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
1727        'F m m 2','F m 2 m','F 2 m m',
1728        'F d d 2','F d 2 d','F 2 d d',),
1729    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
1730        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
1731        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
1732        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
1733        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
1734        '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',
1735        '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',
1736        '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',
1737        'P 42/n c m',),
1738    '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',
1739        'I 4 c m','I 41 m d','I 41 c d',
1740        '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',
1741        'I 41/a m d','I 41/a c d'),
1742    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
1743    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
1744        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
1745        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
1746        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
1747        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
1748        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
1749        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
1750    'Pm3m': ('P 2 3','P 21 3','P m 3','P n 3','P a 3','P 4 3 2','P 42 3 2',
1751        'P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n','P m 3 m','P n 3 n',
1752        'P m 3 n','P n 3 m',),
1753    'Im3m':('I 2 3','I 21 3','I m -3','I a -3', 'I 4 3 2','I 41 3 2',
1754        'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a -3 d',),
1755    'Fm3m':('F 2 3','F m -3','F d -3','F 4 3 2','F 41 3 2','F -4 3 m',
1756        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
1757}
1758
1759ssdict = {}
1760'''A dictionary of superspace group symbols allowed for each entry in spglist
1761(except cubics). Monoclinics are all b-unique setting.
1762'''
1763ssdict = {
1764    'P 1':['(abg)',],'P -1':['(abg)',],
1765    #monoclinic - done
1766    'P 2':['(a0g)','(a1/2g)','(0b0)','(0b0)s','(1/2b0)','(0b1/2)',],
1767    'P 21':['(a0g)','(0b0)','(0b0)s','(1/2b0)','(0b1/2)','(1/2b0)s','(0b1/2)s',],
1768    'P m':['(a0g)','(a0g)s','(a1/2g)','(a1/2g)s','(0b0)','(1/2b0)','(0b1/2)',],
1769    'P a':['(a0g)','(a1/2g)','(a0g)s','(a1/2g)s','(0b0)','(0b1/2)',],
1770    'P c':['(a0g)','(a1/2g)','(a0g)s','(a1/2g)s','(0b0)','(1/2b0)',],
1771    'P n':['(a0g)','(a1/2g)','(a0g)s','(a1/2g)s','(0b0)','(1/2b1/2)',],
1772    'P 2/m':['(a0g)','(a1/2g)','(a0g)0s','(a1/2g)0s',
1773        '(0b0)','(0b0)s0','(1/2b0)','(0b1/2)','(1/2b0)s0','(0b1/2)s0',],
1774    'P 21/m':['(a0g)','(a0g)0s','(0b0)','(0b0)s',
1775        '(1/2b0)','(0b1/2)','(1/2b0)s0','(0b1/2)s0'],
1776    'P 2/c':['(a0g)','(a1/2g)','(a0g)0s','(a1/2g)0s',
1777        '(0b0)','(0b0)s0','(1/2b0)','(1/2b0)s0',],
1778    'P 2/a':['(a0g)','(a1/2g)','(a0g)0s','(a1/2g)0s',
1779        '(0b0)','(0b0)s0','(0b1/2)','(0b1/2)s0',],
1780    'P 2/n':['(a0g)','(a1/2g)','(a0g)0s','(a1/2g)0s',
1781        '(0b0)','(0b0)s0','(1/2b1/2)','(1/2b1/2)s0',],
1782    'P 21/c':['(a0g)','(a0g)0s','(0b0)','(0b0)s0','(1/2b0)','(1/2b0)s0',],
1783    'P 21/a':['(a0g)','(a0g)0s','(0b0)','(0b0)s0','(0b1/2)','(0b1/2)s0',],
1784    'P 21/n':['(a0g)','(a0g)0s','(0b0)','(0b0)s0','(1/2b1/2)','(1/2b1/2)s0',],
1785    'C 2':['(a0g)','(0b0)','(0b0)s','(0b1/2)','(0b1/2)s',],
1786    'C m':['(a0g)','(a0g)s','(0b0)','(0b1/2)',],
1787    'C c':['(a0g)','(a0g)s','(0b0)',],
1788    'C n':['(a0g)','(a0g)s','(0b0)',],
1789    'C 2/m':['(a0g)','(a0g)0s','(0b0)','(0b0)s0',
1790        '(1/2b0)','(0b1/2)','(1/2b0)s0','(0b1/2)s0',],
1791    'C 2/c':['(a0g)','(a0g)0s','(0b0)','(0b0)s0',],
1792    'C 2/n':['(a0g)','(a0g)0s','(0b0)','(0b0)s0',],
1793    #orthorhombic   
1794    'P 2 2 2':['(a00)','(0b0)','(00g)','(a00)s00','(0b0)0s0','(00g)00s',
1795        '(a1/20)','(a01/2)','(0b1/2)','(1/2b0)','(01/2g)','(1/20g)',
1796        '(a1/21/2)','(1/2b1/2)','(1/21/2g)',],
1797       
1798    'P 2 2 21':['(a00)','(0b0)','(00g)','(a00)s00','(0b0)0s0',
1799        '(a1/20)','(1/2b0)','(01/2g)','(1/20g)','(1/21/2g)',],
1800    'P 2 21 2':['(a00)','(0b0)','(00g)','(a00)s00','(00g)00s',
1801        '(a01/2)','(1/20g)','(0b1/2)','(1/2b0)','(1/2b1/2)',],
1802    'P 21 2 2':['(a00)','(0b0)','(00g)','(0b0)0s0','(00g)00s',
1803        '(01/2g)','(0b1/2)','(a01/2)','(a1/20)','(a1/21/2)',],
1804       
1805    'P 21 21 2':['(a00)','(0b0)','(00g)','(00g)00s','(a01/2)','(0b1/2)',],
1806    'P 21 2 21':['(a00)','(0b0)','(00g)','(0b0)0s0','(a1/20)','(01/2g)',],
1807    'P 2 21 21':['(a00)','(0b0)','(00g)','(a00)s00','(1/2b0)','(1/20g)',],
1808       
1809    'P 21 21 21':['(a00)','(0b0)','(00g)',],
1810       
1811    'P m m 2':['(a00)','(0b0)','(00g)',
1812        '(a00)s00','(0b0)0s0','(00g)s0s','(00g)ss0','(00g)0ss',
1813        '(a1/20)','(a01/2)','(0b1/2)','(1/2b0)','(01/2g)','(1/20g)',
1814        '(a1/21/2)','(1/2b1/2)','(1/21/2g)',],
1815    'P m 2 m':[],
1816    'P 2 m m':[],
1817       
1818    'P m c 21':[],
1819    'P c m 21':[],
1820    'P 21 m a':[],
1821    'P 21 a m':[],
1822    'P b 21 m':[],
1823    'P m 21 b':[],
1824       
1825    'P c c 2':[],
1826    'P 2 a a':[],
1827    'P b 2 b':[],
1828       
1829    'P m a 2':[],
1830    'P b m 2':[],
1831    'P 2 m b':[],
1832    'P 2 c m':[],
1833    'P c 2 m':[],
1834    'P m 2 a':[],
1835       
1836    'P c a 21':[],
1837    'P b c 21':[],
1838    'P 21 a b':[],
1839    'P 21 c a':[],
1840    'P c 21 b':[],
1841    'P b 21 a':[],
1842       
1843    'P n c 2':[],
1844    'P c n 2':[],
1845    'P 2 n a':[],
1846    'P 2 a n':[],
1847    'P b 2 n':[],
1848    'P n 2 b':[],
1849       
1850    'P m n 21':[],
1851    'P n m 21':[],
1852    'P 21 m n':[],
1853    'P 21 n m':[],
1854    'P n 21 m':[],
1855    'P m 21 n':[],
1856       
1857    'P b a 2':[],
1858    'P 2 c b':[],
1859    'P c 2 a':[],
1860       
1861    'P n a 21':[],
1862    'P b n 21':[],
1863    'P 21 n b':[],
1864    'P 21 c n':[],
1865    'P c 21 n':[],
1866    'P n 21 a':[],
1867       
1868    'P n n 2':[],
1869    'P 2 n n':[],
1870    'P n 2 n':[],
1871       
1872    'P m m m':[],
1873       
1874    'P n n n':[],
1875       
1876    'P c c m':[],
1877    'P m a a':[],
1878    'P b m b':[],
1879       
1880    'P b a n':[],
1881    'P n c b':[],
1882    'P c n a':[],
1883       
1884    'P m m a':[],
1885    'P m m b':[],
1886    'P b m m':[],
1887    'P c m m':[],
1888    'P m c m':[],
1889    'P m a m':[],
1890       
1891    'P n n a':[],
1892    'P n n b':[],
1893    'P b n n':[],
1894    'P c n n':[],
1895    'P n c n':[],
1896    'P n a n':[],
1897       
1898    'P m n a':[],
1899    'P n m b':[],
1900    'P b m n':[],
1901    'P c n m':[],
1902    'P n c m':[],
1903    'P m a n':[],
1904       
1905    'P c c a':[],
1906    'P c c b':[],
1907    'P b a a':[],
1908    'P c a a':[],
1909    'P b c b':[],
1910    'P b a b':[],
1911       
1912    'P b a m':[],
1913    'P m c b':[],
1914    'P c m a':[],
1915       
1916    'P c c n':[],
1917    'P n a a':[],
1918    'P b n b':[],
1919       
1920    'P b c m':[],
1921    'P c a m':[],
1922    'P m c a':[],
1923    'P m a b':[],
1924    'P b m a':[],
1925    'P c m b':[],
1926       
1927    'P n n m':[],
1928    'P m n n':[],
1929    'P n m n':[],
1930       
1931    'P m m n':[],
1932    'P n m m':[],
1933    'P m n m':[],
1934       
1935    'P b c n':[],
1936    'P c a n':[],
1937    'P n c a':[],
1938    'P n a b':[],
1939    'P b n a':[],
1940    'P c n b':[],
1941       
1942    'P b c a':[],
1943    'P c a b':[],
1944       
1945    'P n m a':[],
1946    'P m n b':[],
1947    'P b n m':[],
1948    'P c m n':[],
1949    'P m c n':[],
1950    'P n a m':[],
1951       
1952    'C 2 2 21':['(a00)','(0b0)','(00g)','(10g)','(01g)',],
1953    'C 2 2 2':[],
1954    'C m m 2':[],
1955    'C m c 21':[],
1956    'C c c 2':[],
1957       
1958    'C m 2 m':[],
1959    'C 2 m m':[],
1960       
1961    'C m 2 a':[],
1962    'C 2 m b':[],
1963       
1964    'C 2 c m':[],
1965    'C c 2 m':[],
1966    'C 2 c m':[],
1967       
1968    'C m c a':[],
1969    'C m m m':[],
1970    'C c c m':[],
1971    'C m m a':[],
1972    'C c c a':[],
1973    'C m c m':[],
1974    'I 2 2 2':[],
1975    'I 21 21 21':[],
1976    'I m m m':[],
1977       
1978    'I m m 2':[],
1979    'I m 2 m':[],
1980    'I 2 m m':[],
1981       
1982    'I b a 2':[],
1983    'I 2 c b':[],
1984    'I c 2 a':[],
1985       
1986    'I m a 2':[],
1987    'I b m 2':[],
1988    'I 2 m b':[],
1989       
1990    'I 2 c m':[],
1991    'I c 2 m':[],
1992    'I m 2 a':[],
1993       
1994    'I b a m':[],
1995    'I m c b':[],
1996    'I c m a':[],
1997       
1998    'I b c a':[],
1999    'I c a b':[],
2000       
2001    'I m m a':[],
2002    'I m m b':[],
2003    'I b m m ':[],
2004    'I c m m':[],
2005    'I m c m':[],
2006    'I m a m':[],
2007       
2008    'F 2 2 2':[],
2009    'F m m m':[],
2010    'F d d d':[],
2011       
2012    'F m m 2':[],
2013    'F m 2 m':[],
2014    'F 2 m m':[],
2015       
2016    'F d d 2':[],
2017    'F d 2 d':[],
2018    'F 2 d d':[],       
2019    #tetragonal - done
2020    'P 4':['(00g)','(00g)q','(00g)s','(1/21/2g)','(1/21/2g)q',],
2021    'P 41':['(00g)','(00g)q','(00g)s','(1/21/2g)','(1/21/2g)q',],
2022    'P 42':['(00g)','(00g)q','(00g)s','(1/21/2g)','(1/21/2g)q',],
2023    'P 43':['(00g)','(00g)q','(00g)s','(1/21/2g)','(1/21/2g)q',],
2024    'P -4':['(00g)','(1/21/2g)',],
2025    'P 4/m':['(00g)','(00g)s0','(1/21/2g)','(1/21/2g)s0',],
2026    'P 42/m':['(00g)','(00g)s0','(1/21/2g)','(1/21/2g)s0',],
2027    'P 4/n':['(00g)','(00g)s0','(1/21/2g)','(1/21/2g)q0',],
2028    'P 42/n':['(00g)','(00g)s0','(1/21/2g)','(1/21/2g)q0',],
2029    'P 4 2 2':['(00g)','(00g)q00','(00g)s00','(1/21/2g)','(1/21/2g)q00','(1/21/2g)s00',],
2030    'P 4 21 2':['(00g)','(00g)q00','(00g)s00',],
2031    'P 41 2 2':['(00g)','(00g)q00','(00g)s00','(1/21/2g)','(1/21/2g)q00','(1/21/2g)s00',],
2032    'P 41 21 2':['(00g)','(00g)q00','(00g)s00',],
2033    'P 42 2 2':['(00g)','(00g)q00','(00g)s00','(1/21/2g)','(1/21/2g)q00','(1/21/2g)s00',],
2034    'P 42 21 2':['(00g)','(00g)q00','(00g)s00',],
2035    'P 43 2 2':['(00g)','(00g)q00','(00g)s00','(1/21/2g)','(1/21/2g)q00','(1/21/2g)s00',],
2036    'P 43 21 2':['(00g)','(00g)q00','(00g)s00',],
2037    'P 4 m m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',
2038        '(1/21/2g)','(1/21/2g)ss0','(1/21/2g)0ss','(1/21/2g)s0s',],
2039    'P 4 b m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s','(1/21/2g)','(1/21/2g)qq0','(1/21/2g)qqs',],
2040    'P 42 c m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',
2041        '(1/21/2g)','(1/21/2g)ss0','(1/21/2g)0ss','(1/21/2g)s0s',],
2042    'P 42 n m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s','(1/21/2g)','(1/21/2g)qq0','(1/21/2g)qqs',],
2043    'P 4 c c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',
2044        '(1/21/2g)','(1/21/2g)ss0','(1/21/2g)0ss','(1/21/2g)s0s',],
2045    'P 4 n c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s','(1/21/2g)','(1/21/2g)qq0','(1/21/2g)qqs'],
2046    'P 42 m c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',
2047        '(1/21/2g)','(1/21/2g)ss0','(1/21/2g)0ss','(1/21/2g)s0s',],
2048    'P 42 b c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s','(1/21/2g)','(1/21/2g)qq0','(1/21/2g)qqs'],
2049    'P -4 2 m':['(00g)','(00g)00s','(1/21/2g)','(1/21/2g)00s',],
2050    'P -4 2 c':['(00g)','(00g)00s',],
2051    'P -4 21 m':['(00g)','(00g)00s',],
2052    'P -4 21 c':['(00g)','(00g)00s',],
2053    'P -4 m 2':['(00g)','(00g)0s0','(1/21/2g)','(1/21/2g)0s0',],
2054    'P -4 c 2':['(00g)','(00g)0s0','(1/21/2g)','(1/21/2g)0s0',],
2055    'P -4 b 2':['(00g)','(00g)0s0','(1/21/2g)','(1/21/2g)0q0',],
2056    'P -4 n 2':['(00g)','(00g)0s0','(1/21/2g)','(1/21/2g)0q0',],
2057    'P 4/m m m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2058        '(1/21/2g)','(1/21/2g)s0s0','(1/21/2g)00ss','(1/21/2g)s00s',],
2059    'P 4/m c c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2060        '(1/21/2g)','(1/21/2g)s0s0','(1/21/2g)00ss','(1/21/2g)s00s',],
2061    'P 4/n b m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2062        '(1/21/2g)','(1/21/2g)q0q0','(1/21/2g)q0qs',],
2063    'P 4/n n c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2064        '(1/21/2g)','(1/21/2g)q0q0','(1/21/2g)q0qs',],
2065    'P 4/m b m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2066    'P 4/m n c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2067    'P 4/n m m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2068    'P 4/n c c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2069    'P 42/m m c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2070        '(1/21/2g)','(1/21/2g)s0s0','(1/21/2g)00ss','(1/21/2g)s00s',],
2071    'P 42/m c m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2072        '(1/21/2g)','(1/21/2g)s0s0','(1/21/2g)00ss','(1/21/2g)s00s',],
2073    'P 42/n b c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2074        '(1/21/2g)','(1/21/2g)q0q0','(1/21/2g)q0qs',],
2075    'P 42/n n m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',
2076        '(1/21/2g)','(1/21/2g)q0q0','(1/21/2g)q0qs',],
2077    'P 42/m b c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2078    'P 42/m n m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2079    'P 42/n m c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2080    'P 42/n c m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2081    'I 4':['(00g)','(00g)q','(00g)s',],
2082    'I 41':['(00g)','(00g)q','(00g)s',],
2083    'I -4':['(00g)',],
2084    'I 4/m':['(00g)','(00g)s0',],
2085    'I 41/a':['(00g)','(00g)s0',],  #s0?
2086    'I 4 2 2':['(00g)','(00g)q00','(00g)s00',],
2087    'I 41 2 2':['(00g)','(00g)q00','(00g)s00',],
2088    'I 4 m m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2089    'I 4 c m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2090    'I 41 m d':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2091    'I 41 c d':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2092    'I -4 m 2':['(00g)','(00g)0s0',],
2093    'I -4 c 2':['(00g)','(00g)0s0',],
2094    'I -4 2 m':['(00g)','(00g)00s',],
2095    'I -4 2 d':['(00g)','(00g)00s',],
2096    'I 4/m m m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2097    'I 4/m c m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2098    'I 41/a m d':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2099    'I 41/a c d':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2100    #trigonal/rhombahedral - done & checked
2101    'R 3':['(00g)','(00g)t',],
2102    'R -3':['(00g)',],
2103    'R 3 2':['(00g)','(00g)t0',],
2104    'R 3 m':['(00g)','(00g)0s',],
2105    'R 3 c':['(00g)','(00g)0s',],   #not 0s0
2106    'R -3 m':['(00g)','(00g)0s',],
2107    'R -3 c':['(00g)','(00g)0s',],  #not 0s0
2108    'P 3':['(00g)','(00g)t','(1/31/3g)','(1/31/3g)t',],
2109    'P 31':['(00g)','(00g)t','(1/31/3g)','(1/31/3g)t',],
2110    'P 32':['(00g)','(00g)t','(1/31/3g)','(1/31/3g)t',],
2111    'P -3':['(00g)','(1/31/3g)',],
2112    'P 3 1 2':['(00g)','(00g)t00','(1/31/3g)','(1/31/3g)t00',],
2113    'P 3 2 1':['(00g)','(00g)t00',],
2114    'P 31 1 2':['(00g)','(00g)t00','(1/31/3g)','(1/31/3g)t00',],
2115    'P 31 2 1':['(00g)','(00g)t00',],
2116    'P 32 1 2':['(00g)','(00g)t00','(1/31/3g)','(1/31/3g)t00',],
2117    'P 32 2 1':['(00g)','(00g)t00',],
2118    'P 3 m 1':['(00g)','(00g)0s0',],
2119    'P 3 1 m':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
2120    'P 3 c 1':['(00g)','(00g)0s0',],
2121    'P 3 1 c':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
2122    'P -3 1 m':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
2123    'P -3 1 c':['(00g)','(00g)00s','(1/31/3g)','(1/31/3g)00s',],
2124    'P -3 m 1':['(00g)','(00g)0s0',],
2125    'P -3 c 1':['(00g)','(00g)0s0',],
2126    #hexagonal - done & checked
2127    'P 6':['(00g)','(00g)h','(00g)t','(00g)s',],
2128    'P 61':['(00g)','(00g)h','(00g)t','(00g)s',],
2129    'P 65':['(00g)','(00g)h','(00g)t','(00g)s',],
2130    'P 62':['(00g)','(00g)h','(00g)t','(00g)s',],
2131    'P 64':['(00g)','(00g)h','(00g)t','(00g)s',],
2132    'P 63':['(00g)','(00g)h','(00g)t','(00g)s',],
2133    'P -6':['(00g)',],
2134    'P 6/m':['(00g)','(00g)s0',],
2135    'P 63/m':['(00g)','(00g)s0'],
2136    'P 6 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2137    'P 61 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2138    'P 65 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2139    'P 62 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2140    'P 64 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2141    'P 63 2 2':['(00g)','(00g)h00','(00g)t00','(00g)s00',],
2142    'P 6 m m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2143    'P 6 c c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2144    'P 63 c m':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2145    'P 63 m c':['(00g)','(00g)ss0','(00g)0ss','(00g)s0s',],
2146    'P -6 m 2':['(00g)','(00g)0s0',],
2147    'P -6 c 2':['(00g)','(00g)0s0',],
2148    'P -6 2 m':['(00g)','(00g)00s',],
2149    'P -6 2 c':['(00g)','(00g)00s',],
2150    'P 6/m m m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2151    'P 6/m c c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2152    'P 63/m c m':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2153    'P 63/m m c':['(00g)','(00g)s0s0','(00g)00ss','(00g)s00s',],
2154    }
2155
2156#'A few non-standard space groups for test use'
2157nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
2158                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
2159                      'R 3 c r','R -3 c r','R -3 m r',),
2160
2161#A list of orthorhombic space groups that were renamed in the 2002 Volume A,
2162# along with the pre-2002 name. The e designates a double glide-plane'''
2163sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
2164                            ('A e a 2', 'A b a 2',),
2165                            ('C m c e', 'C m c a',),
2166                            ('C m m e', 'C m m a',),
2167                            ('C c c e', 'C c c a'),)
2168#Use the space groups types in this order to list the symbols in the
2169#order they are listed in the International Tables, vol. A'''
2170symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
2171               'trigonal', 'hexagonal', 'cubic')
2172
2173# self-test materials follow. Requires files in directory testinp
2174selftestlist = []
2175'''Defines a list of self-tests'''
2176selftestquiet = True
2177def _ReportTest():
2178    'Report name and doc string of current routine when ``selftestquiet`` is False'
2179    if not selftestquiet:
2180        import inspect
2181        caller = inspect.stack()[1][3]
2182        doc = eval(caller).__doc__
2183        if doc is not None:
2184            print('testing '+__file__+' with '+caller+' ('+doc+')')
2185        else:
2186            print('testing '+__file__()+" with "+caller)
2187def test0():
2188    '''self-test #0: exercise MoveToUnitCell'''
2189    _ReportTest()
2190    msg = "MoveToUnitCell failed"
2191    assert (MoveToUnitCell([1,2,3]) == [0,0,0]).all, msg
2192    assert (MoveToUnitCell([2,-1,-2]) == [0,0,0]).all, msg
2193    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9) < 1e-6, msg
2194    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1) < 1e-6, msg
2195selftestlist.append(test0)
2196
2197def test1():
2198    '''self-test #1: SpcGroup against previous results'''
2199    #'''self-test #1: SpcGroup and SGPrint against previous results'''
2200    _ReportTest()
2201    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
2202    if ospath.exists(testdir):
2203        if testdir not in sys.path: sys.path.insert(0,testdir)
2204    import spctestinp
2205    def CompareSpcGroup(spc, referr, refdict, reflist): 
2206        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
2207        # if an error is reported, the dictionary can be ignored
2208        msg0 = "CompareSpcGroup failed on space group %s" % spc
2209        result = SpcGroup(spc)
2210        if result[0] == referr and referr > 0: return True
2211        keys = result[1].keys()
2212        #print result[1]['SpGrp']
2213        #msg = msg0 + " in list lengths"
2214        #assert len(keys) == len(refdict.keys()), msg
2215        for key in refdict.keys():
2216            if key == 'SGOps' or  key == 'SGCen':
2217                msg = msg0 + (" in key %s length" % key)
2218                assert len(refdict[key]) == len(result[1][key]), msg
2219                for i in range(len(refdict[key])):
2220                    msg = msg0 + (" in key %s level 0" % key)
2221                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
2222                    msg = msg0 + (" in key %s level 1" % key)
2223                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
2224            else:
2225                msg = msg0 + (" in key %s" % key)
2226                assert result[1][key] == refdict[key], msg
2227        msg = msg0 + (" in key %s reflist" % key)
2228        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
2229        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
2230        # for now disable SGPrint testing, output has changed
2231        #assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
2232    for spc in spctestinp.SGdat:
2233        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
2234selftestlist.append(test1)
2235
2236def test2():
2237    '''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
2238    _ReportTest()
2239    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
2240    if ospath.exists(testdir):
2241        if testdir not in sys.path: sys.path.insert(0,testdir)
2242    import sgtbxtestinp
2243    def CompareWcctbx(spcname, cctbx_in, debug=0):
2244        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
2245        cctbx = cctbx_in[:] # make copy so we don't delete from the original
2246        spc = (SpcGroup(spcname))[1]
2247        if debug: print spc['SpGrp']
2248        if debug: print spc['SGCen']
2249        latticetype = spcname.strip().upper()[0]
2250        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
2251        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
2252        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
2253        onebar = [1]
2254        if spc['SGInv']: onebar.append(-1)
2255        for (op,off) in spc['SGOps']:
2256            for inv in onebar:
2257                for cen in spc['SGCen']:
2258                    noff = off + cen
2259                    noff = MoveToUnitCell(noff)
2260                    mult = tuple((op*inv).ravel().tolist())
2261                    if debug: print "\n%s: %s + %s" % (spcname,mult,noff)
2262                    for refop in cctbx:
2263                        if debug: print refop
2264                        # check the transform
2265                        if refop[:9] != mult: continue
2266                        if debug: print "mult match"
2267                        # check the translation
2268                        reftrans = list(refop[-3:])
2269                        reftrans = MoveToUnitCell(reftrans)
2270                        if all(abs(noff - reftrans) < 1.e-5):
2271                            cctbx.remove(refop)
2272                            break
2273                    else:
2274                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
2275    for key in sgtbxtestinp.sgtbx:
2276        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
2277selftestlist.append(test2)
2278
2279def test3(): 
2280    '''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
2281     for selected space groups against info in IT Volume A '''
2282    _ReportTest()
2283    def ExerciseSiteSym (spc, crdlist):
2284        'compare site symmetries and multiplicities for a specified space group'
2285        msg = "failed on site sym test for %s" % spc
2286        (E,S) = SpcGroup(spc)
2287        assert not E, msg
2288        for t in crdlist:
2289            symb, m = SytSym(t[0],S)
2290            if symb.strip() != t[2].strip() or m != t[1]:
2291                print spc,t[0],m,symb,t[2]
2292            assert m == t[1]
2293            #assert symb.strip() == t[2].strip()
2294
2295    ExerciseSiteSym('p 1',[
2296            ((0.13,0.22,0.31),1,'1'),
2297            ((0,0,0),1,'1'),
2298            ])
2299    ExerciseSiteSym('p -1',[
2300            ((0.13,0.22,0.31),2,'1'),
2301            ((0,0.5,0),1,'-1'),
2302            ])
2303    ExerciseSiteSym('C 2/c',[
2304            ((0.13,0.22,0.31),8,'1'),
2305            ((0.0,.31,0.25),4,'2(y)'),
2306            ((0.25,.25,0.5),4,'-1'),
2307            ((0,0.5,0),4,'-1'),
2308            ])
2309    ExerciseSiteSym('p 2 2 2',[
2310            ((0.13,0.22,0.31),4,'1'),
2311            ((0,0.5,.31),2,'2(z)'),
2312            ((0.5,.31,0.5),2,'2(y)'),
2313            ((.11,0,0),2,'2(x)'),
2314            ((0,0.5,0),1,'222'),
2315            ])
2316    ExerciseSiteSym('p 4/n',[
2317            ((0.13,0.22,0.31),8,'1'),
2318            ((0.25,0.75,.31),4,'2(z)'),
2319            ((0.5,0.5,0.5),4,'-1'),
2320            ((0,0.5,0),4,'-1'),
2321            ((0.25,0.25,.31),2,'4(001)'),
2322            ((0.25,.75,0.5),2,'-4(001)'),
2323            ((0.25,.75,0.0),2,'-4(001)'),
2324            ])
2325    ExerciseSiteSym('p 31 2 1',[
2326            ((0.13,0.22,0.31),6,'1'),
2327            ((0.13,0.0,0.833333333),3,'2(100)'),
2328            ((0.13,0.13,0.),3,'2(110)'),
2329            ])
2330    ExerciseSiteSym('R 3 c',[
2331            ((0.13,0.22,0.31),18,'1'),
2332            ((0.0,0.0,0.31),6,'3'),
2333            ])
2334    ExerciseSiteSym('R 3 c R',[
2335            ((0.13,0.22,0.31),6,'1'),
2336            ((0.31,0.31,0.31),2,'3(111)'),
2337            ])
2338    ExerciseSiteSym('P 63 m c',[
2339            ((0.13,0.22,0.31),12,'1'),
2340            ((0.11,0.22,0.31),6,'m(100)'),
2341            ((0.333333,0.6666667,0.31),2,'3m(100)'),
2342            ((0,0,0.31),2,'3m(100)'),
2343            ])
2344    ExerciseSiteSym('I a -3',[
2345            ((0.13,0.22,0.31),48,'1'),
2346            ((0.11,0,0.25),24,'2(x)'),
2347            ((0.11,0.11,0.11),16,'3(111)'),
2348            ((0,0,0),8,'-3(111)'),
2349            ])
2350selftestlist.append(test3)
2351
2352if __name__ == '__main__':
2353    # run self-tests
2354    selftestquiet = False
2355    for test in selftestlist:
2356        test()
2357    print "OK"
Note: See TracBrowser for help on using the repository browser.