source: trunk/GSASIIspc.py @ 1397

Last change on this file since 1397 was 1169, checked in by vondreele, 11 years ago

fix space group ops. cif output
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 56.3 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: 2013-12-17 15:38:23 +0000 (Tue, 17 Dec 2013) $
12# $Author: toby $
13# $Revision: 1169 $
14# $URL: trunk/GSASIIspc.py $
15# $Id: GSASIIspc.py 1169 2013-12-17 15:38:23Z toby $
16########### SVN repository information ###################
17import numpy as np
18import numpy.ma as ma
19import numpy.linalg as nl
20import math
21import sys
22import os.path as ospath
23
24import GSASIIpath
25GSASIIpath.SetVersionNumber("$Revision: 1169 $")
26import pyspg
27
28def SpcGroup(SGSymbol):
29    """
30    Determines cell and symmetry information from a short H-M space group name
31
32    :param SGSymbol: space group symbol (string) with spaces between axial fields
33    :returns: (SGError,SGData)
34       * SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
35       * SGData - is a dict (see :ref:`Space Group object<SGData_table>`) with entries:
36       
37             * 'SpGrp': space group symbol, slightly cleaned up
38             * 'Laue':  one of '-1', '2/m', 'mmm', '4/m', '4/mmm', '3R',
39               '3mR', '3', '3m1', '31m', '6/m', '6/mmm', 'm3', 'm3m'
40             * 'SGInv': boolean; True if centrosymmetric, False if not
41             * 'SGLatt': one of 'P', 'A', 'B', 'C', 'I', 'F', 'R'
42             * 'SGUniq': one of 'a', 'b', 'c' if monoclinic, '' otherwise
43             * 'SGCen': cell centering vectors [0,0,0] at least
44             * 'SGOps': symmetry operations as [M,T] so that M*x+T = x'
45             * 'SGSys': one of 'triclinic', 'monoclinic', 'orthorhombic',
46               'tetragonal', 'rhombohedral', 'trigonal', 'hexagonal', 'cubic'
47             * 'SGPolax': one of '', 'x', 'y', 'x y', 'z', 'x z', 'y z',
48               'xyz', '111' for arbitrary axes
49
50    """
51    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
52    LattSym = ('P','A','B','C','I','F','R')
53    UniqSym = ('','','a','b','c','',)
54    SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
55    SGData = {}
56    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
57    SGInfo = pyspg.sgforpy(SGSymbol)
58    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
59    SGData['SGInv'] = bool(SGInfo[1])
60    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
61    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
62    if SGData['SGLatt'] == 'P':
63        SGData['SGCen'] = np.array(([0,0,0],))
64    elif SGData['SGLatt'] == 'A':
65        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
66    elif SGData['SGLatt'] == 'B':
67        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
68    elif SGData['SGLatt'] == 'C':
69        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
70    elif SGData['SGLatt'] == 'I':
71        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
72    elif SGData['SGLatt'] == 'F':
73        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
74    elif SGData['SGLatt'] == 'R':
75        SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
76    SGData['SGOps'] = []
77    for i in range(SGInfo[5]):
78        Mat = np.array(SGInfo[6][i])
79        Trns = np.array(SGInfo[7][i])
80        SGData['SGOps'].append([Mat,Trns])
81    if SGData['SGLaue'] in '-1':
82        SGData['SGSys'] = SysSym[0]
83    elif SGData['SGLaue'] in '2/m':
84        SGData['SGSys'] = SysSym[1]
85    elif SGData['SGLaue'] in 'mmm':
86        SGData['SGSys'] = SysSym[2]
87    elif SGData['SGLaue'] in ['4/m','4/mmm']:
88        SGData['SGSys'] = SysSym[3]
89    elif SGData['SGLaue'] in ['3R','3mR']:
90        SGData['SGSys'] = SysSym[4]
91    elif SGData['SGLaue'] in ['3','3m1','31m']:
92        SGData['SGSys'] = SysSym[5]
93    elif SGData['SGLaue'] in ['6/m','6/mmm']:
94        SGData['SGSys'] = SysSym[6]
95    elif SGData['SGLaue'] in ['m3','m3m']:
96        SGData['SGSys'] = SysSym[7]
97    SGData['SGPolax'] = SGpolar(SGData)
98    return SGInfo[8],SGData
99
100def SGErrors(IErr):
101    '''
102    Interprets the error message code from SpcGroup. Used in SpaceGroup.
103   
104    :param IErr: see SGError in :func:`SpcGroup`
105    :returns:
106        ErrString - a string with the error message or "Unknown error"
107    '''
108
109    ErrString = [' ',
110        'Less than 2 operator fields were found',
111        'Illegal Lattice type, not P, A, B, C, I, F or R',
112        'Rhombohedral lattice requires a 3-axis',
113        'Minus sign does not preceed 1, 2, 3, 4 or 6',
114        'Either a 5-axis anywhere or a 3-axis in field not allowed',
115        ' ',
116        'I for COMPUTED GO TO out of range.',
117        'An a-glide mirror normal to A not allowed',
118        'A b-glide mirror normal to B not allowed',
119        'A c-glide mirror normal to C not allowed',
120        'D-glide in a primitive lattice not allowed',
121        'A 4-axis not allowed in the 2nd operator field',
122        'A 6-axis not allowed in the 2nd operator field',
123        'More than 24 matrices needed to define group',
124        ' ',
125        'Improper construction of a rotation operator',
126        'Mirror following a / not allowed',
127        'A translation conflict between operators',
128        'The 2bar operator is not allowed',
129        '3 fields are legal only in R & m3 cubic groups',
130        'Syntax error. Expected I -4 3 d at this point',
131        ' ',
132        'A or B centered tetragonal not allowed',
133        ' ','unknown error in sgroup',' ',' ',' ',
134        'Illegal character in the space group symbol',
135        ]
136    try:
137        return ErrString[IErr]
138    except:
139        return "Unknown error"
140
141def SGpolar(SGData):
142    '''
143    Determine identity of polar axes if any
144    '''
145    POL = ('','x','y','x y','z','x z','y z','xyz','111')
146    NP = [1,2,4]
147    NPZ = [0,1]
148    for M,T in SGData['SGOps']:
149        for i in range(3):
150            if M[i][i] <= 0.: NP[i] = 0
151        if M[0][2] > 0: NPZ[0] = 8
152        if M[1][2] > 0: NPZ[1] = 0
153    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
154    return POL[NPol]
155   
156def SGPrint(SGData):
157    '''
158    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
159
160    :param SGData: from :func:`SpcGroup`
161    :returns:
162        SGText - list of strings with the space group details
163    '''
164    Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
165    SGText = []
166    SGText.append(' Space Group: '+SGData['SpGrp'])
167    CentStr = 'centrosymmetric'
168    if not SGData['SGInv']:
169        CentStr = 'non'+CentStr
170    if SGData['SGLatt'] in 'ABCIFR':
171        SGText.append(' The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
172    else:
173        SGText.append(' The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())       
174    SGText.append(' Multiplicity of a general site is '+str(Mult))
175    SGText.append(' The Laue symmetry is '+SGData['SGLaue'])
176    if SGData['SGUniq'] in ['a','b','c']:
177        SGText.append(' The unique monoclinic axis is '+SGData['SGUniq'])
178    if SGData['SGInv']:
179        SGText.append(' The inversion center is located at 0,0,0')
180    if SGData['SGPolax']:
181        SGText.append(' The location of the origin is arbitrary in '+SGData['SGPolax'])
182    SGText.append('\n'+' The equivalent positions are:')
183    if SGData['SGLatt'] != 'P':
184        SGText.append('\n ('+Latt2text(SGData['SGLatt'])+')+\n')
185    Ncol = 2
186    line = ' '
187    col = 0
188    for iop,[M,T] in enumerate(SGData['SGOps']):
189        OPtxt = MT2text(M,T)
190        Fld = '(%2i) '%(iop+1)+OPtxt+'\t'
191        line += Fld
192        if '/' not in Fld:
193            line += '\t'
194        col += 1
195        if col == Ncol:
196            SGText.append(line)       
197            line = ' '
198            col = 0
199    SGText.append(line)       
200    return SGText
201
202def AllOps(SGData):
203    '''
204    Returns a list of all operators for a space group, including those for
205    centering and a center of symmetry
206   
207    :param SGData: from :func:`SpcGroup`
208    :returns: (SGTextList,offsetList,symOpList,G2oprList) where
209
210      * SGTextList: a list of strings with formatted and normalized
211        symmetry operators.
212      * offsetList: a tuple of (dx,dy,dz) offsets that relate the GSAS-II
213        symmetry operation to the operator in SGTextList and symOpList.
214        these dx (etc.) values are added to the GSAS-II generated
215        positions to provide the positions that are generated
216        by the normalized symmetry operators.       
217      * symOpList: a list of tuples with the normalized symmetry
218        operations as (M,T) values
219        (see ``SGOps`` in the :ref:`Space Group object<SGData_table>`)
220      * G2oprList: The GSAS-II operations for each symmetry operation as
221        a tuple with (center,mult,opnum), where center is (0,0,0), (0.5,0,0),
222        (0.5,0.5,0.5),...; where mult is 1 or -1 for the center of symmetry
223        and opnum is the number for the symmetry operation, in ``SGOps``
224        (starting with 0).
225    '''
226    SGTextList = []
227    offsetList = []
228    symOpList = []
229    G2oprList = []
230    onebar = (1,)
231    if SGData['SGInv']:
232        onebar += (-1,)
233    for cen in SGData['SGCen']:
234        for mult in onebar:
235            for j,(M,T) in enumerate(SGData['SGOps']):
236                offset = [0,0,0]
237                Tprime = (mult*T)+cen
238                for i in range(3):
239                    while Tprime[i] < 0:
240                        Tprime[i] += 1
241                        offset[i] += 1
242                    while Tprime[i] >= 1:
243                        Tprime[i] += -1
244                        offset[i] += -1
245                OPtxt = MT2text(mult*M,Tprime)
246                SGTextList.append(OPtxt.replace(' ',''))
247                offsetList.append(tuple(offset))
248                symOpList.append((mult*M,Tprime))
249                G2oprList.append((cen,mult,j))
250    return SGTextList,offsetList,symOpList,G2oprList
251   
252def MT2text(M,T):
253    "From space group matrix/translation operator returns text version"
254    XYZ = ('-Z','-Y','-X','X-Y','ERR','Y-X','X','Y','Z')
255    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
256    Fld = ''
257    for j in range(3):
258        IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4))%12
259        IK = int(round(T[j]*12))%12
260        if IK:
261            if IJ < 3:
262                Fld += TRA[IK]+XYZ[IJ]
263            else:
264                Fld += TRA[IK]+'+'+XYZ[IJ]
265        else:
266            Fld += XYZ[IJ]
267        if j != 2: Fld += ', '
268    return Fld
269   
270def Latt2text(Latt):
271    "From lattice type ('P',A', etc.) returns ';' delimited cell centering vectors"
272    lattTxt = {'A':'0,0,0; 0,1/2,1/2','B':'0,0,0; 1/2,0,1/2',
273        'C':'0,0,0; 1/2,1/2,0','I':'0,0,0; 1/2,1/2,1/2',
274        'F':'0,0,0; 0,1/2,1/2; 1/2,0,1/2; 1/2,1/2,0',
275        'R':'0,0,0; 1/3,2/3,2/3; 2/3,1/3,1/3','P':'0,0,0'}
276    return lattTxt[Latt]   
277       
278def SpaceGroup(SGSymbol):
279    '''
280    Print the output of SpcGroup in a nicely formatted way.
281
282    :param SGSymbol: space group symbol (string) with spaces between axial fields
283    :returns: nothing
284    '''
285    E,A = SpcGroup(SGSymbol)
286    if E > 0:
287        print SGErrors(E)
288        return
289    for l in SGPrint(A):
290        print l
291
292def MoveToUnitCell(xyz):
293    '''
294    Translates a set of coordinates so that all values are >=0 and < 1
295
296    :param xyz: a list or numpy array of fractional coordinates
297    :returns: XYZ - numpy array of new coordinates now 0 or greater and less than 1
298    '''
299    XYZ = np.zeros(3)
300    for i,x in enumerate(xyz):
301        XYZ[i] = (x-int(x))%1.0
302    return XYZ
303       
304def Opposite(XYZ,toler=0.0002):
305    '''
306    Gives opposite corner, edge or face of unit cell for position within tolerance.
307        Result may be just outside the cell within tolerance
308
309    :param XYZ: 0 >= np.array[x,y,z] > 1 as by MoveToUnitCell
310    :param toler: unit cell fraction tolerance making opposite
311    :returns:
312        XYZ: array of opposite positions; always contains XYZ
313    '''
314    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]]
315    TB = np.where(abs(XYZ-1)<toler,-1,0)+np.where(abs(XYZ)<toler,1,0)
316    perm = TB*perm3
317    cperm = ['%d%d%d'%(i,j,k) for i,j,k in perm]
318    D = dict(zip(cperm,perm))
319    new = []
320    for key in D:
321        new.append(np.array(D[key])+np.array(XYZ))
322    return new
323       
324def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
325    '''
326    Generates the equivalent positions for a specified coordinate and space group
327
328    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
329    :param SGData: from :func:`SpcGroup`
330    :param All: True return all equivalent positions including duplicates;
331      False return only unique positions
332    :param Uij: [U11,U22,U33,U12,U13,U23] or [] if no Uij
333    :param Move: True move generated atom positions to be inside cell
334      False do not move atoms       
335    :return: [[XYZEquiv],Idup,[UijEquiv]]
336
337      *  [XYZEquiv] is list of equivalent positions (XYZ is first entry)
338      *  Idup = [-][C]SS where SS is the symmetry operator number (1-24), C (if not 0,0,0)
339      * is centering operator number (1-4) and - is for inversion
340        Cell = unit cell translations needed to put new positions inside cell
341        [UijEquiv] - equivalent Uij; absent if no Uij given
342       
343    '''
344    XYZEquiv = []
345    UijEquiv = []
346    Idup = []
347    Cell = []
348    X = np.array(XYZ)
349    if Move:
350        X = MoveToUnitCell(X)
351    for ic,cen in enumerate(SGData['SGCen']):
352        C = np.array(cen)
353        for invers in range(int(SGData['SGInv']+1)):
354            for io,[M,T] in enumerate(SGData['SGOps']):
355                idup = ((io+1)+100*ic)*(1-2*invers)
356                XT = np.inner(M,X)+T
357                if len(Uij):
358                    U = Uij2U(Uij)
359                    U = np.inner(M,np.inner(U,M).T)
360                    newUij = U2Uij(U)
361                if invers:
362                    XT = -XT
363                XT += C
364                newX = MoveToUnitCell(XT)
365                cell = np.asarray(np.rint(newX-XT),dtype=np.int32)
366                if All:
367                    if np.allclose(newX,X,atol=0.0002):
368                        idup = False
369                else:
370                    if True in [np.allclose(newX,oldX,atol=0.0002) for oldX in XYZEquiv]:
371                        idup = False
372                if All or idup:
373                    XYZEquiv.append(newX)
374                    Idup.append(idup)
375                    Cell.append(cell)
376                    if len(Uij):
377                        UijEquiv.append(newUij)                   
378    if len(Uij):
379        return zip(XYZEquiv,UijEquiv,Idup,Cell)
380    else:
381        return zip(XYZEquiv,Idup,Cell)
382
383def GenHKLf(HKL,SGData):
384    '''
385    Uses old GSAS Fortran routine genhkl.for
386
387    :param HKL:  [h,k,l]
388    :param SGData: space group data obtained from SpcGroup
389    :returns: iabsnt,mulp,Uniq,phi
390
391     *   iabsnt = True if reflection is forbidden by symmetry
392     *   mulp = reflection multiplicity including Friedel pairs
393     *   Uniq = numpy array of equivalent hkl in descending order of h,k,l
394
395    '''
396    hklf = HKL+[0,]
397    Ops = SGData['SGOps']
398    OpM = np.array([op[0] for op in Ops])
399    OpT = np.array([op[1] for op in Ops])
400    Inv = SGData['SGInv']
401    Cen = np.array([cen for cen in SGData['SGCen']])
402   
403    Nuniq,Uniq,iabsnt,mulp = pyspg.genhklpy(hklf,len(Ops),OpM,OpT,SGData['SGInv'],len(Cen),Cen)
404    h,k,l,f = Uniq
405    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
406    phi = f[:Nuniq]
407   
408    return iabsnt,mulp,Uniq,phi
409                                 
410def GetOprPtrName(key):
411    'Needs a doc string'
412    OprPtrName = {
413        '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
414        '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
415        '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
416        '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
417        '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
418        '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
419        '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
420        '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
421        '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
422        '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
423        '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
424        '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
425        '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
426        '2223' :[  68,' 3+++2',39],
427        '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
428        '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
429        '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
430        '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
431        '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
432        '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
433        '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
434        '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
435        '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
436        }
437    return OprPtrName[key]
438
439def GetKNsym(key):
440    'Needs a doc string'
441    KNsym = {
442        '0'         :'    1   ','1'         :'   -1   ','64'        :'  2(100)','32'        :'  m(100)',
443        '97'        :'2/m(100)','16'        :'  2(010)','8'         :'  m(010)','25'        :'2/m(010)',
444        '2'         :'  2(001)','4'         :'  m(001)','7'         :'2/m(001)','134217728' :'  2(011)',
445        '67108864'  :'  m(011)','201326593' :'2/m(011)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
446        '3145729'   :'2/m(0+-)','8388608'   :'  2(101)','4194304'   :'  m(101)','12582913'  :'2/m(101)',
447        '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'  2(110)',
448        '512'       :'  m(110)','1537'      :'2/m(110)','256'       :'  2(+-0)','128'       :'  m(+-0)',
449        '385'       :'2/m(+-0)','76'        :'mm2(100)','52'        :'mm2(010)','42'        :'mm2(001)',
450        '135266336' :'mm2(011)','69206048'  :'mm2(0+-)','8650760'   :'mm2(101)','4718600'   :'mm2(+0-)',
451        '1156'      :'mm2(110)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'222(100)',
452        '8912912'   :'222(010)','1282'      :'222(001)','127'       :'  mmm   ','204472417' :'mmm(100)',
453        '13369369'  :'mmm(010)','1927'      :'mmm(001)','33554496'  :'  4(100)','16777280'  :' -4(100)',
454        '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
455        '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
456        '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
457        '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
458        '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
459        '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
460        '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
461        '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
462        '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
463        '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
464        '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
465        '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
466        '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
467        '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
468        '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
469        '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
470        '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
471        '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
472        '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
473        '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
474        '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
475        '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
476        '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
477        '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
478        '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
479        '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   '
480        }
481    return KNsym[key]       
482
483def GetNXUPQsym(siteSym):       
484    'Needs a doc string'
485    NXUPQsym = {
486        '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'  2(100)':(12,18,12,25),'  m(100)':(25,18,12,25),
487        '2/m(100)':( 1,18, 0,-1),'  2(010)':(13,17,13,24),'  m(010)':(24,17,13,24),'2/m(010)':( 1,17, 0,-1),
488        '  2(001)':(14,16,14,23),'  m(001)':(23,16,14,23),'2/m(001)':( 1,16, 0,-1),'  2(011)':(10,23,10,22),
489        '  m(011)':(22,23,10,22),'2/m(011)':( 1,23, 0,-1),'  2(0+-)':(11,24,11,21),'  m(0+-)':(21,24,11,21),
490        '2/m(0+-)':( 1,24, 0,-1),'  2(101)':( 8,21, 8,20),'  m(101)':(20,21, 8,20),'2/m(101)':( 1,21, 0,-1),
491        '  2(+0-)':( 9,22, 9,19),'  m(+0-)':(19,22, 9,19),'2/m(+0-)':( 1,22, 0,-1),'  2(110)':( 6,19, 6,18),
492        '  m(110)':(18,19, 6,18),'2/m(110)':( 1,19, 0,-1),'  2(+-0)':( 7,20, 7,17),'  m(+-0)':(17,20, 7,17),
493        '2/m(+-0)':( 1,20, 0,-1),'mm2(100)':(12,10, 0,-1),'mm2(010)':(13,10, 0,-1),'mm2(001)':(14,10, 0,-1),
494        'mm2(011)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),'mm2(101)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
495        'mm2(110)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'222(100)':( 1,13, 0,-1),
496        '222(010)':( 1,12, 0,-1),'222(001)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'mmm(100)':( 1,13, 0,-1),
497        'mmm(010)':( 1,12, 0,-1),'mmm(001)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
498        '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),
499        '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),
500        '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
501        '  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),
502        '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
503        ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
504        '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
505        '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
506        ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
507        ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
508        '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
509        ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
510        '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
511        '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
512        '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
513        ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
514        '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
515        'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
516        'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
517        'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
518        '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),
519        '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),
520        ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
521        ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
522        ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
523        ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
524        }
525    return NXUPQsym[siteSym]
526
527def GetCSxinel(siteSym): 
528    'Needs a doc string'
529    CSxinel = [[],                         # 0th empty - indices are Fortran style
530        [[0,0,0],[ 0.0, 0.0, 0.0]],      #  0  0  0
531        [[1,1,1],[ 1.0, 1.0, 1.0]],      #  X  X  X
532        [[1,1,1],[ 1.0, 1.0,-1.0]],      #  X  X -X
533        [[1,1,1],[ 1.0,-1.0, 1.0]],      #  X -X  X
534        [[1,1,1],[ 1.0,-1.0,-1.0]],      # -X  X  X
535        [[1,1,0],[ 1.0, 1.0, 0.0]],      #  X  X  0
536        [[1,1,0],[ 1.0,-1.0, 0.0]],      #  X -X  0
537        [[1,0,1],[ 1.0, 0.0, 1.0]],      #  X  0  X
538        [[1,0,1],[ 1.0, 0.0,-1.0]],      #  X  0 -X
539        [[0,1,1],[ 0.0, 1.0, 1.0]],      #  0  Y  Y
540        [[0,1,1],[ 0.0, 1.0,-1.0]],      #  0  Y -Y
541        [[1,0,0],[ 1.0, 0.0, 0.0]],      #  X  0  0
542        [[0,1,0],[ 0.0, 1.0, 0.0]],      #  0  Y  0
543        [[0,0,1],[ 0.0, 0.0, 1.0]],      #  0  0  Z
544        [[1,1,0],[ 1.0, 2.0, 0.0]],      #  X 2X  0
545        [[1,1,0],[ 2.0, 1.0, 0.0]],      # 2X  X  0
546        [[1,1,2],[ 1.0, 1.0, 1.0]],      #  X  X  Z
547        [[1,1,2],[ 1.0,-1.0, 1.0]],      #  X -X  Z
548        [[1,2,1],[ 1.0, 1.0, 1.0]],      #  X  Y  X
549        [[1,2,1],[ 1.0, 1.0,-1.0]],      #  X  Y -X
550        [[1,2,2],[ 1.0, 1.0, 1.0]],      #  X  Y  Y
551        [[1,2,2],[ 1.0, 1.0,-1.0]],      #  X  Y -Y
552        [[1,2,0],[ 1.0, 1.0, 0.0]],      #  X  Y  0
553        [[1,0,2],[ 1.0, 0.0, 1.0]],      #  X  0  Z
554        [[0,1,2],[ 0.0, 1.0, 1.0]],      #  0  Y  Z
555        [[1,1,2],[ 1.0, 2.0, 1.0]],      #  X 2X  Z
556        [[1,1,2],[ 2.0, 1.0, 1.0]],      # 2X  X  Z
557        [[1,2,3],[ 1.0, 1.0, 1.0]],      #  X  Y  Z
558        ]
559    indx = GetNXUPQsym(siteSym)
560    return CSxinel[indx[0]]
561   
562def GetCSuinel(siteSym):
563    "returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers"
564    CSuinel = [[],                                             # 0th empty - indices are Fortran style
565        [[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]],    #  A  A  A  0  0  0
566        [[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]],    #  A  A  C  0  0  0
567        [[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]],    #  A  B  A  0  0  0
568        [[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]],    #  A  B  B  0  0  0
569        [[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]],    #  A  A  A  D  D  D
570        [[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]],    #  A  A  A  D -D -D
571        [[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]],    #  A  A  A  D -D  D
572        [[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]],    #  A  A  A  D  D -D
573        [[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]],    #  A  A  C A/2 0  0
574        [[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]],    #  A  B  C  0  0  0
575        [[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]],    #  A  A  C  D  0  0
576        [[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]],    #  A  B  A  0  E  0
577        [[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]],    #  A  B  B  0  0  F
578        [[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]],    #  A  B  C B/2 0  0
579        [[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]],    #  A  B  C A/2 0  0
580        [[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]],    #  A  B  C  D  0  0
581        [[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]],    #  A  B  C  0  E  0
582        [[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]],    #  A  B  C  0  0  F
583        [[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]],    #  A  A  C  D  E -E
584        [[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]],    #  A  A  C  D  E  E
585        [[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]],    #  A  B  A  D  E -D
586        [[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]],    #  A  B  A  D  E  D
587        [[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]],    #  A  B  B  D -D  F
588        [[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]],    #  A  B  B  D  D  F
589        [[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]],    #  A  B  C B/2 F/2 F
590        [[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]],    #  A  B  C A/2  0  F
591        [[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]],    #  A  B  C B/2  E  0
592        [[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]],    #  A  B  C A/2  E E/2
593        [[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]],    #  A  B  C  D  E   F
594        ]
595    indx = GetNXUPQsym(siteSym)
596    return CSuinel[indx[1]]
597   
598def MustrainNames(SGData):
599    'Needs a doc string'
600    laue = SGData['SGLaue']
601    uniq = SGData['SGUniq']
602    if laue in ['m3','m3m']:
603        return ['S400','S220']
604    elif laue in ['6/m','6/mmm','3m1']:
605        return ['S400','S004','S202']
606    elif laue in ['31m','3']:
607        return ['S400','S004','S202','S211']
608    elif laue in ['3R','3mR']:
609        return ['S400','S220','S310','S211']
610    elif laue in ['4/m','4/mmm']:
611        return ['S400','S004','S220','S022']
612    elif laue in ['mmm']:
613        return ['S400','S040','S004','S220','S202','S022']
614    elif laue in ['2/m']:
615        SHKL = ['S400','S040','S004','S220','S202','S022']
616        if uniq == 'a':
617            SHKL += ['S013','S031','S211']
618        elif uniq == 'b':
619            SHKL += ['S301','S103','S121']
620        elif uniq == 'c':
621            SHKL += ['S130','S310','S112']
622        return SHKL
623    else:
624        SHKL = ['S400','S040','S004','S220','S202','S022']
625        SHKL += ['S310','S103','S031','S130','S301','S013']
626        SHKL += ['S211','S121','S112']
627        return SHKL
628
629def HStrainNames(SGData):
630    'Needs a doc string'
631    laue = SGData['SGLaue']
632    uniq = SGData['SGUniq']
633    if laue in ['m3','m3m']:
634        return ['D11','eA']         #add cubic strain term
635    elif laue in ['6/m','6/mmm','3m1','31m','3']:
636        return ['D11','D33']
637    elif laue in ['3R','3mR']:
638        return ['D11','D12']
639    elif laue in ['4/m','4/mmm']:
640        return ['D11','D33']
641    elif laue in ['mmm']:
642        return ['D11','D22','D33']
643    elif laue in ['2/m']:
644        Dij = ['D11','D22','D33']
645        if uniq == 'a':
646            Dij += ['D23']
647        elif uniq == 'b':
648            Dij += ['D13']
649        elif uniq == 'c':
650            Dij += ['D12']
651        return Dij
652    else:
653        Dij = ['D11','D22','D33','D12','D13','D23']
654        return Dij
655   
656def MustrainCoeff(HKL,SGData):
657    'Needs a doc string'
658    #NB: order of terms is the same as returned by MustrainNames
659    laue = SGData['SGLaue']
660    uniq = SGData['SGUniq']
661    h,k,l = HKL
662    Strm = []
663    if laue in ['m3','m3m']:
664        Strm.append(h**4+k**4+l**4)
665        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
666    elif laue in ['6/m','6/mmm','3m1']:
667        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
668        Strm.append(l**4)
669        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
670    elif laue in ['31m','3']:
671        Strm.append(h**4+k**4+2.0*k*h**3+2.0*h*k**3+3.0*(h*k)**2)
672        Strm.append(l**4)
673        Strm.append(3.0*((h*l)**2+(k*l)**2+h*k*l**2))
674        Strm.append(4.0*h*k*l*(h+k))
675    elif laue in ['3R','3mR']:
676        Strm.append(h**4+k**4+l**4)
677        Strm.append(3.0*((h*k)**2+(h*l)**2+(k*l)**2))
678        Strm.append(2.0*(h*l**3+l*k**3+k*h**3)+2.0*(l*h**3+k*l**3+l*k**3))
679        Strm.append(4.0*(k*l*h**2+h*l*k**2+h*k*l**2))
680    elif laue in ['4/m','4/mmm']:
681        Strm.append(h**4+k**4)
682        Strm.append(l**4)
683        Strm.append(3.0*(h*k)**2)
684        Strm.append(3.0*((h*l)**2+(k*l)**2))
685    elif laue in ['mmm']:
686        Strm.append(h**4)
687        Strm.append(k**4)
688        Strm.append(l**4)
689        Strm.append(3.0*(h*k)**2)
690        Strm.append(3.0*(h*l)**2)
691        Strm.append(3.0*(k*l)**2)
692    elif laue in ['2/m']:
693        Strm.append(h**4)
694        Strm.append(k**4)
695        Strm.append(l**4)
696        Strm.append(3.0*(h*k)**2)
697        Strm.append(3.0*(h*l)**2)
698        Strm.append(3.0*(k*l)**2)
699        if uniq == 'a':
700            Strm.append(2.0*k*l**3)
701            Strm.append(2.0*l*k**3)
702            Strm.append(4.0*k*l*h**2)
703        elif uniq == 'b':
704            Strm.append(2.0*l*h**3)
705            Strm.append(2.0*h*l**3)
706            Strm.append(4.0*h*l*k**2)
707        elif uniq == 'c':
708            Strm.append(2.0*h*k**3)
709            Strm.append(2.0*k*h**3)
710            Strm.append(4.0*h*k*l**2)
711    else:
712        Strm.append(h**4)
713        Strm.append(k**4)
714        Strm.append(l**4)
715        Strm.append(3.0*(h*k)**2)
716        Strm.append(3.0*(h*l)**2)
717        Strm.append(3.0*(k*l)**2)
718        Strm.append(2.0*k*h**3)
719        Strm.append(2.0*h*l**3)
720        Strm.append(2.0*l*k**3)
721        Strm.append(2.0*h*k**3)
722        Strm.append(2.0*l*h**3)
723        Strm.append(2.0*k*l**3)
724        Strm.append(4.0*k*l*h**2)
725        Strm.append(4.0*h*l*k**2)
726        Strm.append(4.0*k*h*l**2)
727    return Strm
728   
729def Muiso2Shkl(muiso,SGData,cell):
730    "this is to convert isotropic mustrain to generalized Shkls - doesn't work just now"
731    import GSASIIlattice as G2lat
732    from scipy.optimize import fmin
733    A = G2lat.cell2AB(cell)[0]
734    def minMus(Shkl,H,muiso,SGData,A):
735        U = np.inner(A.T,H)
736        S = np.array(MustrainCoeff(H.T,SGData))
737        sum = np.sqrt(np.sum(np.multiply(S,Shkl)))
738        return abs(muiso-sum*H)
739    laue = SGData['SGLaue']
740    if laue in ['m3','m3m']:
741        H = [[1,0,0],[1,1,0]]
742        S0 = [0.01,0.01]
743    elif laue in ['6/m','6/mmm','3m1']:
744        H = [[1,0,0],[0,0,1],[1,0,1]]
745        S0 = [0.01,0.01,0.01]
746    elif laue in ['31m','3']:
747        H = [[1,0,0],[0,0,1],[1,0,1],[1,1,1]]
748        S0 = [0.01,0.01,0.01,0.01]
749    elif laue in ['3R','3mR']:
750        H = [[1,0,0],[1,1,0],[1,0,1],[1,1,1]]
751        S0 = [0.01,0.01,0.01,0.01]
752    elif laue in ['4/m','4/mmm']:
753        H = [[1,0,0],[0,0,1],[1,1,0],[1,0,1]]
754        S0 = [0.01,0.01,0.01,0.01]
755    elif laue in ['mmm']:
756        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1]]
757        S0 = [0.01,0.01,0.01,0.01,0.01,0.01]
758    elif laue in ['2/m']:
759        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1]]
760        if uniq == 'a':
761            H.append([0,1,-1])
762            H.append([0,-2,1])
763        elif uniq == 'b':
764            H.append([1,0,-1])
765            H.append([-2,0,1])
766        elif uniq == 'c':
767            H.append([1,-1,0])
768            H.append([-2,1,0])
769        H.append([1,1,1])
770        S0 = [9*[0.01,]]
771    else:
772        H = [[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],
773            [-1,1,0],[1,0,-1],[0,-1,1],[1,-2,0],[-2,0,1],[0,1,-2],
774            [1,-1,1],[-1, 1, 1],[1,-1,1]]
775        S0 = [15*[0.01,]]
776    H = np.array(H)
777    S0 = np.array(S0)
778    return fmin(minMus,S0,(H,muiso,SGData,A))
779       
780def SytSym(XYZ,SGData):
781    '''
782    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
783
784    :param XYZ: an array, tuple or list containing 3 elements: x, y & z
785    :param SGData: from SpcGroup
786    :Returns: a two element tuple:
787
788     * The 1st element is a code for the site symmetry (see GetKNsym)
789     * The 2nd element is the site multiplicity
790
791    '''
792    def PackRot(SGOps):
793        IRT = []
794        for ops in SGOps:
795            M = ops[0]
796            irt = 0
797            for j in range(2,-1,-1):
798                for k in range(2,-1,-1):
799                    irt *= 3
800                    irt += M[k][j]
801            IRT.append(int(irt))
802        return IRT
803       
804    SymName = ''
805    Mult = 1
806    Isym = 0
807    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
808        Isym = 1073741824
809    Jdup = 0
810    Xeqv = GenAtom(XYZ,SGData,True)
811    IRT = PackRot(SGData['SGOps'])
812    L = -1
813    for ic,cen in enumerate(SGData['SGCen']):
814        for invers in range(int(SGData['SGInv']+1)):
815            for io,ops in enumerate(SGData['SGOps']):
816                irtx = (1-2*invers)*IRT[io]
817                L += 1
818                if not Xeqv[L][1]:
819                    Jdup += 1
820                    jx = GetOprPtrName(str(irtx))
821                    if jx[2] < 39:
822                        Isym += 2**(jx[2]-1)
823    if Isym == 1073741824: Isym = 0
824    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)/Jdup
825         
826    return GetKNsym(str(Isym)),Mult
827   
828def ElemPosition(SGData):
829    ''' Under development.
830    Object here is to return a list of symmetry element types and locations suitable
831    for say drawing them.
832    So far I have the element type... getting all possible locations without lookup may be impossible!
833    '''
834    SymElements = []
835    Inv = SGData['SGInv']
836    Cen = SGData['SGCen']
837    eleSym = {-3:['','-1'],-2:['',-6],-1:['2','-4'],0:['3','-3'],1:['4','m'],2:['6',''],3:['1','']}
838    # get operators & expand if centrosymmetric
839    Ops = SGData['SGOps']
840    opM = np.array([op[0].T for op in Ops])
841    opT = np.array([op[1] for op in Ops])
842    if Inv:
843        opM = np.concatenate((opM,-opM))
844        opT = np.concatenate((opT,-opT))
845    opMT = zip(opM,opT)
846    for M,T in opMT[1:]:        #skip I
847        Dt = int(nl.det(M))
848        Tr = int(np.trace(M))
849        Dt = -(Dt-1)/2
850        Es = eleSym[Tr][Dt]
851        if Dt:              #rotation-inversion
852            I = np.eye(3)
853            if Tr == 1:     #mirrors/glides
854                if np.any(T):       #glide
855                    M2 = np.inner(M,M)
856                    MT = np.inner(M,T)+T
857                    print 'glide',Es,MT
858                    print M2
859                else:               #mirror
860                    print 'mirror',Es,T
861                    print I-M
862                X = [-1,-1,-1]
863            elif Tr == -3:  # pure inversion
864                X = np.inner(nl.inv(I-M),T)
865                print 'inversion',Es,X
866            else:           #other rotation-inversion
867                M2 = np.inner(M,M)
868                MT = np.inner(M,T)+T
869                print 'rot-inv',Es,MT
870                print M2
871                X = [-1,-1,-1]
872               
873           
874           
875        else:               #rotations
876            print 'rotation',Es
877            X = [-1,-1,-1]
878        #SymElements.append([Es,X])
879       
880    return #SymElements
881   
882def ApplyStringOps(A,SGData,X,Uij=[]):
883    'Needs a doc string'
884    SGOps = SGData['SGOps']
885    SGCen = SGData['SGCen']
886    Ax = A.split('+')
887    Ax[0] = int(Ax[0])
888    iC = 0
889    if Ax[0] < 0:
890        iC = 1
891    Ax[0] = abs(Ax[0])
892    nA = Ax[0]%100-1
893    cA = Ax[0]/100
894    Cen = SGCen[cA]
895    M,T = SGOps[nA]
896    if len(Ax)>1:
897        cellA = Ax[1].split(',')
898        cellA = np.array([int(a) for a in cellA])
899    else:
900        cellA = np.zeros(3)
901    newX = (1-2*iC)*(Cen+np.inner(M,X)+T)+cellA
902    if len(Uij):
903        U = Uij2U(Uij)
904        U = np.inner(M,np.inner(U,M).T)
905        newUij = U2Uij(U)
906        return [newX,newUij]
907    else:
908        return newX
909       
910def StringOpsProd(A,B,SGData):
911    """
912    Find A*B where A & B are in strings '-' + '100*c+n' + '+ijk'
913    where '-' indicates inversion, c(>0) is the cell centering operator,
914    n is operator number from SgOps and ijk are unit cell translations (each may be <0).
915    Should return resultant string - C. SGData - dictionary using entries:
916
917       *  'SGCen': cell centering vectors [0,0,0] at least
918       *  'SGOps': symmetry operations as [M,T] so that M*x+T = x'
919
920    """
921    SGOps = SGData['SGOps']
922    SGCen = SGData['SGCen']
923    #1st split out the cell translation part & work on the operator parts
924    Ax = A.split('+'); Bx = B.split('+')
925    Ax[0] = int(Ax[0]); Bx[0] = int(Bx[0])
926    iC = 0
927    if Ax[0]*Bx[0] < 0:
928        iC = 1
929    Ax[0] = abs(Ax[0]); Bx[0] = abs(Bx[0])
930    nA = Ax[0]%100-1;  nB = Bx[0]%100-1
931    cA = Ax[0]/100;  cB = Bx[0]/100
932    Cen = (SGCen[cA]+SGCen[cB])%1.0
933    cC = np.nonzero([np.allclose(C,Cen) for C in SGCen])[0][0]
934    Ma,Ta = SGOps[nA]; Mb,Tb = SGOps[nB]
935    Mc = np.inner(Ma,Mb.T)
936#    print Ma,Mb,Mc
937    Tc = (np.add(np.inner(Mb,Ta)+1.,Tb))%1.0
938#    print Ta,Tb,Tc
939#    print [np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps]
940    nC = np.nonzero([np.allclose(M,Mc)&np.allclose(T,Tc) for M,T in SGOps])[0][0]
941    #now the cell translation part
942    if len(Ax)>1:
943        cellA = Ax[1].split(',')
944        cellA = [int(a) for a in cellA]
945    else:
946        cellA = [0,0,0]
947    if len(Bx)>1:
948        cellB = Bx[1].split(',')
949        cellB = [int(b) for b in cellB]
950    else:
951        cellB = [0,0,0]
952    cellC = np.add(cellA,cellB)
953    C = str(((nC+1)+(100*cC))*(1-2*iC))+'+'+ \
954        str(int(cellC[0]))+','+str(int(cellC[1]))+','+str(int(cellC[2]))
955    return C
956           
957def U2Uij(U):
958    #returns the UIJ vector U11,U22,U33,U12,U13,U23 from tensor U
959    return [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
960   
961def Uij2U(Uij):
962    #returns the thermal motion tensor U from Uij as numpy array
963    return np.array([[Uij[0],Uij[3],Uij[4]],[Uij[3],Uij[1],Uij[5]],[Uij[4],Uij[5],Uij[2]]])
964
965def StandardizeSpcName(spcgroup):
966    '''Accept a spacegroup name where spaces may have not been used
967    in the names according to the GSAS convention (spaces between symmetry
968    for each axis) and return the space group name as used in GSAS
969    '''
970    rspc = spcgroup.replace(' ','').upper()
971    # deal with rhombohedral and hexagonal setting designations
972    rhomb = ''
973    if rspc[-1:] == 'R':
974        rspc = rspc[:-1]
975        rhomb = ' R'
976    if rspc[-1:] == 'H': # hexagonal is assumed and thus can be ignored
977        rspc = rspc[:-1]
978    # look for a match in the spacegroup lists
979    for i in spglist.values():
980        for spc in i:
981            if rspc == spc.replace(' ','').upper():
982                return spc + rhomb
983    # how about the post-2002 orthorhombic names?
984    for i,spc in sgequiv_2002_orthorhombic:
985        if rspc == i.replace(' ','').upper():
986            return spc
987    # not found
988    return ''
989
990   
991'''A dictionary of space groups as ordered and named in the pre-2002 International
992Tables Volume A, except that spaces are used following the GSAS convention to
993separate the different crystallographic directions.
994Note that the symmetry codes here will recognize many non-standard space group
995symbols with different settings. They are ordered by Laue group
996'''
997spglist = {
998    'P1' : ('P 1','P -1',), # 1-2
999    'P2/m': ('P 2','P 21','P m','P a','P c','P n',
1000        '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
1001    'C2/m':('C 2','C m','C c','C n',
1002        'C 2/m','C 2/c','C 2/n',),
1003    'Pmmm':('P 2 2 2',
1004        'P 2 2 21','P 2 21 2','P 21 2 2',
1005        'P 21 21 2','P 21 2 21','P 2 21 21',
1006        'P 21 21 21',
1007        'P m m 2','P m 2 m','P 2 m m',
1008        'P m c 21','P c m 21','P 21 m a','P 21 a m','P b 21 m','P m 21 b',
1009        'P c c 2','P 2 a a','P b 2 b',
1010        'P m a 2','P b m 2','P 2 m b','P 2 c m','P c 2 m','P m 2 a',
1011        'P c a 21','P b c 21','P 21 a b','P 21 c a','P c 21 b','P b 21 a',
1012        'P n c 2','P c n 2','P 2 n a','P 2 a n','P b 2 n','P n 2 b',
1013        'P m n 21','P n m 21','P 21 m n','P 21 n m','P n 21 m','P m 21 n',
1014        'P b a 2','P 2 c b','P c 2 a',
1015        'P n a 21','P b n 21','P 21 n b','P 21 c n','P c 21 n','P n 21 a',
1016        'P n n 2','P 2 n n','P n 2 n',
1017        'P m m m','P n n n',
1018        'P c c m','P m a a','P b m b',
1019        'P b a n','P n c b','P c n a',
1020        'P m m a','P m m b','P b m m','P c m m','P m c m','P m a m',
1021        'P n n a','P n n b','P b n n','P c n n','P n c n','P n a n',
1022        'P m n a','P n m b','P b m n','P c n m','P n c m','P m a n',
1023        'P c c a','P c c b','P b a a','P c a a','P b c b','P b a b',
1024        'P b a m','P m c b','P c m a',
1025        'P c c n','P n a a','P b n b',
1026        'P b c m','P c a m','P m c a','P m a b','P b m a','P c m b',
1027        'P n n m','P m n n','P n m n',
1028        'P m m n','P n m m','P m n m',
1029        'P b c n','P c a n','P n c a','P n a b','P b n a','P c n b',
1030        'P b c a','P c a b',
1031        'P n m a','P m n b','P b n m','P c m n','P m c n','P n a m',
1032        ),
1033    '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',
1034        '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
1035        'C m c a','C m m m','C c c m','C m m a','C c c a','C m c m',),
1036    'Immm':('I 2 2 2','I 21 21 21','I m m m',
1037        'I m m 2','I m 2 m','I 2 m m',
1038        'I b a 2','I 2 c b','I c 2 a',
1039        'I m a 2','I b m 2','I 2 m b','I 2 c m','I c 2 m','I m 2 a',
1040        'I b a m','I m c b','I c m a',
1041        'I b c a','I c a b',
1042        'I m m a','I m m b','I b m m ','I c m m','I m c m','I m a m',),
1043    'Fmmm':('F 2 2 2','F m m m', 'F d d d',
1044        'F m m 2','F m 2 m','F 2 m m',
1045        'F d d 2','F d 2 d','F 2 d d',),
1046    'P4/mmm':('P 4','P 41','P 42','P 43','P -4','P 4/m','P 42/m','P 4/n','P 42/n',
1047        'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
1048        'P 42 21 2','P 43 2 2','P 43 21 2','P 4 m m','P 4 b m','P 42 c m',
1049        'P 42 n m','P 4 c c','P 4 n c','P 42 m c','P 42 b c','P -4 2 m',
1050        'P -4 2 c','P -4 21 m','P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2',
1051        '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',
1052        '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',
1053        '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',
1054        'P 42/n c m',),
1055    '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',
1056        'I 4 c m','I 41 m d','I 41 c d',
1057        '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',
1058        'I 41/a m d','I 41/a c d'),
1059    'R3-H':('R 3','R -3','R 3 2','R 3 m','R 3 c','R -3 m','R -3 c',),
1060    'P6/mmm': ('P 3','P 31','P 32','P -3','P 3 1 2','P 3 2 1','P 31 1 2',
1061        'P 31 2 1','P 32 1 2','P 32 2 1', 'P 3 m 1','P 3 1 m','P 3 c 1',
1062        'P 3 1 c','P -3 1 m','P -3 1 c','P -3 m 1','P -3 c 1','P 6','P 61',
1063        'P 65','P 62','P 64','P 63','P -6','P 6/m','P 63/m','P 6 2 2',
1064        'P 61 2 2','P 65 2 2','P 62 2 2','P 64 2 2','P 63 2 2','P 6 m m',
1065        'P 6 c c','P 63 c m','P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m',
1066        'P -6 2 c','P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',),
1067    'Pm3m': ('P 2 3','P 21 3','P m 3','P n 3','P a 3','P 4 3 2','P 42 3 2',
1068        'P 43 3 2','P 41 3 2','P -4 3 m','P -4 3 n','P m 3 m','P n 3 n',
1069        'P m 3 n','P n 3 m',),
1070    'Im3m':('I 2 3','I 21 3','I m -3','I a -3', 'I 4 3 2','I 41 3 2',
1071        'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a -3 d',),
1072    'Fm3m':('F 2 3','F m -3','F d -3','F 4 3 2','F 41 3 2','F -4 3 m',
1073        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
1074}
1075
1076#'A few non-standard space groups for test use'
1077nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
1078                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
1079                      'R 3 c r','R -3 c r','R -3 m r',),
1080
1081#A list of orthorhombic space groups that were renamed in the 2002 Volume A,
1082# along with the pre-2002 name. The e designates a double glide-plane'''
1083sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
1084                            ('A e a 2', 'A b a 2',),
1085                            ('C m c e', 'C m c a',),
1086                            ('C m m e', 'C m m a',),
1087                            ('C c c e', 'C c c a'),)
1088#Use the space groups types in this order to list the symbols in the
1089#order they are listed in the International Tables, vol. A'''
1090symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
1091               'trigonal', 'hexagonal', 'cubic')
1092
1093# self-test materials follow. Requires files in directory testinp
1094selftestlist = []
1095'''Defines a list of self-tests'''
1096selftestquiet = True
1097def _ReportTest():
1098    'Report name and doc string of current routine when ``selftestquiet`` is False'
1099    if not selftestquiet:
1100        import inspect
1101        caller = inspect.stack()[1][3]
1102        doc = eval(caller).__doc__
1103        if doc is not None:
1104            print('testing '+__file__+' with '+caller+' ('+doc+')')
1105        else:
1106            print('testing '+__file__()+" with "+caller)
1107def test0():
1108    '''self-test #0: exercise MoveToUnitCell'''
1109    _ReportTest()
1110    msg = "MoveToUnitCell failed"
1111    assert (MoveToUnitCell([1,2,3]) == [0,0,0]).all, msg
1112    assert (MoveToUnitCell([2,-1,-2]) == [0,0,0]).all, msg
1113    assert abs(MoveToUnitCell(np.array([-.1]))[0]-0.9) < 1e-6, msg
1114    assert abs(MoveToUnitCell(np.array([.1]))[0]-0.1) < 1e-6, msg
1115selftestlist.append(test0)
1116
1117def test1():
1118    '''self-test #1: SpcGroup and SGPrint against previous results'''
1119    _ReportTest()
1120    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
1121    if ospath.exists(testdir):
1122        if testdir not in sys.path: sys.path.insert(0,testdir)
1123    import spctestinp
1124    def CompareSpcGroup(spc, referr, refdict, reflist): 
1125        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
1126        # if an error is reported, the dictionary can be ignored
1127        msg0 = "CompareSpcGroup failed on space group %s" % spc
1128        result = SpcGroup(spc)
1129        if result[0] == referr and referr > 0: return True
1130        keys = result[1].keys()
1131        #print result[1]['SpGrp']
1132        msg = msg0 + " in list lengths"
1133        assert len(keys) == len(refdict.keys()), msg
1134        for key in keys:
1135            if key == 'SGOps' or  key == 'SGCen':
1136                msg = msg0 + (" in key %s length" % key)
1137                assert len(refdict[key]) == len(result[1][key]), msg
1138                for i in range(len(refdict[key])):
1139                    msg = msg0 + (" in key %s level 0" % key)
1140                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
1141                    msg = msg0 + (" in key %s level 1" % key)
1142                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
1143            else:
1144                msg = msg0 + (" in key %s" % key)
1145                assert result[1][key] == refdict[key], msg
1146        msg = msg0 + (" in key %s reflist" % key)
1147        #for (l1,l2) in zip(reflist, SGPrint(result[1])):
1148        #    assert l2.replace('\t','').replace(' ','') == l1.replace(' ',''), 'SGPrint ' +msg
1149        assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
1150    for spc in spctestinp.SGdat:
1151        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
1152selftestlist.append(test1)
1153
1154def test2():
1155    '''self-test #2: SpcGroup against cctbx (sgtbx) computations'''
1156    _ReportTest()
1157    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
1158    if ospath.exists(testdir):
1159        if testdir not in sys.path: sys.path.insert(0,testdir)
1160    import sgtbxtestinp
1161    def CompareWcctbx(spcname, cctbx_in, debug=0):
1162        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
1163        cctbx = cctbx_in[:] # make copy so we don't delete from the original
1164        spc = (SpcGroup(spcname))[1]
1165        if debug: print spc['SpGrp']
1166        if debug: print spc['SGCen']
1167        latticetype = spcname.strip().upper()[0]
1168        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
1169        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
1170        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
1171        onebar = [1]
1172        if spc['SGInv']: onebar.append(-1)
1173        for (op,off) in spc['SGOps']:
1174            for inv in onebar:
1175                for cen in spc['SGCen']:
1176                    noff = off + cen
1177                    noff = MoveToUnitCell(noff)
1178                    mult = tuple((op*inv).ravel().tolist())
1179                    if debug: print "\n%s: %s + %s" % (spcname,mult,noff)
1180                    for refop in cctbx:
1181                        if debug: print refop
1182                        # check the transform
1183                        if refop[:9] != mult: continue
1184                        if debug: print "mult match"
1185                        # check the translation
1186                        reftrans = list(refop[-3:])
1187                        reftrans = MoveToUnitCell(reftrans)
1188                        if all(abs(noff - reftrans) < 1.e-5):
1189                            cctbx.remove(refop)
1190                            break
1191                    else:
1192                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
1193    for key in sgtbxtestinp.sgtbx:
1194        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
1195selftestlist.append(test2)
1196
1197def test3(): 
1198    '''self-test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
1199     for selected space groups against info in IT Volume A '''
1200    _ReportTest()
1201    def ExerciseSiteSym (spc, crdlist):
1202        'compare site symmetries and multiplicities for a specified space group'
1203        msg = "failed on site sym test for %s" % spc
1204        (E,S) = SpcGroup(spc)
1205        assert not E, msg
1206        for t in crdlist:
1207            symb, m = SytSym(t[0],S)
1208            if symb.strip() != t[2].strip() or m != t[1]:
1209                print spc,t[0],m,symb
1210            assert m == t[1]
1211            #assert symb.strip() == t[2].strip()
1212
1213    ExerciseSiteSym('p 1',[
1214            ((0.13,0.22,0.31),1,'1'),
1215            ((0,0,0),1,'1'),
1216            ])
1217    ExerciseSiteSym('p -1',[
1218            ((0.13,0.22,0.31),2,'1'),
1219            ((0,0.5,0),1,'-1'),
1220            ])
1221    ExerciseSiteSym('C 2/c',[
1222            ((0.13,0.22,0.31),8,'1'),
1223            ((0.0,.31,0.25),4,'2(010)'),
1224            ((0.25,.25,0.5),4,'-1'),
1225            ((0,0.5,0),4,'-1'),
1226            ])
1227    ExerciseSiteSym('p 2 2 2',[
1228            ((0.13,0.22,0.31),4,'1'),
1229            ((0,0.5,.31),2,'2(001)'),
1230            ((0.5,.31,0.5),2,'2(010)'),
1231            ((.11,0,0),2,'2(100)'),
1232            ((0,0.5,0),1,'222'),
1233            ])
1234    ExerciseSiteSym('p 4/n',[
1235            ((0.13,0.22,0.31),8,'1'),
1236            ((0.25,0.75,.31),4,'2(001)'),
1237            ((0.5,0.5,0.5),4,'-1'),
1238            ((0,0.5,0),4,'-1'),
1239            ((0.25,0.25,.31),2,'4(001)'),
1240            ((0.25,.75,0.5),2,'-4(001)'),
1241            ((0.25,.75,0.0),2,'-4(001)'),
1242            ])
1243    ExerciseSiteSym('p 31 2 1',[
1244            ((0.13,0.22,0.31),6,'1'),
1245            ((0.13,0.0,0.833333333),3,'2(100)'),
1246            ((0.13,0.13,0.),3,'2(110)'),
1247            ])
1248    ExerciseSiteSym('R 3 c',[
1249            ((0.13,0.22,0.31),18,'1'),
1250            ((0.0,0.0,0.31),6,'3'),
1251            ])
1252    ExerciseSiteSym('R 3 c R',[
1253            ((0.13,0.22,0.31),6,'1'),
1254            ((0.31,0.31,0.31),2,'3(111)'),
1255            ])
1256    ExerciseSiteSym('P 63 m c',[
1257            ((0.13,0.22,0.31),12,'1'),
1258            ((0.11,0.22,0.31),6,'m(100)'),
1259            ((0.333333,0.6666667,0.31),2,'3m(100)'),
1260            ((0,0,0.31),2,'3m(100)'),
1261            ])
1262    ExerciseSiteSym('I a -3',[
1263            ((0.13,0.22,0.31),48,'1'),
1264            ((0.11,0,0.25),24,'2(100)'),
1265            ((0.11,0.11,0.11),16,'3(111)'),
1266            ((0,0,0),8,'-3(111)'),
1267            ])
1268selftestlist.append(test3)
1269
1270if __name__ == '__main__':
1271    # run self-tests
1272    selftestquiet = False
1273    for test in selftestlist:
1274        test()
1275    print "OK"
Note: See TracBrowser for help on using the repository browser.