source: trunk/GSASIIspc.py @ 1530

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

supersymmetry modifications & revise space group display/print to give better formatted columns of operators

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