source: trunk/GSASIIspc.py @ 1489

Last change on this file since 1489 was 1489, checked in by vondreele, 7 years ago

correct Uij to U6 conversions Uij = Uk/2 for i != j; k > 2 in G2lattice & G2spc
This fixes derivative issues for aniso tropic atoms in Al2O3
change names of sytsyms e.g. 2(100) to 2(x) for non-hex/trig space groups. Removes possible ambiguity with ortho, etc. cases.
Still checking extinction issues

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