source: trunk/GSASIIspc.py @ 973

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