source: trunk/GSASIIlattice.py @ 296

Last change on this file since 296 was 296, checked in by vondreele, 14 years ago

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