source: trunk/GSASIIspc.py @ 956

Last change on this file since 956 was 956, checked in by toby, 10 years ago

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