source: trunk/GSASIIspc.py @ 1471

Last change on this file since 1471 was 1471, checked in by vondreele, 9 years ago

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