source: trunk/GSASIIspc.py @ 1534

Last change on this file since 1534 was 1534, checked in by toby, 9 years ago

fix unit tests
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 81.8 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-24 21:11:20 +0000 (Fri, 24 Oct 2014) $
12# $Author: toby $
13# $Revision: 1534 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 1534 2014-10-24 21:11:20Z toby $
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: 1534 $")
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' and '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            try:
499                if SGData['SGLatt'] == 'R' and modsym.index('1/3'):
500                    return False
501            except ValueError:
502                pass
503            return True
504        elif LaueId in [10,11,] and LaueModId in [19,]:
505            return True
506        return False
507       
508    def fixMonoOrtho():
509        mod = ''.join(modsym).replace('1/2','0').replace('1','0')
510        if SGData['SGPtGrp'] in ['2','m']:  #OK
511            if mod in ['a00','0b0','00g']:
512                result = [i*-1 for i in SGData['SSGKl']]
513            else:
514                result = SGData['SSGKl'][:]
515            if '/' in mod:
516                return [i*-1 for i in result]
517            else:
518                return result
519        elif SGData['SGPtGrp'] == '2/m':    #OK
520            if mod in ['a00','0b0','00g']:
521                result =  SGData['SSGKl'][:]
522            else:
523                result = [i*-1 for i in SGData['SSGKl']]
524            if '/' in mod:
525                return [i*-1 for i in result]
526            else:
527                return result
528        else:   #orthorhombic
529            if SGData['SGPtGrp'] == '222':
530                return [1 if i in ['a','b','g'] else -1 for i in mod]
531            elif SGData['SGPtGrp'] == 'mm2':
532                if 'g' in mod:
533                    return [1,1,1]
534                elif 'b' in mod:
535                    return [1,-1,-1]
536                else:
537                    return [-1,1,-1]
538            elif SGData['SGPtGrp'] == 'm2m':
539                if 'b' in mod:
540                    return [1,1,1]
541                elif 'g' in mod:
542                    return [1,-1,-1]
543                else:
544                    return [-1,-1,1]               
545            elif SGData['SGPtGrp'] == '2mm':
546                if 'a' in mod:
547                    return [1,1,1]
548                elif 'b' in mod:
549                    return [-1,-1,1]
550                else:
551                    return [-1,1,-1]
552            else:
553                return [-1 if i in ['a','b','g'] else 1 for i in mod]
554               
555    def extendSSGOps(SSGOps):
556        nOps = len(SSGOps)
557        for i in range(nOps):
558            if np.allclose(SSGOps[i][0][3],np.zeros(4)):
559                continue
560            for j in range(nOps):
561                if np.allclose(SSGOps[j][0][3],np.zeros(4)):
562                    continue
563                OpC = list(SGProd(SSGOps[j],SSGOps[i]))
564                OpC[1] %= 1.
565                for k in range(nOps):
566                    OpD = SSGOps[k]
567                    if SSMT2text(OpC) == SSMT2text(OpD):
568                        continue
569                    elif np.allclose(OpC[0][:3,:3],OpD[0][:3,:3]):
570                        if np.allclose(OpD[0][3],np.zeros(4)):
571                            SSGOps[k] = OpC
572                        elif np.any([np.allclose(OpC[0][3][:3],cen) for cen in SGData['SGCen']]):   #?
573                            continue
574                        else:
575                            OpCtxt = SSMT2text(OpC).replace(' ','')
576                            OpDtxt = SSMT2text(OpD).replace(' ','')
577                            print 'OpC',OpCtxt,'OpD',OpDtxt
578                            return False,OpCtxt+' conflict with '+OpDtxt
579        return True,SSGOps
580       
581               
582    def genSSGOps():
583        SSGOps = SSGData['SSGOps'][:]
584        iFrac = {}
585        for i,frac in enumerate(SSGData['modSymb']):
586            if frac in ['1/2','1/3','1/4','1/6','1']:
587                iFrac[i] = frac
588        print SGData['SpGrp']+SSymbol
589        print 'SSGKl',SSGKl,'genQ',genQ,'iFrac',iFrac
590# set identity & 1,-1; triclinic
591        SSGOps[0][0][3,3] = 1.
592# expand if centrosymmetric
593        if SGData['SGInv']:
594            SSGOps += [[-1*M,V] for M,V in SSGOps[:]]
595# monoclinic
596        if SGData['SGPtGrp'] in ['2','m']:  #OK
597            SSGOps[1][0][3,3] = SSGKl[0]
598            SSGOps[1][1][3] = genQ[0]
599            for i in iFrac:
600                SSGOps[1][0][3,i] = -SSGKl[0]
601        elif SGData['SGPtGrp'] == '2/m':
602            for i,j in enumerate([1,3]):
603                SSGOps[j][0][3,3] = -SSGKl[i]
604                if genQ[i]:
605                    SSGOps[j][1][3] = genQ[i]
606                for k in iFrac:
607                    SSGOps[j][0][3,k] = SSGKl[i]
608                E,SSGOps = extendSSGOps(SSGOps)
609            print E,SSMT2text(SSGOps[1]).replace(' ',''),SSMT2text(SSGOps[3]).replace(' ','')
610           
611# orthorhombic
612        elif SGData['SGPtGrp'] in ['222','mm2','m2m','2mm','mmm']:
613            for i in [0,1,2]:
614                SSGOps[i+1][0][3,3] = SSGKl[i]
615                SSGOps[i+1][1][3] = genQ[i]
616            for i in iFrac:
617                SSGOps[1][0][3,i] = -1
618            print SSMT2text(SSGOps[1]).replace(' ',''),SSMT2text(SSGOps[2]).replace(' ',''), \
619                SSMT2text(SSGOps[3]).replace(' ','')
620# tetragonal
621        elif SGData['SGPtGrp'] == '4':  #OK
622            SSGOps[1][0][3,3] = SSGKl[0]
623            SSGOps[1][1][3] = genQ[0]
624            if '1/2' in SSGData['modSymb']:
625                SSGOps[1][0][3,1] = -1
626        elif SGData['SGPtGrp'] == '-4': #OK
627            SSGOps[1][0][3,3] = SSGKl[0]
628            if '1/2' in SSGData['modSymb']:
629                SSGOps[1][0][3,1] = 1
630        elif SGData['SGPtGrp'] in ['4/m',]:
631            SSGOps[1][0][3,3] = SSGKl[0]
632            SSGOps[5][0][3,3] = SSGKl[1]
633            SSGOps[1][1][3] = genQ[0]
634            if '1/2' in SSGData['modSymb']:
635                SSGOps[1][0][3,1] = -1
636        elif SGData['SGPtGrp'] in ['422','4mm','-42m','-4m2',]:
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            SSGOps[1][0][3,3] = SSGKl[1]
644            if '1/2' in SSGData['modSymb']:
645                SSGOps[1][0][3,3] *= -1
646                SSGOps[1][0][3,1] = -1
647            SSGOps[2][0][3,3] = SSGKl[1]
648            SSGOps[3][0][3,3] = SSGKl[2]
649            SSGOps[4][0][3,3] = SSGKl[3]
650# trigonal
651        elif SGData['SGPtGrp'] == '3':
652            SSGOps[1][0][3,3] = SSGKl[0]
653            SSGOps[1][1][3] = genQ[0]
654        elif SGData['SGPtGrp'] == '-3':
655            SSGOps[1][0][3,3] = -SSGKl[0]
656        elif SGData['SGPtGrp'] in ['32','3m','-3m',]:
657            SSGOps[1][0][3,3] = SSGKl[1]
658        elif SGData['SGPtGrp'] in ['312','321','3m1','31m','-3m1','-31m',]:
659            SSGOps[1][0][3,3] = SSGKl[0]
660# hexagonal
661        elif SGData['SGPtGrp'] == '6':  #OK
662            SSGOps[1][0][3,3] = SSGKl[0]
663            SSGOps[1][1][3] = genQ[0]
664        elif SGData['SGPtGrp'] == '-6': #OK
665            SSGOps[1][0][3,3] = SSGKl[0]
666        elif SGData['SGPtGrp'] in ['6/m',]: #OK
667            SSGOps[1][0][3,3] = -SSGKl[1]
668            SSGOps[1][1][3] = genQ[0]
669            SSGOps[2][1][3] = genQ[1]
670        elif SGData['SGPtGrp'] in ['622','6mm','-62m','-6m2',]: #OK
671            for i,j in enumerate([1,10,11]):
672                SSGOps[j][0][3,3] = SSGKl[i]
673                if genQ[i]:
674                    SSGOps[j][1][3] = genQ[i]
675                E,SSGOps = extendSSGOps(SSGOps)
676        elif SGData['SGPtGrp'] in ['6/mmm',]: #not OK
677            for i,j in enumerate([1,15,19,11]):
678                SSGOps[j][0][3,3] = SSGKl[i]
679                if genQ[i]:
680                    SSGOps[j][1][3] = genQ[i]
681                E,SSGOps = extendSSGOps(SSGOps)
682        if SGData['SGPtGrp'] in ['1','-1']: #triclinic - done
683            return True,SSGOps
684        E,SSGOps = extendSSGOps(SSGOps)
685        return E,SSGOps
686       
687    def specialGen(gensym):
688        sym = ''.join(gensym)
689        if SGData['SGPtGrp'] in ['-62m',] and sym == '00s':
690            gensym = '0ss'
691        return gensym
692                   
693    def checkGen(gensym):
694        sym = ''.join(gensym)
695# monoclinic - all done
696        if str(SSGKl) == '[-1]' and sym == 's':
697            return False
698        elif str(SSGKl) == '[-1, 1]' and sym == 's0':
699            return False
700        elif str(SSGKl) == '[1, -1]' and sym == '0s':
701            return False
702#orthorhombic - all
703        elif SGData['SGPtGrp'] in ['222',] and sym not in ['','s00','0s0','00s']:
704            return False 
705        elif SGData['SGPtGrp'] in ['2mm','m2m','mm2','mmm'] and sym not in GenSymList[4:15]:
706            return False 
707#tetragonal - all done
708        elif SGData['SGPtGrp'] in ['4',] and sym not in ['','s','q']:
709            return False 
710        elif SGData['SGPtGrp'] in ['-4',] and sym not in ['',]:
711            return False             
712        elif SGData['SGPtGrp'] in ['4/m',] and sym not in ['','s0','q0']:
713            return False
714        elif SGData['SGPtGrp'] in ['422',] and sym not in ['','q00','s00']:
715            return False         
716        elif SGData['SGPtGrp'] in ['4mm',] and sym not in ['','ss0','s0s','0ss','qq0','qqs']:
717            return False
718        elif SGData['SGPtGrp'] in ['-4m2',] and sym not in ['','00s','00q']:
719            return False
720        elif SGData['SGPtGrp'] in ['-42m',] and sym not in ['','0s0','0q0']:
721            return False
722        elif SGData['SGPtGrp'] in ['4/mmm',] and sym not in ['','s00s','s0s0','00ss','q0q0','q0qs']:
723            return False
724#trigonal/rhombohedral - all done
725        elif SGData['SGPtGrp'] in ['3',] and sym not in ['','t']:
726            return False 
727        elif SGData['SGPtGrp'] in ['-3',] and sym not in ['',]:
728            return False 
729        elif SGData['SGPtGrp'] in ['32',] and sym not in ['','t0']:
730            return False 
731        elif SGData['SGPtGrp'] in ['321',] and sym not in ['','t00']:
732            return False 
733        elif SGData['SGPtGrp'] in ['312',] and sym not in ['',]:
734            return False 
735        elif SGData['SGPtGrp'] in ['3m','-3m'] and sym not in ['','0s']:
736            return False 
737        elif SGData['SGPtGrp'] in ['3m1','-3m1'] and sym not in ['','0s0']:
738            return False 
739        elif SGData['SGPtGrp'] in ['31m','-31m'] and sym not in ['','00s']:
740            return False 
741#hexagonal - all done
742        elif SGData['SGPtGrp'] in ['6',] and sym not in ['','s','h','t']:
743            return False 
744        elif SGData['SGPtGrp'] in ['-6',] and sym not in ['',]:
745            return False
746        elif SGData['SGPtGrp'] in ['6/m',] and sym not in ['','s0']:
747            return False
748        elif SGData['SGPtGrp'] in ['622',] and sym not in ['','h00','t00','s00']:
749            return False         
750        elif SGData['SGPtGrp'] in ['6mm',] and sym not in ['','ss0','s0s','0ss']:
751            return False
752        elif SGData['SGPtGrp'] in ['-6m2',] and sym not in ['','0s0']:
753            return False
754        elif SGData['SGPtGrp'] in ['-62m',] and sym not in ['','0ss']:
755            return False
756        elif SGData['SGPtGrp'] in ['6/mmm',] and sym not in ['','s00s','s0s0','00ss']:
757            return False
758        return True
759       
760    LaueModList = ['abg', 'ab0', 'ab1/2', 'a0g', 'a1/2g','0bg', '1/2bg',
761               'a00', 'a01/2', 'a1/20', 'a1/21/2', 'a01', 'a10', 
762               '0b0', '0b1/2', '1/2b0', '1/2b1/2', '0b1', '1b0',
763               '00g', '01/2g', '1/20g', '1/21/2g', '01g', '10g','1/31/3g']
764    LaueList = ['-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m']
765    GenSymList = ['','s','0s','s0','00s','0s0','s00','s0s','ss0','0ss','q00','0q0','00q','qq0','q0q','0qq',
766        'q','q0','0q','qqs','s0s0','00ss','s00s','q0q0','q0qs','t','t00','t0','h','h00']
767    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.}
768    LaueId = LaueList.index(SGData['SGLaue'])
769    if SGData['SGLaue'] in ['m3','m3m']:
770        return '(3+1) superlattices not defined for cubic space groups',None
771    elif SGData['SGLaue'] in ['3R','3mR']:
772        return '(3+1) superlattices not defined for rhombohedral settings - use hexagonal setting',None
773    try:
774        modsym,gensym = splitSSsym(SSymbol)
775    except ValueError:
776        return 'Error in superspace symbol '+SSymbol,None
777    if ''.join(gensym) not in GenSymList:
778        return 'unknown generator symbol '+''.join(gensym),None
779    try:
780        LaueModId = LaueModList.index(''.join(modsym))
781    except ValueError:
782        return 'Unknown modulation symbol '+''.join(modsym),None
783    if not checkModSym():
784        return 'Modulation '+''.join(modsym)+' not consistent with space group '+SGData['SpGrp'],None
785    modQ = [Fracs[mod] for mod in modsym]
786    SSGKl = SGData['SSGKl'][:]
787    if SGData['SGLaue'] in ['2/m','mmm']:
788        SSGKl = fixMonoOrtho()
789    if len(gensym) and len(gensym) != len(SSGKl):
790        return 'Wrong number of items in generator symbol '+''.join(gensym),None
791    gensym = specialGen(gensym)
792    if not checkGen(gensym):
793        return 'Generator '+''.join(gensym)+' not consistent with space group '+SGData['SpGrp'],None
794    genQ = [Fracs[mod] for mod in gensym]
795    if not genQ:
796        genQ = [0,0,0,0]
797    SSGData = {'SSpGrp':SGData['SpGrp']+SSymbol,'modQ':modQ,'modSymb':modsym}
798    SSCen = np.ones((len(SGData['SGCen']),4))*0.5
799    for icen,cen in enumerate(SGData['SGCen']):
800        SSCen[icen,0:3] = cen
801    SSGData['SSGCen'] = SSCen
802    SSGData['SSGOps'] = []
803    for iop,op in enumerate(SGData['SGOps']):
804        T = np.zeros(4)
805        ssop = np.zeros((4,4))
806        ssop[:3,:3] = op[0]
807        T[:3] = op[1]
808        SSGData['SSGOps'].append([ssop,T])
809    E,Result = genSSGOps()
810    if E:
811        SSGData['SSGOps'] = Result                     
812        return None,SSGData
813    else:
814        return Result+'\nOperator conflict - incorrect superspace symbol',None
815
816def SSGPrint(SGData,SSGData):
817    '''
818    Print the output of SSpcGroup in a nicely formatted way. Used in SSpaceGroup
819
820    :param SGData: space group data structure as defined in SpcGroup above.
821    :param SSGData: from :func:`SSpcGroup`
822    :returns:
823        SSGText - list of strings with the superspace group details
824        SGTable - list of strings for each of the operations
825    '''
826    Mult = len(SSGData['SSGCen'])*len(SSGData['SSGOps'])
827    SSGText = []
828    SSGText.append(' Superspace Group: '+SSGData['SSpGrp'])
829    CentStr = 'centrosymmetric'
830    if not SGData['SGInv']:
831        CentStr = 'non'+CentStr
832    if SGData['SGLatt'] in 'ABCIFR':
833        SSGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
834    else:
835        SSGText.append(' The superlattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())       
836    SSGText.append(' The Laue symmetry is '+SGData['SGLaue'])
837    SSGText.append(' The superlattice point group is '+SGData['SGPtGrp']+','+''.join([str(i) for i in SGData['SSGKl']]))
838    SSGText.append(' The number of superspace group generators is '+str(len(SGData['SSGKl'])))
839    SSGText.append(' Multiplicity of a general site is '+str(Mult))
840    if SGData['SGUniq'] in ['a','b','c']:
841        SSGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
842    if SGData['SGInv']:
843        SSGText.append(' The inversion center is located at 0,0,0')
844    if SGData['SGPolax']:
845        SSGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
846    SSGText.append(' ')
847    if len(SSGData['SSGCen']) > 1:
848        SSGText.append(' The equivalent positions are:')
849        SSGText.append(' ('+SSLatt2text(SSGData['SSGCen'])+')+\n')
850    else:
851        SSGText.append(' The equivalent positions are:\n')
852    SSGTable = []
853    for i,Opr in enumerate(SSGData['SSGOps']):
854        SSGTable.append('(%2d) %s'%(i+1,SSMT2text(Opr)))
855    return SSGText,SSGTable
856   
857def SSGModCheck(Vec,SSGData):
858    ''' Checks modulation vector compatibility with supersymmetry space group symbol.
859    Superspace group symbol takes precidence & the vector will be modified accordingly
860    '''
861    modQ = SSGData['modQ']
862    modSymb = SSGData['modSymb']
863    Vec = [0.1 if (vec == 0.0 and mod in ['a','b','g']) else vec for [vec,mod] in zip(Vec,modSymb)]
864    return [Q if mod not in ['a','b','g'] and vec != Q else vec for [vec,mod,Q] in zip(Vec,modSymb,modQ)]
865
866def SSMT2text(Opr):
867    "From superspace group matrix/translation operator returns text version"
868    XYZS = ('x','y','z','t')    #Stokes, Campbell & van Smaalen notation
869    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
870    Fld = ''
871    M,T = Opr
872    for j in range(4):
873        IJ = ''
874        for k in range(4):
875            txt = str(int(round(M[j][k])))
876            txt = txt.replace('1',XYZS[k]).replace('0','')
877            if IJ and M[j][k] > 0:
878                IJ += '+'+txt
879            else:
880                IJ += txt
881        IK = int(round(T[j]*12))%12
882        if IK:
883            if not IJ:
884                break
885            if IJ[0] == '-':
886                Fld += (TRA[IK]+IJ).rjust(8)
887            else:
888                Fld += (TRA[IK]+'+'+IJ).rjust(8)
889        else:
890            Fld += IJ.rjust(8)
891        if j != 3: Fld += ', '
892    return Fld
893   
894def SSLatt2text(SSGCen):
895    "Lattice centering vectors to text"
896    lattTxt = ''
897    for vec in SSGCen:
898        lattTxt += ' '
899        for item in vec:
900            if int(item*12.):
901                lattTxt += '1/%d,'%(12/int(item*12))
902            else:
903                lattTxt += '0,'
904        lattTxt = lattTxt.rstrip(',')
905        lattTxt += ';'
906    lattTxt = lattTxt.rstrip(';').lstrip(' ')
907    return lattTxt
908       
909def SSpaceGroup(SGSymbol,SSymbol):
910    '''
911    Print the output of SSpcGroup in a nicely formatted way.
912
913    :param SGSymbol: space group symbol with spaces between axial fields.
914    :param SSymbol: superspace group symbol extension (string).
915    :returns: nothing
916    '''
917
918    E,A = SpcGroup(SGSymbol)
919    if E > 0:
920        print SGErrors(E)
921        return
922    E,B = SSpcGroup(A,SSymbol)   
923    if E > 0:
924        print E
925        return
926    for l in SSGPrint(B):
927        print l
928       
929def SGProd(OpA,OpB):
930    '''
931    Form space group operator product. OpA & OpB are [M,V] pairs;
932        both must be of same dimension (3 or 4). Returns [M,V] pair
933    '''
934    A,U = OpA
935    B,V = OpB
936    M = np.inner(B.T,A)
937    W = np.inner(B,U)+V
938    return M.T,W
939       
940def MoveToUnitCell(xyz):
941    '''
942    Translates a set of coordinates so that all values are >=0 and < 1
943
944    :param xyz: a list or numpy array of fractional coordinates
945    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
946    '''
947    XYZ = np.zeros(3)
948    for i,x in enumerate(xyz):
949        XYZ[i] = (x-int(x))%1.0
950    return XYZ
951       
952def Opposite(XYZ,toler=0.0002):
953    '''
954    Gives opposite corner, edge or face of unit cell for position within tolerance.
955        Result may be just outside the cell within tolerance
956
957    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
958    :param toler: unit cell fraction tolerance making opposite
959    :returns:
960        XYZ: array of opposite positions; always contains XYZ
961    '''
962    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]]
963    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
964    perm = TB*perm3
965    cperm = ['%d%d%d'%(i,j,k) for i,j,k in perm]
966    D = dict(zip(cperm,perm))
967    new = []
968    for key in D:
969        new.append(np.array(D[key])+np.array(XYZ))
970    return new
971       
972def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
973    '''
974    Generates the equivalent positions for a specified coordinate and space group
975
976    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
977    :param SGData: from :func:`SpcGroup`
978    :param All: True return all equivalent positions including duplicates;
979      False return only unique positions
980    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
981    :param Move: True move generated atom positions to be inside cell
982      False do not move atoms       
983    :return: [[XYZEquiv],Idup,[UijEquiv]]
984
985      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
986      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
987      * is centering operator number (1-4) and - is for inversion
988        Cell = unit cell translations needed to put new positions inside cell
989        [UijEquiv] - equivalent Uij; absent if no Uij given
990       
991    '''
992    XYZEquiv = []
993    UijEquiv = []
994    Idup = []
995    Cell = []
996    X = np.array(XYZ)
997    if Move:
998        X = MoveToUnitCell(X)
999    for ic,cen in enumerate(SGData['SGCen']):
1000        C = np.array(cen)
1001        for invers in range(int(SGData['SGInv']+1)):
1002            for io,[M,T] in enumerate(SGData['SGOps']):
1003                idup = ((io+1)+100*ic)*(1-2*invers)
1004                XT = np.inner(M,X)+T
1005                if len(Uij):
1006                    U = Uij2U(Uij)
1007                    U = np.inner(M,np.inner(U,M).T)
1008                    newUij = U2Uij(U)
1009                if invers:
1010                    XT = -XT
1011                XT += C
1012                if Move:
1013                    newX = MoveToUnitCell(XT)
1014                else:
1015                    newX = XT
1016                cell = np.asarray(np.rint(newX-XT),dtype=np.int32)
1017                if All:
1018                    if np.allclose(newX,X,atol=0.0002):
1019                        idup = False
1020                else:
1021                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
1022                        idup = False
1023                if All or idup:
1024                    XYZEquiv.append(newX)
1025                    Idup.append(idup)
1026                    Cell.append(cell)
1027                    if len(Uij):
1028                        UijEquiv.append(newUij)                   
1029    if len(Uij):
1030        return zip(XYZEquiv,UijEquiv,Idup,Cell)
1031    else:
1032        return zip(XYZEquiv,Idup,Cell)
1033
1034def GenHKLf(HKL,SGData):
1035    '''
1036    Uses old GSAS Fortran routine genhkl.for
1037
1038    :param HKL:  [h,k,l]
1039    :param SGData: space group data obtained from SpcGroup
1040    :returns: iabsnt,mulp,Uniq,phi
1041
1042     *   iabsnt = True if reflection is forbidden by symmetry
1043     *   mulp = reflection multiplicity including Friedel pairs
1044     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
1045
1046    '''
1047    hklf = HKL+[0,]
1048    Ops = SGData['SGOps']
1049    OpM = np.array([op[0] for op in Ops])
1050    OpT = np.array([op[1] for op in Ops])
1051    Inv = SGData['SGInv']
1052    Cen = np.array([cen for cen in SGData['SGCen']])
1053   
1054    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
1055    h,k,l,f = Uniq
1056    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
1057    phi = f[:Nuniq]
1058   
1059    return iabsnt,mulp,Uniq,phi
1060                                 
1061def GetOprPtrName(key):
1062    'Needs a doc string'
1063    OprPtrName = {
1064        '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
1065        '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
1066        '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
1067        '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
1068        '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
1069        '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
1070        '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
1071        '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
1072        '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
1073        '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
1074        '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
1075        '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
1076        '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
1077        '2223' :[  68,' 3+++2',39],
1078        '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
1079        '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
1080        '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
1081        '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
1082        '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
1083        '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
1084        '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
1085        '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
1086        '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
1087        }
1088    return OprPtrName[key]
1089
1090def GetKNsym(key):
1091    'Needs a doc string'
1092    KNsym = {
1093        '0'         :'    1   ','1'         :'   -1   ','64'        :'    2(x)','32'        :'    m(x)',
1094        '97'        :'  2/m(x)','16'        :'    2(y)','8'         :'    m(y)','25'        :'  2/m(y)',
1095        '2'         :'    2(z)','4'         :'    m(z)','7'         :'  2/m(z)','134217728' :'   2(yz)',
1096        '67108864'  :'   m(yz)','201326593' :' 2/m(yz)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
1097        '3145729'   :'2/m(0+-)','8388608'   :'   2(xz)','4194304'   :'   m(xz)','12582913'  :' 2/m(xz)',
1098        '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'   2(xy)',
1099        '512'       :'   m(xy)','1537'      :' 2/m(xy)','256'       :'  2(+-0)','128'       :'  m(+-0)',
1100        '385'       :'2/m(+-0)','76'        :'  mm2(x)','52'        :'  mm2(y)','42'        :'  mm2(z)',
1101        '135266336' :' mm2(yz)','69206048'  :'mm2(0+-)','8650760'   :' mm2(xz)','4718600'   :'mm2(+0-)',
1102        '1156'      :' mm2(xy)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'  222(x)',
1103        '8912912'   :'  222(y)','1282'      :'  222(z)','127'       :'  mmm   ','204472417' :'  mmm(x)',
1104        '13369369'  :'  mmm(y)','1927'      :'  mmm(z)','33554496'  :'  4(100)','16777280'  :' -4(100)',
1105        '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
1106        '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
1107        '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
1108        '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
1109        '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
1110        '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
1111        '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
1112        '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
1113        '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
1114        '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
1115        '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
1116        '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
1117        '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
1118        '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
1119        '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
1120        '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
1121        '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
1122        '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
1123        '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
1124        '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
1125        '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
1126        '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
1127        '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
1128        '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
1129        '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
1130        '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   '
1131        }
1132    return KNsym[key]       
1133
1134def GetNXUPQsym(siteSym):       
1135    'Needs a doc string'
1136    NXUPQsym = {
1137        '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'    2(x)':(12,18,12,25),'    m(x)':(25,18,12,25),
1138        '  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),
1139        '    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),
1140        '   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),
1141        '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),
1142        '  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),
1143        '   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),
1144        '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),
1145        ' mm2(yz)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),' mm2(xz)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
1146        ' mm2(xy)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'  222(x)':( 1,13, 0,-1),
1147        '  222(y)':( 1,12, 0,-1),'  222(z)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'  mmm(x)':( 1,13, 0,-1),
1148        '  mmm(y)':( 1,12, 0,-1),'  mmm(z)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
1149        '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),
1150        '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),
1151        '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
1152        '  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),
1153        '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
1154        ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
1155        '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
1156        '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
1157        ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
1158        ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
1159        '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
1160        ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
1161        '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
1162        '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
1163        '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
1164        ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
1165        '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
1166        'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
1167        'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
1168        'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
1169        '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),
1170        '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),
1171        ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
1172        ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
1173        ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
1174        ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
1175        }
1176    return NXUPQsym[siteSym]
1177
1178def GetCSxinel(siteSym): 
1179    'Needs a doc string'
1180    CSxinel = [[],                         # 0th empty - indices are Fortran style
1181        [[0,0,0],[ 0.0, 0.0, 0.0]],      #1  0  0  0
1182        [[1,1,1],[ 1.0, 1.0, 1.0]],      #2  X  X  X
1183        [[1,1,1],[ 1.0, 1.0,-1.0]],      #3  X  X -X
1184        [[1,1,1],[ 1.0,-1.0, 1.0]],      #4  X -X  X
1185        [[1,1,1],[ 1.0,-1.0,-1.0]],      #5 -X  X  X
1186        [[1,1,0],[ 1.0, 1.0, 0.0]],      #6  X  X  0
1187        [[1,1,0],[ 1.0,-1.0, 0.0]],      #7  X -X  0
1188        [[1,0,1],[ 1.0, 0.0, 1.0]],      #8  X  0  X
1189        [[1,0,1],[ 1.0, 0.0,-1.0]],      #9  X  0 -X
1190        [[0,1,1],[ 0.0, 1.0, 1.0]],      #10  0  Y  Y
1191        [[0,1,1],[ 0.0, 1.0,-1.0]],      #11 0  Y -Y
1192        [[1,0,0],[ 1.0, 0.0, 0.0]],      #12  X  0  0
1193        [[0,1,0],[ 0.0, 1.0, 0.0]],      #13  0  Y  0
1194        [[0,0,1],[ 0.0, 0.0, 1.0]],      #14  0  0  Z
1195        [[1,1,0],[ 1.0, 2.0, 0.0]],      #15  X 2X  0
1196        [[1,1,0],[ 2.0, 1.0, 0.0]],      #16 2X  X  0
1197        [[1,1,2],[ 1.0, 1.0, 1.0]],      #17  X  X  Z
1198        [[1,1,2],[ 1.0,-1.0, 1.0]],      #18  X -X  Z
1199        [[1,2,1],[ 1.0, 1.0, 1.0]],      #19  X  Y  X
1200        [[1,2,1],[ 1.0, 1.0,-1.0]],      #20  X  Y -X
1201        [[1,2,2],[ 1.0, 1.0, 1.0]],      #21  X  Y  Y
1202        [[1,2,2],[ 1.0, 1.0,-1.0]],      #22  X  Y -Y
1203        [[1,2,0],[ 1.0, 1.0, 0.0]],      #23  X  Y  0
1204        [[1,0,2],[ 1.0, 0.0, 1.0]],      #24  X  0  Z
1205        [[0,1,2],[ 0.0, 1.0, 1.0]],      #25  0  Y  Z
1206        [[1,1,2],[ 1.0, 2.0, 1.0]],      #26  X 2X  Z
1207        [[1,1,2],[ 2.0, 1.0, 1.0]],      #27 2X  X  Z
1208        [[1,2,3],[ 1.0, 1.0, 1.0]],      #28  X  Y  Z
1209        ]
1210    indx = GetNXUPQsym(siteSym)
1211    return CSxinel[indx[0]]
1212   
1213def GetCSuinel(siteSym):
1214    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
1215    CSuinel = [[],                                             # 0th empty - indices are Fortran style
1216        [[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
1217        [[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
1218        [[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
1219        [[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
1220        [[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
1221        [[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
1222        [[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
1223        [[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
1224        [[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
1225        [[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
1226        [[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
1227        [[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
1228        [[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
1229        [[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
1230        [[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
1231        [[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
1232        [[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
1233        [[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
1234        [[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
1235        [[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
1236        [[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
1237        [[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
1238        [[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
1239        [[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
1240        [[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
1241        [[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
1242        [[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
1243        [[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
1244        [[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
1245        ]
1246    indx = GetNXUPQsym(siteSym)
1247    return CSuinel[indx[1]]
1248   
1249def MustrainNames(SGData):
1250    'Needs a doc string'
1251    laue = SGData['SGLaue']
1252    uniq = SGData['SGUniq']
1253    if laue in ['m3','m3m']:
1254        return ['S400','S220']
1255    elif laue in ['6/m','6/mmm','3m1']:
1256        return ['S400','S004','S202']
1257    elif laue in ['31m','3']:
1258        return ['S400','S004','S202','S211']
1259    elif laue in ['3R','3mR']:
1260        return ['S400','S220','S310','S211']
1261    elif laue in ['4/m','4/mmm']:
1262        return ['S400','S004','S220','S022']
1263    elif laue in ['mmm']:
1264        return ['S400','S040','S004','S220','S202','S022']
1265    elif laue in ['2/m']:
1266        SHKL = ['S400','S040','S004','S220','S202','S022']
1267        if uniq == 'a':
1268            SHKL += ['S013','S031','S211']
1269        elif uniq == 'b':
1270            SHKL += ['S301','S103','S121']
1271        elif uniq == 'c':
1272            SHKL += ['S130','S310','S112']
1273        return SHKL
1274    else:
1275        SHKL = ['S400','S040','S004','S220','S202','S022']
1276        SHKL += ['S310','S103','S031','S130','S301','S013']
1277        SHKL += ['S211','S121','S112']
1278        return SHKL
1279
1280def HStrainNames(SGData):
1281    'Needs a doc string'
1282    laue = SGData['SGLaue']
1283    uniq = SGData['SGUniq']
1284    if laue in ['m3','m3m']:
1285        return ['D11','eA']         #add cubic strain term
1286    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1287        return ['D11','D33']
1288    elif laue in ['3R','3mR']:
1289        return ['D11','D12']
1290    elif laue in ['4/m','4/mmm']:
1291        return ['D11','D33']
1292    elif laue in ['mmm']:
1293        return ['D11','D22','D33']
1294    elif laue in ['2/m']:
1295        Dij = ['D11','D22','D33']
1296        if uniq == 'a':
1297            Dij += ['D23']
1298        elif uniq == 'b':
1299            Dij += ['D13']
1300        elif uniq == 'c':
1301            Dij += ['D12']
1302        return Dij
1303    else:
1304        Dij = ['D11','D22','D33','D12','D13','D23']
1305        return Dij
1306   
1307def MustrainCoeff(HKL,SGData):
1308    'Needs a doc string'
1309    #NB: order of terms is the same as returned by MustrainNames
1310    laue = SGData['SGLaue']
1311    uniq = SGData['SGUniq']
1312    h,k,l = HKL
1313    Strm = []
1314    if laue in ['m3','m3m']:
1315        Strm.append(h**4+k**4+l**4)
1316        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
1317    elif laue in ['6/m','6/mmm','3m1']:
1318        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
1319        Strm.append(l**4)
1320        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
1321    elif laue in ['31m','3']:
1322        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
1323        Strm.append(l**4)
1324        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
1325        Strm.append(4.0*h*k*l*(h+k))
1326    elif laue in ['3R','3mR']:
1327        Strm.append(h**4+k**4+l**4)
1328        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
1329        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
1330        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
1331    elif laue in ['4/m','4/mmm']:
1332        Strm.append(h**4+k**4)
1333        Strm.append(l**4)
1334        Strm.append(3.0*(h*k)**2)
1335        Strm.append(3.0*((h*l)**2+(k*l)**2))
1336    elif laue in ['mmm']:
1337        Strm.append(h**4)
1338        Strm.append(k**4)
1339        Strm.append(l**4)
1340        Strm.append(3.0*(h*k)**2)
1341        Strm.append(3.0*(h*l)**2)
1342        Strm.append(3.0*(k*l)**2)
1343    elif laue in ['2/m']:
1344        Strm.append(h**4)
1345        Strm.append(k**4)
1346        Strm.append(l**4)
1347        Strm.append(3.0*(h*k)**2)
1348        Strm.append(3.0*(h*l)**2)
1349        Strm.append(3.0*(k*l)**2)
1350        if uniq == 'a':
1351            Strm.append(2.0*k*l**3)
1352            Strm.append(2.0*l*k**3)
1353            Strm.append(4.0*k*l*h**2)
1354        elif uniq == 'b':
1355            Strm.append(2.0*l*h**3)
1356            Strm.append(2.0*h*l**3)
1357            Strm.append(4.0*h*l*k**2)
1358        elif uniq == 'c':
1359            Strm.append(2.0*h*k**3)
1360            Strm.append(2.0*k*h**3)
1361            Strm.append(4.0*h*k*l**2)
1362    else:
1363        Strm.append(h**4)
1364        Strm.append(k**4)
1365        Strm.append(l**4)
1366        Strm.append(3.0*(h*k)**2)
1367        Strm.append(3.0*(h*l)**2)
1368        Strm.append(3.0*(k*l)**2)
1369        Strm.append(2.0*k*h**3)
1370        Strm.append(2.0*h*l**3)
1371        Strm.append(2.0*l*k**3)
1372        Strm.append(2.0*h*k**3)
1373        Strm.append(2.0*l*h**3)
1374        Strm.append(2.0*k*l**3)
1375        Strm.append(4.0*k*l*h**2)
1376        Strm.append(4.0*h*l*k**2)
1377        Strm.append(4.0*k*h*l**2)
1378    return Strm
1379   
1380def Muiso2Shkl(muiso,SGData,cell):
1381    "this is to convert isotropic mustrain to generalized Shkls"
1382    import GSASIIlattice as G2lat
1383    A = G2lat.cell2AB(cell)[0]
1384   
1385    def minMus(Shkl,muiso,H,SGData,A):
1386        U = np.inner(A.T,H)
1387        S = np.array(MustrainCoeff(U,SGData))
1388        Sum = np.sqrt(np.sum(np.multiply(S,Shkl[:,np.newaxis]),axis=0))
1389        rad = np.sqrt(np.sum((Sum[:,np.newaxis]*H)**2,axis=1))
1390        return (muiso-rad)**2
1391       
1392    laue = SGData['SGLaue']
1393    PHI = np.linspace(0.,360.,60,True)
1394    PSI = np.linspace(0.,180.,60,True)
1395    X = np.outer(npsind(PHI),npsind(PSI))
1396    Y = np.outer(npcosd(PHI),npsind(PSI))
1397    Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1398    HKL = np.dstack((X,Y,Z))
1399    if laue in ['m3','m3m']:
1400        S0 = [1000.,1000.]
1401    elif laue in ['6/m','6/mmm','3m1']:
1402        S0 = [1000.,1000.,1000.]
1403    elif laue in ['31m','3']:
1404        S0 = [1000.,1000.,1000.,1000.]
1405    elif laue in ['3R','3mR']:
1406        S0 = [1000.,1000.,1000.,1000.]
1407    elif laue in ['4/m','4/mmm']:
1408        S0 = [1000.,1000.,1000.,1000.]
1409    elif laue in ['mmm']:
1410        S0 = [1000.,1000.,1000.,1000.,1000.,1000.]
1411    elif laue in ['2/m']:
1412        S0 = [1000.,1000.,1000.,0.,0.,0.,0.,0.,0.]
1413    else:
1414        S0 = [1000.,1000.,1000.,1000.,1000., 1000.,1000.,1000.,1000.,1000., 
1415            1000.,1000.,0.,0.,0.]
1416    S0 = np.array(S0)
1417    HKL = np.reshape(HKL,(-1,3))
1418    result = so.leastsq(minMus,S0,(np.ones(HKL.shape[0])*muiso,HKL,SGData,A))
1419    return result[0]
1420       
1421def SytSym(XYZ,SGData):
1422    '''
1423    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
1424
1425    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
1426    :param SGData: from SpcGroup
1427    :Returns: a two element tuple:
1428
1429     * The 1st element is a code for the site symmetry (see GetKNsym)
1430     * The 2nd element is the site multiplicity
1431
1432    '''
1433    def PackRot(SGOps):
1434        IRT = []
1435        for ops in SGOps:
1436            M = ops[0]
1437            irt = 0
1438            for j in range(2,-1,-1):
1439                for k in range(2,-1,-1):
1440                    irt *= 3
1441                    irt += M[k][j]
1442            IRT.append(int(irt))
1443        return IRT
1444       
1445    SymName = ''
1446    Mult = 1
1447    Isym = 0
1448    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1449        Isym = 1073741824
1450    Jdup = 0
1451    Xeqv = GenAtom(XYZ,SGData,True)
1452    IRT = PackRot(SGData['SGOps'])
1453    L = -1
1454    for ic,cen in enumerate(SGData['SGCen']):
1455        for invers in range(int(SGData['SGInv']+1)):
1456            for io,ops in enumerate(SGData['SGOps']):
1457                irtx = (1-2*invers)*IRT[io]
1458                L += 1
1459                if not Xeqv[L][1]:
1460                    Jdup += 1
1461                    jx = GetOprPtrName(str(irtx))
1462                    if jx[2] < 39:
1463                        Isym += 2**(jx[2]-1)
1464    if Isym == 1073741824: Isym = 0
1465    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)/Jdup
1466         
1467    return GetKNsym(str(Isym)),Mult
1468   
1469def ElemPosition(SGData):
1470    ''' Under development.
1471    Object here is to return a list of symmetry element types and locations suitable
1472    for say drawing them.
1473    So far I have the element type... getting all possible locations without lookup may be impossible!
1474    '''
1475    SymElements = []
1476    Inv = SGData['SGInv']
1477    Cen = SGData['SGCen']
1478    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
1479    # get operators & expand if centrosymmetric
1480    Ops = SGData['SGOps']
1481    opM = np.array([op[0].T for op in Ops])
1482    opT = np.array([op[1] for op in Ops])
1483    if Inv:
1484        opM = np.concatenate((opM,-opM))
1485        opT = np.concatenate((opT,-opT))
1486    opMT = zip(opM,opT)
1487    for M,T in opMT[1:]:        #skip I
1488        Dt = int(nl.det(M))
1489        Tr = int(np.trace(M))
1490        Dt = -(Dt-1)/2
1491        Es = eleSym[Tr][Dt]
1492        if Dt:              #rotation-inversion
1493            I = np.eye(3)
1494            if Tr == 1:     #mirrors/glides
1495                if np.any(T):       #glide
1496                    M2 = np.inner(M,M)
1497                    MT = np.inner(M,T)+T
1498                    print 'glide',Es,MT
1499                    print M2
1500                else:               #mirror
1501                    print 'mirror',Es,T
1502                    print I-M
1503                X = [-1,-1,-1]
1504            elif Tr == -3:  # pure inversion
1505                X = np.inner(nl.inv(I-M),T)
1506                print 'inversion',Es,X
1507            else:           #other rotation-inversion
1508                M2 = np.inner(M,M)
1509                MT = np.inner(M,T)+T
1510                print 'rot-inv',Es,MT
1511                print M2
1512                X = [-1,-1,-1]
1513        else:               #rotations
1514            print 'rotation',Es
1515            X = [-1,-1,-1]
1516        #SymElements.append([Es,X])
1517       
1518    return #SymElements
1519   
1520def ApplyStringOps(A,SGData,X,Uij=[]):
1521    'Needs a doc string'
1522    SGOps = SGData['SGOps']
1523    SGCen = SGData['SGCen']
1524    Ax = A.split('+')
1525    Ax[0] = int(Ax[0])
1526    iC = 0
1527    if Ax[0] < 0:
1528        iC = 1
1529    Ax[0] = abs(Ax[0])
1530    nA = Ax[0]%100-1
1531    cA = Ax[0]/100
1532    Cen = SGCen[cA]
1533    M,T = SGOps[nA]
1534    if len(Ax)>1:
1535        cellA = Ax[1].split(',')
1536        cellA = np.array([int(a) for a in cellA])
1537    else:
1538        cellA = np.zeros(3)
1539    newX = (1-2*iC)*(Cen+np.inner(M,X)+T)+cellA
1540    if len(Uij):
1541        U = Uij2U(Uij)
1542        U = np.inner(M,np.inner(U,M).T)
1543        newUij = U2Uij(U)
1544        return [newX,newUij]
1545    else:
1546        return newX
1547       
1548def StringOpsProd(A,B,SGData):
1549    """
1550    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
1551    where '-' indicates inversion, c(>0) is the cell centering operator,
1552    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
1553    Should return resultant string - C. SGData - dictionary using entries:
1554
1555       *  'SGCen': cell centering vectors [0,0,0] at least
1556       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
1557
1558    """
1559    SGOps = SGData['SGOps']
1560    SGCen = SGData['SGCen']
1561    #1st split out the cell translation part & work on the operator parts
1562    Ax = A.split('+'); Bx = B.split('+')
1563    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
1564    iC = 0
1565    if Ax[0]*Bx[0] < 0:
1566        iC = 1
1567    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
1568    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
1569    cA = Ax[0]/100;  cB = Bx[0]/100
1570    Cen = (SGCen[cA]+SGCen[cB])%1.0
1571    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
1572    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
1573    Mc = np.inner(Ma,Mb.T)
1574#    print Ma,Mb,Mc
1575    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
1576#    print Ta,Tb,Tc
1577#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
1578    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
1579    #now the cell translation part
1580    if len(Ax)>1:
1581        cellA = Ax[1].split(',')
1582        cellA = [int(a) for a in cellA]
1583    else:
1584        cellA = [0,0,0]
1585    if len(Bx)>1:
1586        cellB = Bx[1].split(',')
1587        cellB = [int(b) for b in cellB]
1588    else:
1589        cellB = [0,0,0]
1590    cellC = np.add(cellA,cellB)
1591    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
1592        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
1593    return C
1594           
1595def U2Uij(U):
1596    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
1597    return [U[0][0],U[1][1],U[2][2],2.*U[0][1],2.*U[0][2],2.*U[1][2]]
1598   
1599def Uij2U(Uij):
1600    #returns the thermal motion tensor U from Uij as numpy array
1601    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]]])
1602
1603def StandardizeSpcName(spcgroup):
1604    '''Accept a spacegroup name where spaces may have not been used
1605    in the names according to the GSAS convention (spaces between symmetry
1606    for each axis) and return the space group name as used in GSAS
1607    '''
1608    rspc = spcgroup.replace(' ','').upper()
1609    # deal with rhombohedral and hexagonal setting designations
1610    rhomb = ''
1611    if rspc[-1:] == 'R':
1612        rspc = rspc[:-1]
1613        rhomb = ' R'
1614    if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
1615        rspc = rspc[:-1]
1616    # look for a match in the spacegroup lists
1617    for i in spglist.values():
1618        for spc in i:
1619            if rspc == spc.replace(' ','').upper():
1620                return spc + rhomb
1621    # how about the post-2002 orthorhombic names?
1622    for i,spc in sgequiv_2002_orthorhombic:
1623        if rspc == i.replace(' ','').upper():
1624            return spc
1625    # not found
1626    return ''
1627
1628   
1629spglist = {}
1630'''A dictionary of space groups as ordered and named in the pre-2002 International
1631Tables Volume A, except that spaces are used following the GSAS convention to
1632separate the different crystallographic directions.
1633Note that the symmetry codes here will recognize many non-standard space group
1634symbols with different settings. They are ordered by Laue group
1635'''
1636spglist = {
1637    'P1' : ('P 1','P -1',), # 1-2
1638    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
1639        '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
1640    'C2/m':('C 2','C m','C c','C n',
1641        'C 2/m','C 2/c','C 2/n',),
1642    'Pmmm':('P 2 2 2',
1643        'P 2 2 21','P 2 21 2','P 21 2 2',
1644        'P 21 21 2','P 21 2 21','P 2 21 21',
1645        'P 21 21 21',
1646        'P m m 2','P m 2 m','P 2 m m',
1647        'P m c 21','P c m 21','P 21 m a','P 21 a m','P b 21 m','P m 21 b',
1648        'P c c 2','P 2 a a','P b 2 b',
1649        'P m a 2','P b m 2','P 2 m b','P 2 c m','P c 2 m','P m 2 a',
1650        'P c a 21','P b c 21','P 21 a b','P 21 c a','P c 21 b','P b 21 a',
1651        'P n c 2','P c n 2','P 2 n a','P 2 a n','P b 2 n','P n 2 b',
1652        'P m n 21','P n m 21','P 21 m n','P 21 n m','P n 21 m','P m 21 n',
1653        'P b a 2','P 2 c b','P c 2 a',
1654        'P n a 21','P b n 21','P 21 n b','P 21 c n','P c 21 n','P n 21 a',
1655        'P n n 2','P 2 n n','P n 2 n',
1656        'P m m m','P n n n',
1657        'P c c m','P m a a','P b m b',
1658        'P b a n','P n c b','P c n a',
1659        'P m m a','P m m b','P b m m','P c m m','P m c m','P m a m',
1660        'P n n a','P n n b','P b n n','P c n n','P n c n','P n a n',
1661        'P m n a','P n m b','P b m n','P c n m','P n c m','P m a n',
1662        'P c c a','P c c b','P b a a','P c a a','P b c b','P b a b',
1663        'P b a m','P m c b','P c m a',
1664        'P c c n','P n a a','P b n b',
1665        'P b c m','P c a m','P m c a','P m a b','P b m a','P c m b',
1666        'P n n m','P m n n','P n m n',
1667        'P m m n','P n m m','P m n m',
1668        'P b c n','P c a n','P n c a','P n a b','P b n a','P c n b',
1669        'P b c a','P c a b',
1670        'P n m a','P m n b','P b n m','P c m n','P m c n','P n a m',
1671        ),
1672    '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',
1673        'C m 2 a','C 2 m b','C 2 c m','C c 2 m','C 2 c m','C c 2 m', # check: C c 2 m & C c 2 m twice
1674        'C m c a','C m m m','C c c m','C m m a','C c c a','C m c m',),
1675    'Immm':('I 2 2 2','I 21 21 21','I m m m',
1676        'I m m 2','I m 2 m','I 2 m m',
1677        'I b a 2','I 2 c b','I c 2 a',
1678        'I m a 2','I b m 2','I 2 m b','I 2 c m','I c 2 m','I m 2 a',
1679        'I b a m','I m c b','I c m a',
1680        'I b c a','I c a b',
1681        'I m m a','I m m b','I b m m ','I c m m','I m c m','I m a m',),
1682    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
1683        'F m m 2','F m 2 m','F 2 m m',
1684        'F d d 2','F d 2 d','F 2 d d',),
1685    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
1686        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
1687        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
1688        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
1689        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
1690        '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',
1691        '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',
1692        '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',
1693        'P 42/n c m',),
1694    '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',
1695        'I 4 c m','I 41 m d','I 41 c d',
1696        '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',
1697        'I 41/a m d','I 41/a c d'),
1698    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
1699    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
1700        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
1701        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
1702        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
1703        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
1704        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
1705        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
1706    'Pm3m': ('P 2 3','P 21 3','P m 3','P n 3','P a 3','P 4 3 2','P 42 3 2',
1707        'P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n','P m 3 m','P n 3 n',
1708        'P m 3 n','P n 3 m',),
1709    'Im3m':('I 2 3','I 21 3','I m -3','I a -3', 'I 4 3 2','I 41 3 2',
1710        'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a -3 d',),
1711    'Fm3m':('F 2 3','F m -3','F d -3','F 4 3 2','F 41 3 2','F -4 3 m',
1712        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
1713}
1714
1715#'A few non-standard space groups for test use'
1716nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
1717                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
1718                      'R 3 c r','R -3 c r','R -3 m r',),
1719
1720#A list of orthorhombic space groups that were renamed in the 2002 Volume A,
1721# along with the pre-2002 name. The e designates a double glide-plane'''
1722sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
1723                            ('A e a 2', 'A b a 2',),
1724                            ('C m c e', 'C m c a',),
1725                            ('C m m e', 'C m m a',),
1726                            ('C c c e', 'C c c a'),)
1727#Use the space groups types in this order to list the symbols in the
1728#order they are listed in the International Tables, vol. A'''
1729symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
1730               'trigonal', 'hexagonal', 'cubic')
1731
1732# self-test materials follow. Requires files in directory testinp
1733selftestlist = []
1734'''Defines a list of self-tests'''
1735selftestquiet = True
1736def _ReportTest():
1737    'Report name and doc string of current routine when ``selftestquiet`` is False'
1738    if not selftestquiet:
1739        import inspect
1740        caller = inspect.stack()[1][3]
1741        doc = eval(caller).__doc__
1742        if doc is not None:
1743            print('testing '+__file__+' with '+caller+' ('+doc+')')
1744        else:
1745            print('testing '+__file__()+" with "+caller)
1746def test0():
1747    '''self-test #0: exercise MoveToUnitCell'''
1748    _ReportTest()
1749    msg = "MoveToUnitCell failed"
1750    assert (MoveToUnitCell([1,2,3]) == [0,0,0]).all, msg
1751    assert (MoveToUnitCell([2,-1,-2]) == [0,0,0]).all, msg
1752    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9) < 1e-6, msg
1753    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1) < 1e-6, msg
1754selftestlist.append(test0)
1755
1756def test1():
1757    '''self-test #1: SpcGroup against previous results'''
1758    #'''self-test #1: SpcGroup and SGPrint against previous results'''
1759    _ReportTest()
1760    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
1761    if ospath.exists(testdir):
1762        if testdir not in sys.path: sys.path.insert(0,testdir)
1763    import spctestinp
1764    def CompareSpcGroup(spc, referr, refdict, reflist): 
1765        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
1766        # if an error is reported, the dictionary can be ignored
1767        msg0 = "CompareSpcGroup failed on space group %s" % spc
1768        result = SpcGroup(spc)
1769        if result[0] == referr and referr > 0: return True
1770        keys = result[1].keys()
1771        #print result[1]['SpGrp']
1772        #msg = msg0 + " in list lengths"
1773        #assert len(keys) == len(refdict.keys()), msg
1774        for key in refdict.keys():
1775            if key == 'SGOps' or  key == 'SGCen':
1776                msg = msg0 + (" in key %s length" % key)
1777                assert len(refdict[key]) == len(result[1][key]), msg
1778                for i in range(len(refdict[key])):
1779                    msg = msg0 + (" in key %s level 0" % key)
1780                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
1781                    msg = msg0 + (" in key %s level 1" % key)
1782                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
1783            else:
1784                msg = msg0 + (" in key %s" % key)
1785                assert result[1][key] == refdict[key], msg
1786        msg = msg0 + (" in key %s reflist" % key)
1787        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
1788        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
1789        # for now disable SGPrint testing, output has changed
1790        #assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
1791    for spc in spctestinp.SGdat:
1792        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
1793selftestlist.append(test1)
1794
1795def test2():
1796    '''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
1797    _ReportTest()
1798    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
1799    if ospath.exists(testdir):
1800        if testdir not in sys.path: sys.path.insert(0,testdir)
1801    import sgtbxtestinp
1802    def CompareWcctbx(spcname, cctbx_in, debug=0):
1803        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
1804        cctbx = cctbx_in[:] # make copy so we don't delete from the original
1805        spc = (SpcGroup(spcname))[1]
1806        if debug: print spc['SpGrp']
1807        if debug: print spc['SGCen']
1808        latticetype = spcname.strip().upper()[0]
1809        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
1810        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
1811        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
1812        onebar = [1]
1813        if spc['SGInv']: onebar.append(-1)
1814        for (op,off) in spc['SGOps']:
1815            for inv in onebar:
1816                for cen in spc['SGCen']:
1817                    noff = off + cen
1818                    noff = MoveToUnitCell(noff)
1819                    mult = tuple((op*inv).ravel().tolist())
1820                    if debug: print "\n%s: %s + %s" % (spcname,mult,noff)
1821                    for refop in cctbx:
1822                        if debug: print refop
1823                        # check the transform
1824                        if refop[:9] != mult: continue
1825                        if debug: print "mult match"
1826                        # check the translation
1827                        reftrans = list(refop[-3:])
1828                        reftrans = MoveToUnitCell(reftrans)
1829                        if all(abs(noff - reftrans) < 1.e-5):
1830                            cctbx.remove(refop)
1831                            break
1832                    else:
1833                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
1834    for key in sgtbxtestinp.sgtbx:
1835        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
1836selftestlist.append(test2)
1837
1838def test3(): 
1839    '''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
1840     for selected space groups against info in IT Volume A '''
1841    _ReportTest()
1842    def ExerciseSiteSym (spc, crdlist):
1843        'compare site symmetries and multiplicities for a specified space group'
1844        msg = "failed on site sym test for %s" % spc
1845        (E,S) = SpcGroup(spc)
1846        assert not E, msg
1847        for t in crdlist:
1848            symb, m = SytSym(t[0],S)
1849            if symb.strip() != t[2].strip() or m != t[1]:
1850                print spc,t[0],m,symb,t[2]
1851            assert m == t[1]
1852            #assert symb.strip() == t[2].strip()
1853
1854    ExerciseSiteSym('p 1',[
1855            ((0.13,0.22,0.31),1,'1'),
1856            ((0,0,0),1,'1'),
1857            ])
1858    ExerciseSiteSym('p -1',[
1859            ((0.13,0.22,0.31),2,'1'),
1860            ((0,0.5,0),1,'-1'),
1861            ])
1862    ExerciseSiteSym('C 2/c',[
1863            ((0.13,0.22,0.31),8,'1'),
1864            ((0.0,.31,0.25),4,'2(y)'),
1865            ((0.25,.25,0.5),4,'-1'),
1866            ((0,0.5,0),4,'-1'),
1867            ])
1868    ExerciseSiteSym('p 2 2 2',[
1869            ((0.13,0.22,0.31),4,'1'),
1870            ((0,0.5,.31),2,'2(z)'),
1871            ((0.5,.31,0.5),2,'2(y)'),
1872            ((.11,0,0),2,'2(x)'),
1873            ((0,0.5,0),1,'222'),
1874            ])
1875    ExerciseSiteSym('p 4/n',[
1876            ((0.13,0.22,0.31),8,'1'),
1877            ((0.25,0.75,.31),4,'2(z)'),
1878            ((0.5,0.5,0.5),4,'-1'),
1879            ((0,0.5,0),4,'-1'),
1880            ((0.25,0.25,.31),2,'4(001)'),
1881            ((0.25,.75,0.5),2,'-4(001)'),
1882            ((0.25,.75,0.0),2,'-4(001)'),
1883            ])
1884    ExerciseSiteSym('p 31 2 1',[
1885            ((0.13,0.22,0.31),6,'1'),
1886            ((0.13,0.0,0.833333333),3,'2(100)'),
1887            ((0.13,0.13,0.),3,'2(110)'),
1888            ])
1889    ExerciseSiteSym('R 3 c',[
1890            ((0.13,0.22,0.31),18,'1'),
1891            ((0.0,0.0,0.31),6,'3'),
1892            ])
1893    ExerciseSiteSym('R 3 c R',[
1894            ((0.13,0.22,0.31),6,'1'),
1895            ((0.31,0.31,0.31),2,'3(111)'),
1896            ])
1897    ExerciseSiteSym('P 63 m c',[
1898            ((0.13,0.22,0.31),12,'1'),
1899            ((0.11,0.22,0.31),6,'m(100)'),
1900            ((0.333333,0.6666667,0.31),2,'3m(100)'),
1901            ((0,0,0.31),2,'3m(100)'),
1902            ])
1903    ExerciseSiteSym('I a -3',[
1904            ((0.13,0.22,0.31),48,'1'),
1905            ((0.11,0,0.25),24,'2(x)'),
1906            ((0.11,0.11,0.11),16,'3(111)'),
1907            ((0,0,0),8,'-3(111)'),
1908            ])
1909selftestlist.append(test3)
1910
1911if __name__ == '__main__':
1912    # run self-tests
1913    selftestquiet = False
1914    for test in selftestlist:
1915        test()
1916    print "OK"
Note: See TracBrowser for help on using the repository browser.