source: trunk/GSASIIlattice.py @ 189

Last change on this file since 189 was 189, checked in by vondreele, 13 years ago

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