source: trunk/GSASIIlattice.py @ 110

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

resolved
File size: 28.2 KB
Line 
1'''Perform lattice-related computations'''
2
3import numpy as np
4import numpy.linalg as nl
5
6# trig functions in degrees
7sind = lambda x: np.sin(x*np.pi/180.)
8asind = lambda x: 180.*np.arcsin(x)/np.pi
9tand = lambda x: np.tan(x*np.pi/180.)
10atand = lambda x: 180.*np.arctan(x)/np.pi
11atan2d = lambda y,x: 180.*np.atan2(y,x)/np.pi
12cosd = lambda x: np.cos(x*np.pi/180.)
13acosd = lambda x: 180.*np.arccos(x)/np.pi
14rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
15
16def sec2HMS(sec):
17    H = int(sec/3600)
18    M = int(sec/60-H*60)
19    S = sec-3600*H-60*M
20    return '%d:%2d:%.2f'%(H,M,S)
21
22def fillgmat(cell):
23    '''Compute lattice metric tensor from unit cell constants
24    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
25    returns 3x3 numpy array
26    '''
27    a,b,c,alp,bet,gam = cell
28    g = np.array([
29        [a*a,  a*b*cosd(gam),  a*c*cosd(bet)],
30        [a*b*cosd(gam),  b*b,  b*c*cosd(alp)],
31        [a*c*cosd(bet) ,b*c*cosd(alp),   c*c]])
32    return g
33           
34def cell2Gmat(cell):
35    '''Compute real and reciprocal lattice metric tensor from unit cell constants
36    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
37    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
38    '''
39    g = fillgmat(cell)
40    G = nl.inv(g)       
41    return G,g
42
43def A2Gmat(A):
44    '''Fill reciprocal metric tensor (G) from A
45    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
46    '''
47    G = np.zeros(shape=(3,3))
48    G = [
49        [A[0],  A[3]/2.,  A[4]/2.], 
50        [A[3]/2.,A[1],    A[5]/2.], 
51        [A[4]/2.,A[5]/2.,    A[2]]]
52    g = nl.inv(G)
53    return G,g
54
55def Gmat2A(G):
56    'Extract A from reciprocal metric tensor (G)'
57    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
58   
59def cell2A(cell):
60    G,g = cell2Gmat(cell)
61    return Gmat2A(G)
62
63def A2cell(A):
64    '''Compute unit cell constants from A tensor
65    returns tuple with a,b,c,alpha, beta, gamma (degrees)
66    '''
67    G,g = A2Gmat(A)
68    return Gmat2cell(g)
69
70def Gmat2cell(g):
71    '''Compute lattice parameters from real metric tensor (g)
72    returns tuple with a,b,c,alpha, beta, gamma (degrees)
73    Alternatively,compute reciprocal lattice parameters from inverse metric tensor (G)
74    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
75    '''
76    a = np.sqrt(max(0,g[0][0]))
77    b = np.sqrt(max(0,g[1][1]))
78    c = np.sqrt(max(0,g[2][2]))
79    alp = acosd(g[2][1]/(b*c))
80    bet = acosd(g[2][0]/(a*c))
81    gam = acosd(g[0][1]/(a*b))
82    return a,b,c,alp,bet,gam
83
84def invcell2Gmat(invcell):
85    '''Compute real and reciprocal lattice metric tensor from reciprocal
86       unit cell constants
87    invcell is tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
88    returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays)
89    '''
90    G = fillgmat(invcell)
91    g = nl.inv(G)
92    return G,g
93       
94def calc_rVsq(A):
95    'Compute the square of the reciprocal lattice volume (V* **2) from A'
96    G,g = A2Gmat(A)
97    rVsq = nl.det(G)
98    if rVsq < 0:
99        return 1
100    return rVsq
101   
102def calc_rV(A):
103    'Compute the reciprocal lattice volume (V*) from A'
104    return np.sqrt(calc_rVsq(A))
105   
106def calc_V(A):
107    'Compute the real lattice volume (V) from A'
108    return 1./calc_rV(A)
109
110def A2invcell(A):
111    '''Compute reciprocal unit cell constants from A
112    returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees)
113    '''
114    G,g = A2Gmat(A)
115    return Gmat2cell(G)
116
117def cell2AB(cell):
118    '''Computes orthogonalization matrix from unit cell constants
119    cell is tuple with a,b,c,alpha, beta, gamma (degrees)
120    returns list of two 3x3 numpy arrays
121       A for crystal to Cartesian transformations A*x = X
122       B (inverse) for Cartesian to crystal transformation B*X = x
123    '''
124    G,g = cell2Gmat(cell) 
125    cellstar = Gmat2cell(G)
126    A = np.zeros(shape=(3,3))
127    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
128    A[0][0] = cell[0]                # a
129    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
130    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
131    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
132    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
133    A[2][2] = 1/cellstar[2]         # 1/c*
134    B = nl.inv(A)
135    return A,B
136
137#reflection generation routines
138#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
139#           cell - a,b,c,alp,bet,gam in A & deg
140   
141def calc_rDsq(H,A):
142    rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
143    return rdsq
144   
145def calc_rDsqZ(H,A,Z,tth,lam):
146    rpd = math.pi/180.
147    rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
148    return rdsq
149       
150def MaxIndex(dmin,A):
151    Hmax = [0,0,0]
152    try:
153        cell = A2cell(A)
154    except:
155        cell = [1,1,1,90,90,90]
156    for i in range(3):
157        Hmax[i] = int(round(cell[i]/dmin))
158    return Hmax
159   
160def sortHKLd(HKLd,ifreverse,ifdup):
161    #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
162    #ifdup = True if duplicate d-spacings allowed
163    T = []
164    for i,H in enumerate(HKLd):
165        if ifdup:
166            T.append((H[3],i))
167        else:
168            T.append(H[3])           
169    D = dict(zip(T,HKLd))
170    T.sort()
171    if ifreverse:
172        T.reverse()
173    X = []
174    okey = ''
175    for key in T: 
176        if key != okey: X.append(D[key])    #remove duplicate d-spacings
177        okey = key
178    return X
179   
180def SwapIndx(Axis,H):
181    if Axis in [1,-1]:
182        return H
183    elif Axis in [2,-3]:
184        return [H[1],H[2],H[0]]
185    else:
186        return [H[2],H[0],H[1]]
187       
188def Rh2Hx(Rh):
189    Hx = [0,0,0]
190    Hx[0] = Rh[0]-Rh[1]
191    Hx[1] = Rh[1]-Rh[2]
192    Hx[2] = np.sum(Rh)
193    return Hx
194   
195def Hx2Rh(Hx):
196        Rh = [0,0,0]
197        itk = -Hx[0]+Hx[1]+Hx[2]
198        if itk%3 != 0:
199            return 0        #error - not rhombohedral reflection
200        else:
201            Rh[1] = itk/3
202            Rh[0] = Rh[1]+Hx[0]
203            Rh[2] = Rh[1]-Hx[1]
204            if Rh[0] < 0:
205                for i in range(3):
206                    Rh[i] = -Rh[i]
207            return Rh
208       
209def CentCheck(Cent,H):
210    h,k,l = H
211    if Cent == 'A' and (k+l)%2:
212        return False
213    elif Cent == 'B' and (h+l)%2:
214        return False
215    elif Cent == 'C' and (h+k)%2:
216        return False
217    elif Cent == 'I' and (h+k+l)%2:
218        return False
219    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
220        return False
221    elif Cent == 'R' and (-h+k+l)%3:
222        return False
223    else:
224        return True
225                                   
226def GetBraviasNum(center,system):
227    '''Determine the Bravais lattice number, as used in GenHBravais
228         center = one of: P, C, I, F, R (see SGLatt from GSASIIspc.SpcGroup)
229         lattice = is cubic, hexagonal, tetragonal, orthorhombic, trigonal (R)
230             monoclinic, triclinic (see SGSys from GSASIIspc.SpcGroup)
231       Returns a number between 0 and 13
232          or throws an exception if the setting is non-standard
233       '''
234    if center.upper() == 'F' and system.lower() == 'cubic':
235        return 0
236    elif center.upper() == 'I' and system.lower() == 'cubic':
237        return 1
238    elif center.upper() == 'P' and system.lower() == 'cubic':
239        return 2
240    elif center.upper() == 'R' and system.lower() == 'trigonal':
241        return 3
242    elif center.upper() == 'P' and system.lower() == 'hexagonal':
243        return 4
244    elif center.upper() == 'I' and system.lower() == 'tetragonal':
245        return 5
246    elif center.upper() == 'P' and system.lower() == 'tetragonal':
247        return 6
248    elif center.upper() == 'F' and system.lower() == 'orthorhombic':
249        return 7
250    elif center.upper() == 'I' and system.lower() == 'orthorhombic':
251        return 8
252    elif center.upper() == 'C' and system.lower() == 'orthorhombic':
253        return 9
254    elif center.upper() == 'P' and system.lower() == 'orthorhombic':
255        return 10
256    elif center.upper() == 'C' and system.lower() == 'monoclinic':
257        return 11
258    elif center.upper() == 'P' and system.lower() == 'monoclinic':
259        return 12
260    elif center.upper() == 'P' and system.lower() == 'triclinic':
261        return 13
262    raise ValueError,'non-standard Bravais lattice center=%s, cell=%s' % (center,system)
263
264def GenHBravais(dmin,Bravais,A):
265    '''Generate the positionally unique powder diffraction reflections
266    input:
267       dmin is minimum d-space
268       Bravais is 0-13 to indicate lattice type (see GetBraviasNum)
269       A is reciprocal cell tensor (see Gmat2A or cell2A)
270    returns:
271       a list of tuples containing: h,k,l,d-space,-1   
272    '''
273# Bravais in range(14) to indicate Bravais lattice:
274#   0 F cubic
275#   1 I cubic
276#   2 P cubic
277#   3 R hexagonal (trigonal not rhombohedral)
278#   4 P hexagonal
279#   5 I tetragonal
280#   6 P tetragonal
281#   7 F orthorhombic
282#   8 I orthorhombic
283#   9 C orthorhombic
284#  10 P orthorhombic
285#  11 C monoclinic
286#  12 P monoclinic
287#  13 P triclinic
288# A - as defined in calc_rDsq
289# returns HKL = [h,k,l,d,0] sorted so d largest first
290    import math
291    if Bravais in [9,11]:
292        Cent = 'C'
293    elif Bravais in [1,5,8]:
294        Cent = 'I'
295    elif Bravais in [0,7]:
296        Cent = 'F'
297    elif Bravais in [3]:
298        Cent = 'R'
299    else:
300        Cent = 'P'
301    Hmax = MaxIndex(dmin,A)
302    dminsq = 1./(dmin**2)
303    HKL = []
304    if Bravais == 13:                       #triclinic
305        for l in range(-Hmax[2],Hmax[2]+1):
306            for k in range(-Hmax[1],Hmax[1]+1):
307                hmin = 0
308                if (k < 0): hmin = 1
309                if (k ==0 and l < 0): hmin = 1
310                for h in range(hmin,Hmax[0]+1):
311                    H=[h,k,l]
312                    rdsq = calc_rDsq(H,A)
313                    if 0 < rdsq <= dminsq:
314                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
315    elif Bravais in [11,12]:                #monoclinic - b unique
316        Hmax = SwapIndx(2,Hmax)
317        for h in range(Hmax[0]+1):
318            for k in range(-Hmax[1],Hmax[1]+1):
319                lmin = 0
320                if k < 0:lmin = 1
321                for l in range(lmin,Hmax[2]+1):
322                    [h,k,l] = SwapIndx(-2,[h,k,l])
323                    H = []
324                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
325                    if H:
326                        rdsq = calc_rDsq(H,A)
327                        if 0 < rdsq <= dminsq:
328                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
329                    [h,k,l] = SwapIndx(2,[h,k,l])
330    elif Bravais in [7,8,9,10]:            #orthorhombic
331        for h in range(Hmax[0]+1):
332            for k in range(Hmax[1]+1):
333                for l in range(Hmax[2]+1):
334                    H = []
335                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
336                    if H:
337                        rdsq = calc_rDsq(H,A)
338                        if 0 < rdsq <= dminsq:
339                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
340    elif Bravais in [5,6]:                  #tetragonal
341        for l in range(Hmax[2]+1):
342            for k in range(Hmax[1]+1):
343                for h in range(k,Hmax[0]+1):
344                    H = []
345                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
346                    if H:
347                        rdsq = calc_rDsq(H,A)
348                        if 0 < rdsq <= dminsq:
349                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
350    elif Bravais in [3,4]:
351        lmin = 0
352        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
353        for l in range(lmin,Hmax[2]+1):
354            for k in range(Hmax[1]+1):
355                hmin = k
356                if l < 0: hmin += 1
357                for h in range(hmin,Hmax[0]+1):
358                    H = []
359                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
360                    if H:
361                        rdsq = calc_rDsq(H,A)
362                        if 0 < rdsq <= dminsq:
363                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
364
365    else:                                   #cubic
366        for l in range(Hmax[2]+1):
367            for k in range(l,Hmax[1]+1):
368                for h in range(k,Hmax[0]+1):
369                    H = []
370                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
371                    if H:
372                        rdsq = calc_rDsq(H,A)
373                        if 0 < rdsq <= dminsq:
374                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
375    return sortHKLd(HKL,True,False)
376   
377def GenHLaue(dmin,Laue,Cent,Axis,A):
378    '''Generate the crystallographically unique powder diffraction reflections
379    for a lattice and Bravais type
380    '''
381# dmin - minimum d-spacing
382# Laue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm',
383#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
384# Cent - lattice centering = 'P','A','B','C','I','F'
385# Axis - code for unique monoclinic axis = 'a','b','c'
386# A - 6 terms as defined in calc_rDsq
387# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
388# part of reciprocal space ignoring anomalous dispersion
389    import math
390    #finds maximum allowed hkl for given A within dmin
391    if Laue in ['3R','3mR']:        #Rhombohedral axes
392        Hmax = [0,0,0]
393        cell = A2cell(A)
394        aHx = cell[0]*math.sqrt(2.0*(1.0-cosd(cell[3])))
395        cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3])))
396        Hmax[0] = Hmax[1] = int(round(aHx/dmin))
397        Hmax[2] = int(round(cHx/dmin))
398        #print Hmax,aHx,cHx
399    else:                           # all others
400        Hmax = MaxIndex(dmin,A)
401       
402    dminsq = 1./(dmin**2)
403    HKL = []
404    if Laue == '-1':                       #triclinic
405        for l in range(-Hmax[2],Hmax[2]+1):
406            for k in range(-Hmax[1],Hmax[1]+1):
407                hmin = 0
408                if (k < 0) or (k ==0 and l < 0): hmin = 1
409                for h in range(hmin,Hmax[0]+1):
410                    H = []
411                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
412                    rdsq = calc_rDsq(H,A)
413                    if 0 < rdsq <= dminsq:
414                        HKL.append([h,k,l,1/math.sqrt(rdsq)])
415    elif Laue == '2/m':                #monoclinic
416        axisnum = 1 + ['a','b','c'].index(Axis)
417        Hmax = SwapIndx(axisnum,Hmax)
418        for h in range(Hmax[0]+1):
419            for k in range(-Hmax[1],Hmax[1]+1):
420                lmin = 0
421                if k < 0:lmin = 1
422                for l in range(lmin,Hmax[2]+1):
423                    [h,k,l] = SwapIndx(-axisnum,[h,k,l])
424                    H = []
425                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
426                    if H:
427                        rdsq = calc_rDsq(H,A)
428                        if 0 < rdsq <= dminsq:
429                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
430                    [h,k,l] = SwapIndx(axisnum,[h,k,l])
431    elif Laue in ['mmm','4/m','6/m']:            #orthorhombic
432        for l in range(Hmax[2]+1):
433            for h in range(Hmax[0]+1):
434                kmin = 1
435                if Laue == 'mmm' or h ==0: kmin = 0
436                for k in range(kmin,Hmax[1]+1):
437                    H = []
438                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
439                    if H:
440                        rdsq = calc_rDsq(H,A)
441                        if 0 < rdsq <= dminsq:
442                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
443    elif Laue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
444        for l in range(Hmax[2]+1):
445            for h in range(Hmax[0]+1):
446                for k in range(h+1):
447                    H = []
448                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
449                    if H:
450                        rdsq = calc_rDsq(H,A)
451                        if 0 < rdsq <= dminsq:
452                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
453    elif Laue in ['3m1','31m','3','3R','3mR']:                  #trigonals
454        for l in range(-Hmax[2],Hmax[2]+1):
455            hmin = 0
456            if l < 0: hmin = 1
457            for h in range(hmin,Hmax[0]+1):
458                if Laue in ['3R','3']:
459                    kmax = h
460                    kmin = -int((h-1.)/2.)
461                else:
462                    kmin = 0
463                    kmax = h
464                    if Laue in ['3m1','3mR'] and l < 0: kmax = h-1
465                    if Laue == '31m' and l < 0: kmin = 1
466                for k in range(kmin,kmax+1):
467                    H = []
468                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
469                    if Laue in ['3R','3mR']:
470                        H = Hx2Rh(H)
471                    if H:
472                        rdsq = calc_rDsq(H,A)
473                        if 0 < rdsq <= dminsq:
474                            HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)])
475    else:                                   #cubic
476        for h in range(Hmax[0]+1):
477            for k in range(h+1):
478                lmin = 0
479                lmax = k
480                if Laue =='m3':
481                    lmax = h-1
482                    if h == k: lmax += 1
483                for l in range(lmin,lmax+1):
484                    H = []
485                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
486                    if H:
487                        rdsq = calc_rDsq(H,A)
488                        if 0 < rdsq <= dminsq:
489                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
490    return sortHKLd(HKL,True,True)
491   
492# output from uctbx computed on platform darwin on 2010-05-28
493NeedTestData = True
494def TestData():
495    array = np.array
496    global NeedTestData
497    NeedTestData = False
498    global CellTestData
499    CellTestData = [
500# cell, g, G, cell*, V, V*
501  [(4, 4, 4, 90, 90, 90), 
502   array([[  1.60000000e+01,   9.79717439e-16,   9.79717439e-16],
503       [  9.79717439e-16,   1.60000000e+01,   9.79717439e-16],
504       [  9.79717439e-16,   9.79717439e-16,   1.60000000e+01]]), array([[  6.25000000e-02,   3.82702125e-18,   3.82702125e-18],
505       [  3.82702125e-18,   6.25000000e-02,   3.82702125e-18],
506       [  3.82702125e-18,   3.82702125e-18,   6.25000000e-02]]), (0.25, 0.25, 0.25, 90.0, 90.0, 90.0), 64.0, 0.015625],
507# cell, g, G, cell*, V, V*
508  [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), 
509   array([[ 16.81      , -13.70423184,   4.48533243],
510       [-13.70423184,  27.04      ,  -5.6887143 ],
511       [  4.48533243,  -5.6887143 ,  39.69      ]]), array([[ 0.10206349,  0.05083339, -0.00424823],
512       [ 0.05083339,  0.06344997,  0.00334956],
513       [-0.00424823,  0.00334956,  0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445],
514# cell, g, G, cell*, V, V*
515  [(3.5, 3.5, 6, 90, 90, 120), 
516   array([[  1.22500000e+01,  -6.12500000e+00,   1.28587914e-15],
517       [ -6.12500000e+00,   1.22500000e+01,   1.28587914e-15],
518       [  1.28587914e-15,   1.28587914e-15,   3.60000000e+01]]), array([[  1.08843537e-01,   5.44217687e-02,   3.36690552e-18],
519       [  5.44217687e-02,   1.08843537e-01,   3.36690552e-18],
520       [  3.36690552e-18,   3.36690552e-18,   2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
521  ]
522    global CoordTestData
523    CoordTestData = [
524# cell, ((frac, ortho),...)
525  ((4,4,4,90,90,90,), [
526 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
527 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
528 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
529 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
530 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
531 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
532 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
533]),
534# cell, ((frac, ortho),...)
535  ((4.1,5.2,6.3,100,80,130,), [
536 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
537 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
538 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
539 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
540 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
541 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
542 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
543]),
544# cell, ((frac, ortho),...)
545  ((3.5,3.5,6,90,90,120,), [
546 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
547 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
548 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
549 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
550 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
551 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
552 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
553]),
554]
555
556def test0():
557    if NeedTestData: TestData()
558    msg = 'test cell2Gmat, fillgmat, Gmat2cell'
559    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
560        G, g = cell2Gmat(cell)
561        assert np.allclose(G,tG),msg
562        assert np.allclose(g,tg),msg
563        tcell = Gmat2cell(g)
564        assert np.allclose(cell,tcell),msg
565        tcell = Gmat2cell(G)
566        assert np.allclose(tcell,trcell),msg
567
568def test1():
569    if NeedTestData: TestData()
570    msg = 'test cell2A and A2Gmat'
571    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
572        G, g = A2Gmat(cell2A(cell))
573        assert np.allclose(G,tG),msg
574        assert np.allclose(g,tg),msg
575
576def test2():
577    if NeedTestData: TestData()
578    msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell'
579    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
580        G, g = cell2Gmat(cell)
581        tcell = A2cell(Gmat2A(G))
582        assert np.allclose(cell,tcell),msg
583
584def test3():
585    if NeedTestData: TestData()
586    msg = 'test invcell2Gmat'
587    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
588        G, g = invcell2Gmat(trcell)
589        assert np.allclose(G,tG),msg
590        assert np.allclose(g,tg),msg
591
592def test4():
593    if NeedTestData: TestData()
594    msg = 'test calc_rVsq, calc_rV, calc_V'
595    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
596        assert np.allclose(calc_rV(cell2A(cell)),trV), msg
597        assert np.allclose(calc_V(cell2A(cell)),tV), msg
598
599def test5():
600    if NeedTestData: TestData()
601    msg = 'test A2invcell'
602    for (cell, tg, tG, trcell, tV, trV) in CellTestData:
603        rcell = A2invcell(cell2A(cell))
604        assert np.allclose(rcell,trcell),msg
605
606def test6():
607    if NeedTestData: TestData()
608    msg = 'test cell2AB'
609    for (cell,coordlist) in CoordTestData:
610        A,B = cell2AB(cell)
611        for (frac,ortho) in coordlist:
612            to = np.inner(A,frac)
613            tf = np.inner(B,to)
614            assert np.allclose(ortho,to), msg
615            assert np.allclose(frac,tf), msg
616            to = np.sum(A*frac,axis=1)
617            tf = np.sum(B*to,axis=1)
618            assert np.allclose(ortho,to), msg
619            assert np.allclose(frac,tf), msg
620
621# test GetBraviasNum(...) and GenHBravais(...)
622def test7():
623    import os.path
624    import sys
625    import GSASIIspc as spc
626    testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp')
627    if os.path.exists(testdir):
628        if testdir not in sys.path: sys.path.insert(0,testdir)
629    import sgtbxlattinp
630    derror = 1e-4
631    def indexmatch(hklin, hkllist, system):
632        for hklref in hkllist:
633            hklref = list(hklref)
634            # these permutations are far from complete, but are sufficient to
635            # allow the test to complete
636            if system == 'cubic':
637                permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
638            elif system == 'monoclinic':
639                permlist = [(1,2,3),(-1,2,-3)]
640            else:
641                permlist = [(1,2,3)]
642
643            for perm in permlist:
644                hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
645                if hkl == hklref: return True
646                if [-i for i in hkl] == hklref: return True
647        else:
648            return False
649
650    for key in sgtbxlattinp.sgtbx7:
651        spdict = spc.SpcGroup(key)
652        cell = sgtbxlattinp.sgtbx7[key][0]
653        system = spdict[1]['SGSys']
654        center = spdict[1]['SGLatt']
655
656        bravcode = GetBraviasNum(center, system)
657
658        g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell))
659
660        assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key
661        for h,k,l,d,num in g2list:
662            for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: 
663                if abs(d-dref) < derror:
664                    if indexmatch((h,k,l,), hkllist, system):
665                        break
666            else:
667                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
668
669def test8():
670    import GSASIIspc as spc
671    import sgtbxlattinp
672    derror = 1e-4
673    dmin = sgtbxlattinp.dmin
674
675    def indexmatch(hklin, hklref, system, axis):
676        # these permutations are far from complete, but are sufficient to
677        # allow the test to complete
678        if system == 'cubic':
679            permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),]
680        elif system == 'monoclinic' and axis=='b':
681            permlist = [(1,2,3),(-1,2,-3)]
682        elif system == 'monoclinic' and axis=='a':
683            permlist = [(1,2,3),(1,-2,-3)]
684        elif system == 'monoclinic' and axis=='c':
685            permlist = [(1,2,3),(-1,-2,3)]
686        elif system == 'trigonal':
687            permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)]
688        elif system == 'rhombohedral':
689            permlist = [(1,2,3),(2,3,1),(3,1,2),(-1,-2,-3),(-2,-3,-1),(-3,-1,-2)]
690        else:
691            permlist = [(1,2,3)]
692
693        hklref = list(hklref)
694        for perm in permlist:
695            hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm]
696            if hkl == hklref: return True
697            if [-i for i in hkl] == hklref: return True
698        return False
699
700    for key in sgtbxlattinp.sgtbx8:
701        spdict = spc.SpcGroup(key)[1]
702        cell = sgtbxlattinp.sgtbx8[key][0]
703        center = spdict['SGLatt']
704        Laue = spdict['SGLaue']
705        Axis = spdict['SGUniq']
706        system = spdict['SGSys']
707
708        g2list = GenHLaue(dmin,Laue,center,Axis,cell2A(cell))
709        #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]):
710        #    print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1])
711        #    print 'GSAS-II:'
712        #    for h,k,l,d in g2list: print '  ',(h,k,l),d
713        #    print 'SGTBX:'
714        #    for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
715        assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), (
716            'Reflection lists differ for %s' % key
717            )
718        #match = True
719        for h,k,l,d in g2list:
720            for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: 
721                if abs(d-dref) < derror:
722                    if indexmatch((h,k,l,), hkllist, system, Axis): break
723            else:
724                assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key)
725                #match = False
726        #if not match:
727            #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print '  ',hkllist,dref
728            #print center, Laue, Axis, system
729
730if __name__ == '__main__':
731    test0()
732    test1()
733    test2()
734    test3()
735    test4()
736    test5()
737    test6()
738    test7()
739    test8()
740    print "OK"
Note: See TracBrowser for help on using the repository browser.