source: trunk/GSASIIspc.py @ 142

Last change on this file since 142 was 142, checked in by vondreel, 13 years ago

major refactor of phase editing, EXP & PDB import
replace "general" list with dictionary,
make Uij obey site symmetry for user access
File size: 38.5 KB
Line 
1"GSASII - Space group interpretion routines"
2
3import numpy as np
4import sys
5import os.path as ospath
6
7import GSASIIpath
8import pyspg
9
10def SpcGroup(SGSymbol):
11    '''
12   Determines cell and symmetry information from a short H-M space group name
13   input: space group symbol (string) with spaces between axial fields
14   returns [SGError,SGData]
15       SGError = 0 for no errors; >0 for errors (see SGErrors below for details)
16       returns dictionary SGData with entries:
17         'SpGrp': space group symbol slightly cleaned up
18         'Laue':  one of '-1','2/m','mmm','4/m','4/mmm','3R','3mR','3',
19                  '3m1','31m','6/m','6/mmm','m3','m3m'
20         'SGInv': boolean; True if centrosymmetric, False if not
21         'SGLatt': one of 'P','A','B','C','I','F','R'
22         'SGUniq': one of 'a','b','c' if monoclinic, '' otherwise
23         'SGCen': cell centering vectors [0,0,0] at least
24         'SGOps': symmetry operations as [M,T] so that M*x+T = x'
25         'SGSys': one of 'triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic'
26       '''
27    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
28    LattSym = ('P','A','B','C','I','F','R')
29    UniqSym = ('','','a','b','c','',)
30    SysSym = ('triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic')
31    SGData = {}
32    SGData['SpGrp'] = SGSymbol.strip().lower().capitalize()
33    SGInfo = pyspg.sgforpy(SGSymbol)
34    SGData['SGLaue'] = LaueSym[SGInfo[0]-1]
35    SGData['SGInv'] = bool(SGInfo[1])
36    SGData['SGLatt'] = LattSym[SGInfo[2]-1]
37    SGData['SGUniq'] = UniqSym[SGInfo[3]+1]
38    if SGData['SGLatt'] == 'P':
39        SGData['SGCen'] = np.array(([0,0,0],))
40    elif SGData['SGLatt'] == 'A':
41        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5]))
42    elif SGData['SGLatt'] == 'B':
43        SGData['SGCen'] = np.array(([0,0,0],[.5,0,.5]))
44    elif SGData['SGLatt'] == 'C':
45        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,0,]))
46    elif SGData['SGLatt'] == 'I':
47        SGData['SGCen'] = np.array(([0,0,0],[.5,.5,.5]))
48    elif SGData['SGLatt'] == 'F':
49        SGData['SGCen'] = np.array(([0,0,0],[0,.5,.5],[.5,0,.5],[.5,.5,0,]))
50    elif SGData['SGLatt'] == 'R':
51        SGData['SGCen'] = np.array(([0,0,0],[1./3.,2./3.,2./3.],[2./3.,1./3.,1./3.]))
52    SGData['SGOps'] = []
53    for i in range(SGInfo[5]):
54        Mat = np.array(SGInfo[6][i])
55        Trns = np.array(SGInfo[7][i])
56        SGData['SGOps'].append([Mat,Trns])
57    if SGData['SGLaue'] in '-1':
58        SGData['SGSys'] = SysSym[0]
59    elif SGData['SGLaue'] in '2/m':
60        SGData['SGSys'] = SysSym[1]
61    elif SGData['SGLaue'] in 'mmm':
62        SGData['SGSys'] = SysSym[2]
63    elif SGData['SGLaue'] in ['4/m','4/mmm']:
64        SGData['SGSys'] = SysSym[3]
65    elif SGData['SGLaue'] in ['3R','3mR']:
66        SGData['SGSys'] = SysSym[4]
67    elif SGData['SGLaue'] in ['3','3m1','31m']:
68        SGData['SGSys'] = SysSym[5]
69    elif SGData['SGLaue'] in ['6/m','6/mmm']:
70        SGData['SGSys'] = SysSym[6]
71    elif SGData['SGLaue'] in ['m3','m3m']:
72        SGData['SGSys'] = SysSym[7]
73    return SGInfo[8],SGData
74
75def SGErrors(IErr):
76    '''Interprets the error message code from SpcGroup. Used in SpaceGroup.
77    input:  SGError, from SpcGroup
78    returns a string with the error message or "Unknown error"
79    '''
80
81    ErrString = [' ',
82        'Less than 2 operator fields were found',
83        'Illegal Lattice type, not P, A, B, C, I, F or R',
84        'Rhombohedral lattice requires a 3-axis',
85        'Minus sign does not preceed 1, 2, 3, 4 or 6',
86        'Either a 5-axis anywhere or a 3-axis in field not allowed',
87        ' ',
88        'I for COMPUTED GO TO out of range.',
89        'An a-glide mirror normal to A not allowed',
90        'A b-glide mirror normal to B not allowed',
91        'A c-glide mirror normal to C not allowed',
92        'D-glide in a primitive lattice not allowed',
93        'A 4-axis not allowed in the 2nd operator field',
94        'A 6-axis not allowed in the 2nd operator field',
95        'More than 24 matrices needed to define group',
96        ' ',
97        'Improper construction of a rotation operator',
98        'Mirror following a / not allowed',
99        'A translation conflict between operators',
100        'The 2bar operator is not allowed',
101        '3 fields are legal only in R & m3 cubic groups',
102        'Syntax error. Expected I -4 3 d at this point',
103        ' ',
104        'A or B centered tetragonal not allowed',
105        ' ','unknown error in sgroup',' ',' ',' ',
106        'Illegal character in the space group symbol',
107        ]
108    try:
109        return ErrString[IErr]
110    except:
111        return "Unknown error"
112   
113def SGPrint(SGData):
114    '''
115    Print the output of SpcGroup in a nicely formatted way. Used in SpaceGroup
116    input:  SGData, from SpcGroup
117    returns a list of strings with the space group details
118    '''
119    XYZ = ('-Z ','-Y ','-X ','X-Y','ERR','Y-X',' X ',' Y ',' Z ','+X ','+Y ','+Z ')
120    TRA = ('   ','ERR','1/6','1/4','1/3','ERR','1/2','ERR','2/3','3/4','5/6','ERR')
121    POL = (' ','x','y','x y','z','x z','y z','xyz','111')
122    Mult = len(SGData['SGCen'])*len(SGData['SGOps'])*(int(SGData['SGInv'])+1)
123    NP = [1,2,4]
124    NPZ = [0,1]
125    for M,T in SGData['SGOps']:
126        for i in range(3):
127            if M[i][i] <= 0.: NP[i] = 0
128        if M[0][2] > 0: NPZ[0] = 8
129        if M[1][2] > 0: NPZ[1] = 0
130    NPol = (NP[0]+NP[1]+NP[2]+NPZ[0]*NPZ[1])*(1-int(SGData['SGInv']))
131    SGText = []
132    SGText.append('Space Group '+SGData['SpGrp'])
133    CentStr = 'centrosymmetric'
134    if not SGData['SGInv']:
135        CentStr = 'non'+CentStr
136    if SGData['SGLatt'] in 'ABCIFR':
137        SGText.append('The lattice is '+CentStr+' '+SGData['SGLatt']+'-centered '+SGData['SGSys'].lower())
138    else:
139        SGText.append('The lattice is '+CentStr+' '+'primitive '+SGData['SGSys'].lower())       
140    SGText.append('Multiplicity of a general site is '+str(Mult))
141    SGText.append('The Laue symmetry is '+SGData['SGLaue'])
142    if SGData['SGUniq'] in ['a','b','c']:
143        SGText.append('The unique monoclinic axis is '+SGData['SGUniq'])
144    if SGData['SGInv']:
145        SGText.append('The inversion center is located at 0,0,0')
146    if NPol:
147        SGText.append('The location of the origin is arbitrary in '+POL[NPol])
148    SGText.append('\n'+'The equivalent positions are:')
149    if SGData['SGLatt'] in 'A':
150        SGText.append('\n'+'    (0,0,0; 0,1/2,1/2)+')
151    elif SGData['SGLatt'] in 'B':
152        SGText.append('\n'+'    (0,0,0; 1/2,0,1/2)+')
153    elif SGData['SGLatt'] in 'C':
154        SGText.append('\n'+'    (0,0,0; 1/2,1/2,0)+')
155    elif SGData['SGLatt'] in 'I':
156        SGText.append('\n'+'    (0,0,0; 1/2,1/2,1/2)+')
157    elif SGData['SGLatt'] in 'F':
158        SGText.append('\n'+'    (0,0,0; 0,1/2,1/2; 1/2,0,1/2; 1/2,1/2,0)+')
159    elif SGData['SGLatt'] in 'R':
160        SGText.append('\n'+'    (0,0,0; 1/3,2/3,2/3; 2/3,1/3,1/3)+')
161    if SGData['SGLaue'] in ['-1','2/m','mmm','4/m','4/mmm']:
162        Ncol = 2
163    else:
164        Ncol = 3
165    line = ''
166    for iop,[M,T] in enumerate(SGData['SGOps']):
167        if iop % Ncol == 0:
168            SGText.append(line)       
169            line = ''
170        Fld = '(%2i) ' % (iop+1)
171        for j in range(3):
172            IJ = int(round(2*M[j][0]+3*M[j][1]+4*M[j][2]+4)) % 12
173            IK = int(round(T[j]*12)) % 12
174            if IK > 0 and IJ > 4: IJ += 3
175            Fld += TRA[IK]+XYZ[IJ]
176            if j != 2: Fld += ','
177        line += Fld
178    SGText.append(line)       
179    return SGText
180   
181def SpaceGroup(SgSym):
182    '''
183    Print the output of SpcGroup in a nicely formatted way.
184      input: space group symbol (string) with spaces between axial fields
185      returns nothing
186    '''
187    E,A = SpcGroup(SgSym)
188    if E > 0:
189        print SGErrors(E)
190        return
191    for l in SGPrint(A):
192        print l
193
194def MoveToUnitCell(XYZ):
195    '''
196    Translates a set of coordinates so that all values are >=0 and < 1
197      input: a list or numpy array of any length. Note that the object is modified  in place.
198      output: none
199    '''
200    for i,x in enumerate(XYZ):
201        x = ((x % 1.0)+1.0) % 1.0
202        if x > 0.9999: x = 0.0
203        XYZ[i] = x
204       
205def GenAtom(XYZ,SGData,ifAll=False):
206    '''
207    Generates the equivalent positions for a specified coordinate and space group
208    input: 
209       XYZ an array, tuple or list containing 3 elements: x, y & z
210       SGData, from SpcGroup
211       ifAll=True causes the return to provide the unique set of
212                  equivalent positions
213            =False causes the input position to be repeated. This is the default,
214                   but why someone would want this, I am not sure.
215    Returns a list of two element tuples:
216       The first element is the coordinate as a three-element array and
217       the second describes the symmetry used to generate the site, of form [-][C]SS
218          C indicates a centering operation was used (omitted if the 1st, [0,0,0])
219          SS is the symmetry operator number (1-24)
220          - indicates the center of symmetry was used (omitted otherwise)     
221    '''
222    XYZEquiv = []
223    Idup = []
224    X = np.array(XYZ)
225    MoveToUnitCell(X)
226    XYZEquiv.append(np.array(X))
227    Idup.append(1)
228    for ic,cen in enumerate(SGData['SGCen']):
229        C = np.array(cen)
230        for invers in range(int(SGData['SGInv']+1)):
231            for io,ops in enumerate(SGData['SGOps']):
232                idup = ((io+1)+100*ic)*(1-2*invers)
233                T = np.array(ops[1])
234                M =  np.array(ops[0])
235                newX = np.sum(M*X,axis=1)+T
236                if invers:
237                    newX = -newX
238                newX += C
239                MoveToUnitCell(newX)
240                New = True
241                if ifAll:
242                    if np.allclose(newX,X,atol=0.0002):
243                        New = False
244                        idup = 0                   
245                    XYZEquiv.append(newX)
246                else:
247                    for oldX in XYZEquiv[:-1]:
248                        if np.allclose(newX,oldX,atol=0.0002):
249                            New = False
250                            idup = 0
251                    if New or ifAll:
252                        XYZEquiv.append(newX)
253                if ifAll and len(XYZEquiv) == 2:
254                    Idup.append(1)
255                else:
256                    Idup.append(idup)
257    if ifAll:
258        return zip(XYZEquiv[1:],Idup[1:])                  #eliminate duplicate initial entry
259    else:
260        return zip(XYZEquiv,Idup)
261                                   
262def GetOprPtrName(key):           
263    OprPtrName = {
264        '-6643':[   2,' 1bar ', 1],'6479' :[  10,'  2z  ', 2],'-6479':[   9,'  mz  ', 3],
265        '6481' :[   7,'  my  ', 4],'-6481':[   6,'  2y  ', 5],'6641' :[   4,'  mx  ', 6],
266        '-6641':[   3,'  2x  ', 7],'6591' :[  28,' m+-0 ', 8],'-6591':[  27,' 2+-0 ', 9],
267        '6531' :[  25,' m110 ',10],'-6531':[  24,' 2110 ',11],'6537' :[  61,'  4z  ',12],
268        '-6537':[  62,' -4z  ',13],'975'  :[  68,' 3+++1',14],'6456' :[ 114,'  3z1 ',15],
269        '-489' :[  73,' 3+-- ',16],'483'  :[  78,' 3-+- ',17],'-969' :[  83,' 3--+ ',18],
270        '819'  :[  22,' m+0- ',19],'-819' :[  21,' 2+0- ',20],'2431' :[  16,' m0+- ',21],
271        '-2431':[  15,' 20+- ',22],'-657' :[  19,' m101 ',23],'657'  :[  18,' 2101 ',24],
272        '1943' :[  48,' -4x  ',25],'-1943':[  47,'  4x  ',26],'-2429':[  13,' m011 ',27],
273        '2429' :[  12,' 2011 ',28],'639'  :[  55,' -4y  ',29],'-639' :[  54,'  4y  ',30],
274        '-6484':[ 146,' 2010 ', 4],'6484' :[ 139,' m010 ', 5],'-6668':[ 145,' 2100 ', 6],
275        '6668' :[ 138,' m100 ', 7],'-6454':[ 148,' 2120 ',18],'6454' :[ 141,' m120 ',19],
276        '-6638':[ 149,' 2210 ',20],'6638' :[ 142,' m210 ',21],              #search ends here
277        '2223' :[  68,' 3+++2',39],
278        '6538' :[ 106,'  6z1 ',40],'-2169':[  83,' 3--+2',41],'2151' :[  73,' 3+--2',42],
279        '2205' :[  79,'-3-+-2',43],'-2205':[  78,' 3-+-2',44],'489'  :[  74,'-3+--1',45],
280        '801'  :[  53,'  4y1 ',46],'1945' :[  47,'  4x3 ',47],'-6585':[  62,' -4z3 ',48],
281        '6585' :[  61,'  4z3 ',49],'6584' :[ 114,'  3z2 ',50],'6666' :[ 106,'  6z5 ',51],
282        '6643' :[   1,' Iden ',52],'-801' :[  55,' -4y1 ',53],'-1945':[  48,' -4x3 ',54],
283        '-6666':[ 105,' -6z5 ',55],'-6538':[ 105,' -6z1 ',56],'-2223':[  69,'-3+++2',57],
284        '-975' :[  69,'-3+++1',58],'-6456':[ 113,' -3z1 ',59],'-483' :[  79,'-3-+-1',60],
285        '969'  :[  84,'-3--+1',61],'-6584':[ 113,' -3z2 ',62],'2169' :[  84,'-3--+2',63],
286        '-2151':[  74,'-3+--2',64],'0':[0,' ????',0]
287        }
288    return OprPtrName[key]
289
290def GetKNsym(key):
291    KNsym = {
292        '0'         :'    1   ','1'         :'   -1   ','64'        :'  2(100)','32'        :'  m(100)',
293        '97'        :'2/m(100)','16'        :'  2(010)','8'         :'  m(010)','25'        :'2/m(010)',
294        '2'         :'  2(001)','4'         :'  m(001)','7'         :'2/m(001)','134217728' :'  2(011)',
295        '67108864'  :'  m(011)','201326593' :'2/m(011)','2097152'   :'  2(0+-)','1048576'   :'  m(0+-)',
296        '3145729'   :'2/m(0+-)','8388608'   :'  2(101)','4194304'   :'  m(101)','12582913'  :'2/m(101)',
297        '524288'    :'  2(+0-)','262144'    :'  m(+0-)','796433'    :'2/m(+0-)','1024'      :'  2(110)',
298        '512'       :'  m(110)','1537'      :'2/m(110)','256'       :'  2(+-0)','128'       :'  m(+-0)',
299        '385'       :'2/m(+-0)','76'        :'mm2(100)','52'        :'mm2(010)','42'        :'mm2(001)',
300        '135266336' :'mm2(011)','69206048'  :'mm2(0+-)','8650760'   :'mm2(101)','4718600'   :'mm2(+0-)',
301        '1156'      :'mm2(110)','772'       :'mm2(+-0)','82'        :'  222   ','136314944' :'222(100)',
302        '8912912'   :'222(010)','1282'      :'222(001)','127'       :'  mmm   ','204472417' :'mmm(100)',
303        '13369369'  :'mmm(010)','1927'      :'mmm(001)','33554496'  :'  4(100)','16777280'  :' -4(100)',
304        '50331745'  :'4/m(100)','169869394' :'422(100)','84934738'  :'-42m 100','101711948' :'4mm(100)',
305        '254804095' :'4/mmm100','536870928 ':'  4(010)','268435472' :' -4(010)','805306393' :'4/m (10)',
306        '545783890' :'422(010)','272891986' :'-42m 010','541327412' :'4mm(010)','818675839' :'4/mmm010',
307        '2050'      :'  4(001)','4098'      :' -4(001)','6151'      :'4/m(001)','3410'      :'422(001)',
308        '4818'      :'-42m 001','2730'      :'4mm(001)','8191'      :'4/mmm001','8192'      :'  3(111)',
309        '8193'      :' -3(111)','2629888'   :' 32(111)','1319040'   :' 3m(111)','3940737'   :'-3m(111)',
310        '32768'     :'  3(+--)','32769'     :' -3(+--)','10519552'  :' 32(+--)','5276160'   :' 3m(+--)',
311        '15762945'  :'-3m(+--)','65536'     :'  3(-+-)','65537'     :' -3(-+-)','134808576' :' 32(-+-)',
312        '67437056'  :' 3m(-+-)','202180097' :'-3m(-+-)','131072'    :'  3(--+)','131073'    :' -3(--+)',
313        '142737664' :' 32(--+)','71434368'  :' 3m(--+)','214040961' :'-3m(--+)','237650'    :'   23   ',
314        '237695'    :'   m3   ','715894098' :'   432  ','358068946' :'  -43m  ','1073725439':'   m3m  ',
315        '68157504'  :' mm2d100','4456464'   :' mm2d010','642'       :' mm2d001','153092172' :'-4m2 100',
316        '277348404' :'-4m2 010','5418'      :'-4m2 001','1075726335':'  6/mmm ','1074414420':'-6m2 100',
317        '1075070124':'-6m2 120','1075069650':'   6mm  ','1074414890':'   622  ','1073758215':'   6/m  ',
318        '1073758212':'   -6   ','1073758210':'    6   ','1073759865':'-3m(100)','1075724673':'-3m(120)',
319        '1073758800':' 3m(100)','1075069056':' 3m(120)','1073759272':' 32(100)','1074413824':' 32(120)',
320        '1073758209':'   -3   ','1073758208':'    3   ','1074135143':'mmm(100)','1075314719':'mmm(010)',
321        '1073743751':'mmm(110)','1074004034':' mm2z100','1074790418':' mm2z010','1073742466':' mm2z110',
322        '1074004004':'mm2(100)','1074790412':'mm2(010)','1073742980':'mm2(110)','1073872964':'mm2(120)',
323        '1074266132':'mm2(210)','1073742596':'mm2(+-0)','1073872930':'222(100)','1074266122':'222(010)',
324        '1073743106':'222(110)','1073741831':'2/m(001)','1073741921':'2/m(100)','1073741849':'2/m(010)',
325        '1073743361':'2/m(110)','1074135041':'2/m(120)','1075314689':'2/m(210)','1073742209':'2/m(+-0)',
326        '1073741828':' m(001) ','1073741888':' m(100) ','1073741840':' m(010) ','1073742336':' m(110) ',
327        '1074003968':' m(120) ','1074790400':' m(210) ','1073741952':' m(+-0) ','1073741826':' 2(001) ',
328        '1073741856':' 2(100) ','1073741832':' 2(010) ','1073742848':' 2(110) ','1073872896':' 2(120) ',
329        '1074266112':' 2(210) ','1073742080':' 2(+-0) ','1073741825':'   -1   '
330        }
331    return KNsym[key]       
332
333def GetNXUPQsym(siteSym):       
334    NXUPQsym = {
335        '    1   ':(28,29,28,28),'   -1   ':( 1,29,28, 0),'  2(100)':(12,18,12,25),'  m(100)':(25,18,12,25),
336        '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),
337        '  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),
338        '  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),
339        '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),
340        '  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),
341        '  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),
342        '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),
343        'mm2(011)':(10,13, 0,-1),'mm2(0+-)':(11,13, 0,-1),'mm2(101)':( 8,12, 0,-1),'mm2(+0-)':( 9,12, 0,-1),
344        'mm2(110)':( 6,11, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'  222   ':( 1,10, 0,-1),'222(100)':( 1,13, 0,-1),
345        '222(010)':( 1,12, 0,-1),'222(001)':( 1,11, 0,-1),'  mmm   ':( 1,10, 0,-1),'mmm(100)':( 1,13, 0,-1),
346        'mmm(010)':( 1,12, 0,-1),'mmm(001)':( 1,11, 0,-1),'  4(100)':(12, 4,12, 0),' -4(100)':( 1, 4,12, 0),
347        '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),
348        '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),
349        '422(010)':( 1, 3, 0,-1),'-42m 010':( 1, 3, 0,-1),'4mm(010)':(13, 3, 0,-1),'4/mmm010':(1, 3, 0,-1,),
350        '  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),
351        '-42m 001':( 1, 2, 0,-1),'4mm(001)':(14, 2, 0,-1),'4/mmm001':( 1, 2, 0,-1),'  3(111)':( 2, 5, 2, 0),
352        ' -3(111)':( 1, 5, 2, 0),' 32(111)':( 1, 5, 0, 2),' 3m(111)':( 2, 5, 0, 2),'-3m(111)':( 1, 5, 0,-1),
353        '  3(+--)':( 5, 8, 5, 0),' -3(+--)':( 1, 8, 5, 0),' 32(+--)':( 1, 8, 0, 5),' 3m(+--)':( 5, 8, 0, 5),
354        '-3m(+--)':( 1, 8, 0,-1),'  3(-+-)':( 4, 7, 4, 0),' -3(-+-)':( 1, 7, 4, 0),' 32(-+-)':( 1, 7, 0, 4),
355        ' 3m(-+-)':( 4, 7, 0, 4),'-3m(-+-)':( 1, 7, 0,-1),'  3(--+)':( 3, 6, 3, 0),' -3(--+)':( 1, 6, 3, 0),
356        ' 32(--+)':( 1, 6, 0, 3),' 3m(--+)':( 3, 6, 0, 3),'-3m(--+)':( 1, 6, 0,-1),'   23   ':( 1, 1, 0, 0),
357        '   m3   ':( 1, 1, 0, 0),'   432  ':( 1, 1, 0, 0),'  -43m  ':( 1, 1, 0, 0),'   m3m  ':( 1, 1, 0, 0),
358        ' mm2d100':(12,13, 0,-1),' mm2d010':(13,12, 0,-1),' mm2d001':(14,11, 0,-1),'-4m2 100':( 1, 4, 0,-1),
359        '-4m2 010':( 1, 3, 0,-1),'-4m2 001':( 1, 2, 0,-1),'  6/mmm ':( 1, 9, 0,-1),'-6m2 100':( 1, 9, 0,-1),
360        '-6m2 120':( 1, 9, 0,-1),'   6mm  ':(14, 9, 0,-1),'   622  ':( 1, 9, 0,-1),'   6/m  ':( 1, 9,14,-1),
361        '   -6   ':( 1, 9,14, 0),'    6   ':(14, 9,14, 0),'-3m(100)':( 1, 9, 0,-1),'-3m(120)':( 1, 9, 0,-1),
362        ' 3m(100)':(14, 9, 0,14),' 3m(120)':(14, 9, 0,14),' 32(100)':( 1, 9, 0,14),' 32(120)':( 1, 9, 0,14),
363        '   -3   ':( 1, 9,14, 0),'    3   ':(14, 9,14, 0),'mmm(100)':( 1,14, 0,-1),'mmm(010)':( 1,15, 0,-1),
364        'mmm(110)':( 1,11, 0,-1),' mm2z100':(14,14, 0,-1),' mm2z010':(14,15, 0,-1),' mm2z110':(14,11, 0,-1),
365        'mm2(100)':(12,14, 0,-1),'mm2(010)':(13,15, 0,-1),'mm2(110)':( 6,11, 0,-1),'mm2(120)':(15,14, 0,-1),
366        'mm2(210)':(16,15, 0,-1),'mm2(+-0)':( 7,11, 0,-1),'222(100)':( 1,14, 0,-1),'222(010)':( 1,15, 0,-1),
367        '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),
368        '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),
369        ' m(001) ':(23,16,14,23),' m(100) ':(26,25,12,26),' m(010) ':(27,28,13,27),' m(110) ':(18,19, 6,18),
370        ' m(120) ':(24,27,15,24),' m(210) ':(25,26,16,25),' m(+-0) ':(17,20, 7,17),' 2(001) ':(14,16,14,23),
371        ' 2(100) ':(12,25,12,26),' 2(010) ':(13,28,13,27),' 2(110) ':( 6,19, 6,18),' 2(120) ':(15,27,15,24),
372        ' 2(210) ':(16,26,16,25),' 2(+-0) ':( 7,20, 7,17),'   -1   ':( 1,29,28, 0)
373        }
374    return NXUPQsym[siteSym]
375
376def GetCSxinel(siteSym): 
377    CSxinel = [[],                         # 0th empty - indices are Fortran style
378        [[0,0,0],[ 0.0, 0.0, 0.0]],      #  0  0  0
379        [[1,1,1],[ 1.0, 1.0, 1.0]],      #  X  X  X
380        [[1,1,1],[ 1.0, 1.0,-1.0]],      #  X  X -X
381        [[1,1,1],[ 1.0,-1.0, 1.0]],      #  X -X  X
382        [[1,1,1],[ 1.0,-1.0,-1.0]],      # -X  X  X
383        [[1,1,0],[ 1.0, 1.0, 0.0]],      #  X  X  0
384        [[1,1,0],[ 1.0,-1.0, 0.0]],      #  X -X  0
385        [[1,0,1],[ 1.0, 0.0, 1.0]],      #  X  0  X
386        [[1,0,1],[ 1.0, 0.0,-1.0]],      #  X  0 -X
387        [[0,1,1],[ 0.0, 1.0, 1.0]],      #  0  Y  Y
388        [[0,1,1],[ 0.0, 1.0,-1.0]],      #  0  Y -Y
389        [[1,0,0],[ 1.0, 0.0, 0.0]],      #  X  0  0
390        [[0,1,0],[ 0.0, 1.0, 0.0]],      #  0  Y  0
391        [[0,0,1],[ 0.0, 0.0, 1.0]],      #  0  0  Z
392        [[1,1,0],[ 1.0, 2.0, 0.0]],      #  X 2X  0
393        [[1,1,0],[ 2.0, 1.0, 0.0]],      # 2X  X  0
394        [[1,1,2],[ 1.0, 1.0, 1.0]],      #  X  X  Z
395        [[1,1,2],[ 1.0,-1.0, 1.0]],      #  X -X  Z
396        [[1,2,1],[ 1.0, 1.0, 1.0]],      #  X  Y  X
397        [[1,2,1],[ 1.0, 1.0,-1.0]],      #  X  Y -X
398        [[1,2,2],[ 1.0, 1.0, 1.0]],      #  X  Y  Y
399        [[1,2,2],[ 1.0, 1.0,-1.0]],      #  X  Y -Y
400        [[1,2,0],[ 1.0, 1.0, 0.0]],      #  X  Y  0
401        [[1,0,2],[ 1.0, 0.0, 1.0]],      #  X  0  Z
402        [[0,1,2],[ 0.0, 1.0, 1.0]],      #  0  Y  Z
403        [[1,1,2],[ 1.0, 2.0, 1.0]],      #  X 2X  Z
404        [[1,1,2],[ 2.0, 1.0, 1.0]],      # 2X  X  Z
405        [[1,2,3],[ 1.0, 1.0, 1.0]],      #  X  Y  Z
406        ]
407    indx = GetNXUPQsym(siteSym)
408    return CSxinel[indx[0]]
409   
410def GetCSuinel(siteSym):
411    # returns Uij terms, multipliers, GUI flags & Uiso2Uij multipliers
412    CSuinel = [[],                                             # 0th empty - indices are Fortran style
413        [[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
414        [[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
415        [[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
416        [[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
417        [[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
418        [[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
419        [[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
420        [[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
421        [[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
422        [[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
423        [[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
424        [[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
425        [[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
426        [[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
427        [[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
428        [[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
429        [[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
430        [[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
431        [[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
432        [[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
433        [[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
434        [[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
435        [[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
436        [[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
437        [[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
438        [[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
439        [[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
440        [[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
441        [[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
442        ]
443    indx = GetNXUPQsym(siteSym)
444    return CSuinel[indx[1]]
445       
446def SytSym(XYZ,SGData):
447    '''
448    Generates the number of equivalent positions and a site symmetry code for a specified coordinate and space group
449    input: 
450       XYZ: an array, tuple or list containing 3 elements: x, y & z
451       SGData: from SpcGroup
452    Returns a two element tuple:
453       The 1st element is a code for the site symmetry (see GetOprPtrName)
454       The 2nd element is the site multiplicity
455    '''
456    def PackRot(SGOps):
457        IRT = []
458        for ops in SGOps:
459            M = ops[0]
460            irt = 0
461            for j in range(2,-1,-1):
462                for k in range(2,-1,-1):
463                    irt *= 3
464                    irt += M[k][j]
465            IRT.append(int(irt))
466        return IRT
467       
468    SymName = ''
469    Mult = 1
470    Isym = 0
471    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
472        Isym = 1073741824
473    Jdup = 1
474    Xeqv = GenAtom(XYZ,SGData,True)
475    IRT = PackRot(SGData['SGOps'])
476    L = -1
477    for ic,cen in enumerate(SGData['SGCen']):
478        for invers in range(int(SGData['SGInv']+1)):
479            for io,ops in enumerate(SGData['SGOps']):
480                irtx = (1-2*invers)*IRT[io]
481                L += 1
482                if not Xeqv[L][1]:
483                    Jdup += 1
484                    jx = GetOprPtrName(str(irtx))
485                    if jx[2] < 39:
486                        Isym += 2**(jx[2]-1)
487    if Isym == 1073741824: Isym = 0
488    Mult = len(SGData['SGOps'])*len(SGData['SGCen'])*(int(SGData['SGInv'])+1)/Jdup
489         
490    return GetKNsym(str(Isym)),Mult
491   
492'''A list of space groups as ordered and named in the pre-2002 International
493Tables Volume A, except that spaces are used following the GSAS convention to
494separate the different crystallographic directions.
495Note that the symmetry codes here will recognize many non-standard space group
496symbols with different settings.
497'''
498spglist = {
499    'triclinic' : ('P 1','P -1',), # 1-2
500    'monoclinic': ('P 2','P 21','C 2','P m','P c','C m','C c','P 2/m','P 21/m',
501                   'C 2/m','P 2/c','P 21/c','C 2/c',), #3-15
502    'orthorhombic': ('P 2 2 2','P 2 2 21','P 21 21 2','P 21 21 21','C 2 2 21',
503                     'C 2 2 2','F 2 2 2','I 2 2 2','I 21 21 21',
504                     'P m m 2','P m c 21','P c c 2','P m a 2','P c a 21',
505                     'P n c 2','P m n 21','P b a 2','P n a 21','P n n 2',
506                     'C m m 2','C m c 21','C c c 2','A m m 2','A b m 2',
507                     'A m a 2','A b a 2','F m m 2','F d d 2','I m m 2',
508                     'I b a 2','I m a 2','P m m m','P n n n','P c c m',
509                     'P b a n','P m m a','P n n a','P m n a','P c c a',
510                     'P b a m','P c c n','P b c m','P n n m','P m m n',
511                     'P b c n','P b c a','P n m a','C m c m','C m c a',
512                     'C m m m','C c c m','C m m a','C c c a','F m m m',
513                     'F d d d','I m m m','I b a m','I b c a','I m m a',), #16-74
514    'tetragonal': ('P 4','P 41','P 42','P 43','I 4','I 41','P -4','I -4',
515                   'P 4/m','P 42/m','P 4/n','P 42/n','I 4/m','I 41/a',
516                   'P 4 2 2','P 4 21 2','P 41 2 2','P 41 21 2','P 42 2 2',
517                   'P 42 21 2','P 43 2 2','P 43 21 2','I 4 2 2','I 41 2 2',
518                   'P 4 m m','P 4 b m','P 42 c m','P 42 n m','P 4 c c',
519                   'P 4 n c','P 42 m c','P 42 b c','I 4 m m','I 4 c m',
520                   'I 41 m d','I 41 c d','P -4 2 m','P -4 2 c','P -4 21 m',
521                   'P -4 21 c','P -4 m 2','P -4 c 2','P -4 b 2','P -4 n 2',
522                   'I -4 m 2','I -4 c 2','I -4 2 m','I -4 2 d','P 4/m m m',
523                   'P 4/m c c','P 4/n b m','P 4/n n c','P 4/m b m','P 4/m n c',
524                   'P 4/n m m','P 4/n c c','P 42/m m c','P 42/m c m',
525                   'P 42/n b c','P 42/n n m','P 42/m b c','P 42/m n m',
526                   'P 42/n m c','P 42/n c m','I 4/m m m','I 4/m c m',
527                   'I 41/a m d','I 41/a c d',), # 75-142
528    'trigonal': ('P 3','P 31','P 32','R 3','P -3','R -3','P 3 1 2','P 3 2 1',
529                 'P 31 1 2','P 31 2 1','P 32 1 2','P 32 2 1','R 3 2', 'P 3 m 1',
530                 'P 3 1 m','P 3 c 1','P 3 1 c','R 3 m','R 3 c','P -3 1 m',
531                 'P -3 1 c','P -3 m 1','P -3 c 1','R -3 m','R -3 c',), #143-167
532    'hexagonal': ('P 6','P 61','P 65','P 62','P 64','P 63','P -6','P 6/m',
533                  'P 63/m','P 6 2 2','P 61 2 2','P 65 2 2','P 62 2 2',
534                  'P 64 2 2','P 63 2 2','P 6 m m','P 6 c c','P 63 c m',
535                  'P 63 m c','P -6 m 2','P -6 c 2','P -6 2 m','P -6 2 c',
536                  'P 6/m m m','P 6/m c c','P 63/m c m','P 63/m m c',), #144-194
537    'cubic': ('P 2 3','F 2 3','I 2 3','P 21 3','I 21 3','P m -3','P n -3',
538              'F m -3','F d -3','I m -3','P a -3','I a -3','P 4 3 2','P 42 3 2',
539              'F 4 3 2','F 41 3 2','I 4 3 2','P 43 3 2','P 41 3 2','I 41 3 2',
540              'P -4 3 m','F -4 3 m','I -4 3 m','P -4 3 n','F -4 3 c','I -4 3 d',
541              'P m -3 m','P n -3 n','P m -3 n','P n -3 m','F m -3 m','F m -3 c',
542              'F d -3 m','F d -3 c','I m -3 m','I a -3 d',), #195-230
543}
544'A few non-standard space groups for test use'
545nonstandard_sglist = ('P 21 1 1','P 1 21 1','P 1 1 21','R 3 r','R 3 2 h', 
546                      'R -3 r', 'R 3 2 r','R 3 m h', 'R 3 m r',
547                      'R 3 c r','R -3 c r','R -3 m r',),
548'''A list of orthorhombic space groups that were renamed in the 2002 Volume A,
549along with the pre-2002 name. The e designates a double glide-plane'''
550sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
551                            ('A e a 2', 'A b a 2',),
552                            ('C m c e', 'C m c a',),
553                            ('C m m e', 'C m m a',),
554                            ('C c c e', 'C c c a'),)
555'''Use the space groups types in this order to list the symbols in the
556order they are listed in the International Tables, vol. A'''
557symtypelist = ('triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 
558               'trigonal', 'hexagonal', 'cubic')
559
560# self-test materials follow. Requires files in directory testinp
561def test0():
562    '''test #0: exercise MoveToUnitCell'''
563    msg = "MoveToUnitCell failed"
564    v = [0,1,2,-1,-2]; MoveToUnitCell(v); assert v==[0,0,0,0,0], msg
565    v = np.array([-.1]); MoveToUnitCell(v); assert abs(v-0.9) < 1e-6, msg
566    v = np.array([.1]); MoveToUnitCell(v); assert abs(v-0.1) < 1e-6, msg
567
568def test1():
569    ''' test #1: SpcGroup and SGPrint against previous results'''
570    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
571    if ospath.exists(testdir):
572        if testdir not in sys.path: sys.path.insert(0,testdir)
573    import spctestinp
574    def CompareSpcGroup(spc, referr, refdict, reflist): 
575        'Compare output from GSASIIspc.SpcGroup with results from a previous run'
576        # if an error is reported, the dictionary can be ignored
577        msg = "failed on space group %s" % spc
578        result = SpcGroup(spc)
579        if result[0] == referr and referr > 0: return True
580        keys = result[1].keys()
581        #print result[1]['SpGrp']
582        assert len(keys) == len(refdict.keys()), msg
583        for key in keys:
584        #print key, type(refdict[key])
585            if key == 'SGOps' or  key == 'SGCen':
586                assert len(refdict[key]) == len(result[1][key]), msg
587                for i in range(len(refdict[key])):
588                    assert np.allclose(result[1][key][i][0],refdict[key][i][0]), msg
589                    assert np.allclose(result[1][key][i][1],refdict[key][i][1]), msg
590            else:
591                assert result[1][key] == refdict[key], msg
592        assert reflist == SGPrint(result[1]), 'SGPrint ' +msg
593    for spc in spctestinp.SGdat:
594        CompareSpcGroup(spc, 0, spctestinp.SGdat[spc], spctestinp.SGlist[spc] )
595
596def test2():
597    ''' test #2: SpcGroup against cctbx (sgtbx) computations'''
598    testdir = ospath.join(ospath.split(ospath.abspath( __file__ ))[0],'testinp')
599    if ospath.exists(testdir):
600        if testdir not in sys.path: sys.path.insert(0,testdir)
601    import sgtbxtestinp
602    def CompareWcctbx(spcname, cctbx_in, debug=0):
603        'Compare output from GSASIIspc.SpcGroup with results from cctbx.sgtbx'
604        cctbx = cctbx_in[:] # make copy so we don't delete from the original
605        spc = (SpcGroup(spcname))[1]
606        if debug: print spc['SpGrp']
607        if debug: print spc['SGCen']
608        latticetype = spcname.strip().upper()[0]
609        # lattice type of R implies Hexagonal centering", fix the rhombohedral settings
610        if latticetype == "R" and len(spc['SGCen']) == 1: latticetype = 'P'
611        assert latticetype == spc['SGLatt'], "Failed: %s does not match Lattice: %s" % (spcname, spc['SGLatt'])
612        onebar = [1]
613        if spc['SGInv']: onebar.append(-1)
614        for (op,off) in spc['SGOps']:
615            for inv in onebar:
616                for cen in spc['SGCen']:
617                    noff = off + cen
618                    MoveToUnitCell(noff)
619                    mult = tuple((op*inv).ravel().tolist())
620                    if debug: print "\n%s: %s + %s" % (spcname,mult,noff)
621                    for refop in cctbx:
622                        if debug: print refop
623                        # check the transform
624                        if refop[:9] != mult: continue
625                        if debug: print "mult match"
626                        # check the translation
627                        reftrans = list(refop[-3:])
628                        MoveToUnitCell(reftrans)
629                        if all(abs(noff - reftrans) < 1.e-5):
630                            cctbx.remove(refop)
631                            break
632                    else:
633                        assert False, "failed on %s:\n\t %s + %s" % (spcname,mult,noff)
634    for key in sgtbxtestinp.sgtbx:
635        CompareWcctbx(key, sgtbxtestinp.sgtbx[key])
636
637def test3(): 
638    ''' test #3: exercise SytSym (includes GetOprPtrName, GenAtom, GetKNsym)
639     for selected space groups against info in IT Volume A '''
640    def ExerciseSiteSym (spc, crdlist):
641        'compare site symmetries and multiplicities for a specified space group'
642        msg = "failed on site sym test for %s" % spc
643        (E,S) = SpcGroup(spc)
644        assert not E, msg
645        for t in crdlist:
646            symb, m = SytSym(t[0],S)
647            if symb.strip() != t[2].strip() or m != t[1]:
648                print spc,t[0],m,symb
649            assert m == t[1]
650            #assert symb.strip() == t[2].strip()
651
652    ExerciseSiteSym('p 1',[
653            ((0.13,0.22,0.31),1,'1'),
654            ((0,0,0),1,'1'),
655            ])
656    ExerciseSiteSym('p -1',[
657            ((0.13,0.22,0.31),2,'1'),
658            ((0,0.5,0),1,'-1'),
659            ])
660    ExerciseSiteSym('C 2/c',[
661            ((0.13,0.22,0.31),8,'1'),
662            ((0.0,.31,0.25),4,'2(010)'),
663            ((0.25,.25,0.5),4,'-1'),
664            ((0,0.5,0),4,'-1'),
665            ])
666    ExerciseSiteSym('p 2 2 2',[
667            ((0.13,0.22,0.31),4,'1'),
668            ((0,0.5,.31),2,'2(001)'),
669            ((0.5,.31,0.5),2,'2(010)'),
670            ((.11,0,0),2,'2(100)'),
671            ((0,0.5,0),1,'222'),
672            ])
673    ExerciseSiteSym('p 4/n',[
674            ((0.13,0.22,0.31),8,'1'),
675            ((0.25,0.75,.31),4,'2(001)'),
676            ((0.5,0.5,0.5),4,'-1'),
677            ((0,0.5,0),4,'-1'),
678            ((0.25,0.25,.31),2,'4(001)'),
679            ((0.25,.75,0.5),2,'-4(001)'),
680            ((0.25,.75,0.0),2,'-4(001)'),
681            ])
682    ExerciseSiteSym('p 31 2 1',[
683            ((0.13,0.22,0.31),6,'1'),
684            ((0.13,0.0,0.833333333),3,'2(100)'),
685            ((0.13,0.13,0.),3,'2(110)'),
686            ])
687    ExerciseSiteSym('R 3 c',[
688            ((0.13,0.22,0.31),18,'1'),
689            ((0.0,0.0,0.31),6,'3'),
690            ])
691    ExerciseSiteSym('R 3 c R',[
692            ((0.13,0.22,0.31),6,'1'),
693            ((0.31,0.31,0.31),2,'3(111)'),
694            ])
695    ExerciseSiteSym('P 63 m c',[
696            ((0.13,0.22,0.31),12,'1'),
697            ((0.11,0.22,0.31),6,'m(100)'),
698            ((0.333333,0.6666667,0.31),2,'3m(100)'),
699            ((0,0,0.31),2,'3m(100)'),
700            ])
701    ExerciseSiteSym('I a -3',[
702            ((0.13,0.22,0.31),48,'1'),
703            ((0.11,0,0.25),24,'2(100)'),
704            ((0.11,0.11,0.11),16,'3(111)'),
705            ((0,0,0),8,'-3(111)'),
706            ])
707
708if __name__ == '__main__':
709    test0()
710    test1()
711    test2()
712    test3()
713    print "OK"
Note: See TracBrowser for help on using the repository browser.