source: trunk/GSASIIspc.py @ 1539

Last change on this file since 1539 was 1539, checked in by vondreele, 8 years ago

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