1 | '''Perform lattice-related computations''' |
---|
2 | |
---|
3 | ########### SVN repository information ################### |
---|
4 | # $Date: 2011-06-20 19:37:07 +0000 (Mon, 20 Jun 2011) $ |
---|
5 | # $Author: vondreele $ |
---|
6 | # $Revision: 303 $ |
---|
7 | # $URL: trunk/GSASIIlattice.py $ |
---|
8 | # $Id: GSASIIlattice.py 303 2011-06-20 19:37:07Z vondreele $ |
---|
9 | ########### SVN repository information ################### |
---|
10 | import math |
---|
11 | import numpy as np |
---|
12 | import numpy.linalg as nl |
---|
13 | |
---|
14 | # trig functions in degrees |
---|
15 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
16 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
17 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
18 | atand = lambda x: 180.*np.arctan(x)/np.pi |
---|
19 | atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
20 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
21 | acosd = lambda x: 180.*np.arccos(x)/np.pi |
---|
22 | rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p) |
---|
23 | |
---|
24 | def 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 | |
---|
30 | def 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 | |
---|
41 | def 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 | |
---|
48 | def 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 | |
---|
60 | def 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 | |
---|
69 | def 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 | |
---|
81 | def 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 | |
---|
85 | def cell2A(cell): |
---|
86 | G,g = cell2Gmat(cell) |
---|
87 | return Gmat2A(G) |
---|
88 | |
---|
89 | def 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 | |
---|
96 | def 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 | |
---|
110 | def 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 | |
---|
120 | def 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 | |
---|
128 | def calc_rV(A): |
---|
129 | 'Compute the reciprocal lattice volume (V*) from A' |
---|
130 | return np.sqrt(calc_rVsq(A)) |
---|
131 | |
---|
132 | def calc_V(A): |
---|
133 | 'Compute the real lattice volume (V) from A' |
---|
134 | return 1./calc_rV(A) |
---|
135 | |
---|
136 | def 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 | |
---|
143 | def 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 | |
---|
163 | def 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 | |
---|
174 | def 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 | |
---|
185 | def 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 | |
---|
203 | def 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 | # |
---|
214 | def _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 |
---|
223 | def 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) |
---|
228 | def 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) |
---|
233 | def 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) |
---|
238 | def 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 | |
---|
246 | def 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 | |
---|
250 | def calc_rDsq2(H,G): |
---|
251 | return np.inner(H,np.inner(G,H)) |
---|
252 | |
---|
253 | def 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 | |
---|
258 | def 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 | |
---|
268 | def 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 | |
---|
288 | def 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 | |
---|
296 | def 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 | |
---|
303 | def 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 | |
---|
317 | def 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 | |
---|
334 | def 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 | |
---|
372 | def 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 | |
---|
485 | def 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 | #Spherical harmonics routines |
---|
601 | def OdfChk(SGLaue,L,M): |
---|
602 | if not L%2 and abs(M) <= L: |
---|
603 | if SGLaue == '0': #cylindrical symmetry |
---|
604 | if M == 0: return True |
---|
605 | elif SGLaue == '-1': |
---|
606 | return True |
---|
607 | elif SGLaue == '2/m': |
---|
608 | if not abs(M)%2: return True |
---|
609 | elif SGLaue == 'mmm': |
---|
610 | if not abs(M)%2 and M >= 0: return True |
---|
611 | elif SGLaue == '4/m': |
---|
612 | if not abs(M)%4: return True |
---|
613 | elif SGLaue == '4/mmm': |
---|
614 | if not abs(M)%4 and M >= 0: return True |
---|
615 | elif SGLaue in ['3R','3']: |
---|
616 | if not abs(M)%3: return True |
---|
617 | elif SGLaue in ['3mR','3m1','31m']: |
---|
618 | if not abs(M)%3 and M >= 0: return True |
---|
619 | elif SGLaue == '6/m': |
---|
620 | if not abs(M)%6: return True |
---|
621 | elif SGLaue == '6/mmm': |
---|
622 | if not abs(M)%6 and M >= 0: return True |
---|
623 | elif SGLaue == 'm3': |
---|
624 | if M > 0: |
---|
625 | if L%12 == 2: |
---|
626 | if M <= L/12: return True |
---|
627 | else: |
---|
628 | if M <= L/12+1: return True |
---|
629 | elif SGLaue == 'm3m': |
---|
630 | if M > 0: |
---|
631 | if L%12 == 2: |
---|
632 | if M <= L/12: return True |
---|
633 | else: |
---|
634 | if M <= L/12+1: return True |
---|
635 | return False |
---|
636 | |
---|
637 | def GenSHCoeff(SGLaue,SamSym,L): |
---|
638 | coeffNames = [] |
---|
639 | for iord in [2*i+2 for i in range(L/2)]: |
---|
640 | for m in [i-iord for i in range(2*iord+1)]: |
---|
641 | if OdfChk(SamSym,iord,m): |
---|
642 | for n in [i-iord for i in range(2*iord+1)]: |
---|
643 | if OdfChk(SGLaue,iord,n): |
---|
644 | coeffNames.append('C(%d,%d,%d)'%(iord,m,n)) |
---|
645 | return coeffNames |
---|
646 | |
---|
647 | def CrsAng(H,cell,SGData): |
---|
648 | a,b,c,al,be,ga = cell |
---|
649 | SQ3 = 1.732050807569 |
---|
650 | H1 = np.array([1,0,0]) |
---|
651 | H2 = np.array([0,1,0]) |
---|
652 | H3 = np.array([0,0,1]) |
---|
653 | H4 = np.array([1,1,1]) |
---|
654 | G,g = cell2Gmat(cell) |
---|
655 | Laue = SGData['SGLaue'] |
---|
656 | Naxis = SGData['SGUniq'] |
---|
657 | DH = np.inner(H,np.inner(G,H)) |
---|
658 | if Laue == '2/m': |
---|
659 | if Naxis == 'a': |
---|
660 | DR = np.inner(H1,np.inner(G,H1)) |
---|
661 | DHR = np.inner(H,np.inner(G,H1)) |
---|
662 | elif Naxis == 'b': |
---|
663 | DR = np.inner(H2,np.inner(G,H2)) |
---|
664 | DHR = np.inner(H,np.inner(G,H2)) |
---|
665 | else: |
---|
666 | DR = np.inner(H3,np.inner(G,H3)) |
---|
667 | DHR = np.inner(H,np.inner(G,H3)) |
---|
668 | elif Laue in ['R3','R3m']: |
---|
669 | DR = np.inner(H4,np.inner(G,H4)) |
---|
670 | DRH = np.inner(H,np.inner(G,H4)) |
---|
671 | |
---|
672 | else: |
---|
673 | DR = np.inner(H3,np.inner(G,H3)) |
---|
674 | DHR = np.inner(H,np.inner(G,H3)) |
---|
675 | DHR /= np.sqrt(DR*DH) |
---|
676 | phi = acosd(DHR) |
---|
677 | if Laue == '-1': |
---|
678 | BA = H[1]*a/(b-H[0]*cosd(ga)) |
---|
679 | BB = H[0]*sind(ga)**2 |
---|
680 | elif Laue == '2/m': |
---|
681 | if Naxis == 'a': |
---|
682 | BA = H[2]*b/(c-H[1]*cods(al)) |
---|
683 | BB = H[1]*sind(al)**2 |
---|
684 | elif Naxis == 'b': |
---|
685 | BA = H[0]*c/(a-H[2]*cosd(be)) |
---|
686 | BB = H[2]*sind(be)**2 |
---|
687 | else: |
---|
688 | BA = H[1]*a/(b-H[0]*cosd(ga)) |
---|
689 | BB = H[0]*sind(ga)**2 |
---|
690 | elif Laue in ['mmm','4/m','4/mmm']: |
---|
691 | BA = H[1]*a |
---|
692 | BB = H[0]*b |
---|
693 | |
---|
694 | elif Laue in ['3R','3mR']: |
---|
695 | BA = H[0]+H[1]-2.0*H[2] |
---|
696 | BB = SQ3*(H[0]-H[1]) |
---|
697 | elif Laue in ['m3','m3m']: |
---|
698 | BA = H[1] |
---|
699 | BB = H[0] |
---|
700 | else: |
---|
701 | BA = H[0]+2.0*H[1] |
---|
702 | BB = SQ3*H[0] |
---|
703 | beta = atan2d(BA,BB) |
---|
704 | return phi,beta |
---|
705 | |
---|
706 | def SamAng(Tth,Gangls,Sangl,IFCoup): |
---|
707 | if IFCoup: |
---|
708 | GSomeg = sind(Gangls[2]+Tth) |
---|
709 | GComeg = cosd(Gangls[2]+Tth) |
---|
710 | else: |
---|
711 | GSomeg = sind(Gangls[2]) |
---|
712 | GComeg = cosd(Gangls[2]) |
---|
713 | GSTth = sind(Tth) |
---|
714 | GCTth = cosd(Tth) |
---|
715 | GSazm = sind(Gangls[3]) |
---|
716 | GCazm = cosd(Gangls[3]) |
---|
717 | GSchi = sind(Gangls[1]) |
---|
718 | GCchi = cosd(Gangls[1]) |
---|
719 | GSphi = sind(Gangls[0]+Sangl[2]) |
---|
720 | GCphi = cosd(Gangls[0]+Sangl[2]) |
---|
721 | SSomeg = sind(Sangl[0]) |
---|
722 | SComeg = cosd(Sangl[0]) |
---|
723 | SSchi = sind(Sangl[1]) |
---|
724 | SCchi = cosd(Sangl[1]) |
---|
725 | AT = -GSTth*GComeg+GCTth*GCazm*GSomeg |
---|
726 | BT = GSTth*GSomeg+GCTth*GCazm*GComeg |
---|
727 | CT = -GCTth*GSazm*GSchi |
---|
728 | DT = -GCTth*GSazm*GCchi |
---|
729 | |
---|
730 | BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi |
---|
731 | BC2 = DT-BT*GSchi |
---|
732 | BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi |
---|
733 | |
---|
734 | BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg |
---|
735 | psi = acosd(BC) |
---|
736 | |
---|
737 | BA = -BC1*SSchi+BC2*SCchi |
---|
738 | BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg |
---|
739 | gam = atand2(BB,BA) |
---|
740 | |
---|
741 | return psi,gam |
---|
742 | |
---|
743 | def Flnh(Start,SHCoef,phi,beta,SGData): |
---|
744 | import pytexture as ptx |
---|
745 | BOH = { |
---|
746 | 'L=2':[[],[],[]], |
---|
747 | 'L=4':[[0.30469720,0.36418281],[],[]], |
---|
748 | 'L=6':[[-0.14104740,0.52775103],[],[]], |
---|
749 | 'L=8':[[0.28646862,0.21545346,0.32826995],[],[]], |
---|
750 | 'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]], |
---|
751 | 'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402], |
---|
752 | [0.09298802,-0.23773812,0.49446631],[]], |
---|
753 | 'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]], |
---|
754 | 'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495], |
---|
755 | [0.12039646,-0.25330128,0.23950998,0.40962508],[]], |
---|
756 | 'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037], |
---|
757 | [-0.06901768,0.16006562,-0.24743528,0.47110273],[]], |
---|
758 | 'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563], |
---|
759 | [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433],[]], |
---|
760 | 'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583], |
---|
761 | [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231],[]], |
---|
762 | 'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505], |
---|
763 | [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497], |
---|
764 | [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956]], |
---|
765 | 'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594], |
---|
766 | [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355],[]], |
---|
767 | 'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501], |
---|
768 | [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598], |
---|
769 | [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602]], |
---|
770 | 'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594], |
---|
771 | [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385], |
---|
772 | [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361]], |
---|
773 | 'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501], |
---|
774 | [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978], |
---|
775 | [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402]], |
---|
776 | 'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508], |
---|
777 | [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152], |
---|
778 | [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318]] |
---|
779 | } |
---|
780 | |
---|
781 | FORPI = 12.5663706143592 |
---|
782 | RSQPI = 0.5641895835478 |
---|
783 | SQ2 = 1.414213562373 |
---|
784 | |
---|
785 | if Start: |
---|
786 | ptx.pyqlmninit() |
---|
787 | Start = False |
---|
788 | Fln = np.zeros(len(SHCoef)) |
---|
789 | for i,term in enumerate(SHCoef): |
---|
790 | if abs(SHCoef[term]) > 1e-6: |
---|
791 | l,m,n = eval(term.strip('C')) |
---|
792 | lNorm = 4.*np.pi*(2.*l+1.) |
---|
793 | if SGData['SGLaue'] in ['m3','m3m']: |
---|
794 | L = [i*4 for i in range(l/4)] |
---|
795 | Kcl = 0.0 |
---|
796 | for i in L: |
---|
797 | im = i/4 |
---|
798 | pcrs = ptx.pyplmpsi(l,i,phi) |
---|
799 | Kcl += BOH['L='+str(l)][n-1][l/2-1]*pcrs*cosd(i*beta) |
---|
800 | else: #all but cubic |
---|
801 | pcrs = ptx.pyplmpsi(l,n,phi)*RSQPI |
---|
802 | if not n: |
---|
803 | pcrs /= SQ2 |
---|
804 | if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm']: |
---|
805 | if n%6 == 3: |
---|
806 | Kcl = pcrs*sind(n*beta) |
---|
807 | else: |
---|
808 | Kcl = pcrs*cosd(n*beta) |
---|
809 | elif SGData['SGLaue'] in ['3mR','3m1','31m']: |
---|
810 | Kcl = pcrs*cosd(n*beta) |
---|
811 | else: |
---|
812 | Kcl = pcrs*(cosd(n*beta)+sind(n*beta)) |
---|
813 | |
---|
814 | Fln[i] = SHCoef[term]*Kcl*lNorm |
---|
815 | ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln)))) |
---|
816 | return ODFln |
---|
817 | |
---|
818 | def polfcal(ODFln,SamSym,psi,gam): |
---|
819 | RSQPI = 0.5641895835478 |
---|
820 | SQ2 = 1.414213562373 |
---|
821 | PolVal = 1.0 |
---|
822 | for term in ODFln: |
---|
823 | if ODFln[term][1] > 1.e-3: |
---|
824 | l,m,n = eval(term.strip('C')) |
---|
825 | psrs = ptx.pyplmpsi(l,m,psi) |
---|
826 | if SamSym in ['-1','2/m']: |
---|
827 | if m: |
---|
828 | Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam)) |
---|
829 | else: |
---|
830 | ksl = RSQPI*psrs/SQ2 |
---|
831 | else: |
---|
832 | if m: |
---|
833 | Ksl = RSQPI*psrs*cosd(m*gam) |
---|
834 | else: |
---|
835 | Ksl = RSQPI*psrs/SQ2 |
---|
836 | PolVal += ODFln[term][1]*Ksl |
---|
837 | return PolVal |
---|
838 | |
---|
839 | # output from uctbx computed on platform darwin on 2010-05-28 |
---|
840 | NeedTestData = True |
---|
841 | def TestData(): |
---|
842 | array = np.array |
---|
843 | global NeedTestData |
---|
844 | NeedTestData = False |
---|
845 | global CellTestData |
---|
846 | CellTestData = [ |
---|
847 | # cell, g, G, cell*, V, V* |
---|
848 | [(4, 4, 4, 90, 90, 90), |
---|
849 | array([[ 1.60000000e+01, 9.79717439e-16, 9.79717439e-16], |
---|
850 | [ 9.79717439e-16, 1.60000000e+01, 9.79717439e-16], |
---|
851 | [ 9.79717439e-16, 9.79717439e-16, 1.60000000e+01]]), array([[ 6.25000000e-02, 3.82702125e-18, 3.82702125e-18], |
---|
852 | [ 3.82702125e-18, 6.25000000e-02, 3.82702125e-18], |
---|
853 | [ 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], |
---|
854 | # cell, g, G, cell*, V, V* |
---|
855 | [(4.0999999999999996, 5.2000000000000002, 6.2999999999999998, 100, 80, 130), |
---|
856 | array([[ 16.81 , -13.70423184, 4.48533243], |
---|
857 | [-13.70423184, 27.04 , -5.6887143 ], |
---|
858 | [ 4.48533243, -5.6887143 , 39.69 ]]), array([[ 0.10206349, 0.05083339, -0.00424823], |
---|
859 | [ 0.05083339, 0.06344997, 0.00334956], |
---|
860 | [-0.00424823, 0.00334956, 0.02615544]]), (0.31947376387537696, 0.25189277536327803, 0.16172643497798223, 85.283666420376008, 94.716333579624006, 50.825714168082683), 100.98576357983838, 0.0099023858863968445], |
---|
861 | # cell, g, G, cell*, V, V* |
---|
862 | [(3.5, 3.5, 6, 90, 90, 120), |
---|
863 | array([[ 1.22500000e+01, -6.12500000e+00, 1.28587914e-15], |
---|
864 | [ -6.12500000e+00, 1.22500000e+01, 1.28587914e-15], |
---|
865 | [ 1.28587914e-15, 1.28587914e-15, 3.60000000e+01]]), array([[ 1.08843537e-01, 5.44217687e-02, 3.36690552e-18], |
---|
866 | [ 5.44217687e-02, 1.08843537e-01, 3.36690552e-18], |
---|
867 | [ 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], |
---|
868 | ] |
---|
869 | global CoordTestData |
---|
870 | CoordTestData = [ |
---|
871 | # cell, ((frac, ortho),...) |
---|
872 | ((4,4,4,90,90,90,), [ |
---|
873 | ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)), |
---|
874 | ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)), |
---|
875 | ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)), |
---|
876 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)), |
---|
877 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)), |
---|
878 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)), |
---|
879 | ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)), |
---|
880 | ]), |
---|
881 | # cell, ((frac, ortho),...) |
---|
882 | ((4.1,5.2,6.3,100,80,130,), [ |
---|
883 | ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)), |
---|
884 | ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)), |
---|
885 | ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)), |
---|
886 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)), |
---|
887 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)), |
---|
888 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)), |
---|
889 | ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)), |
---|
890 | ]), |
---|
891 | # cell, ((frac, ortho),...) |
---|
892 | ((3.5,3.5,6,90,90,120,), [ |
---|
893 | ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)), |
---|
894 | ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)), |
---|
895 | ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)), |
---|
896 | ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)), |
---|
897 | ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)), |
---|
898 | ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)), |
---|
899 | ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)), |
---|
900 | ]), |
---|
901 | ] |
---|
902 | |
---|
903 | def test0(): |
---|
904 | if NeedTestData: TestData() |
---|
905 | msg = 'test cell2Gmat, fillgmat, Gmat2cell' |
---|
906 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
907 | G, g = cell2Gmat(cell) |
---|
908 | assert np.allclose(G,tG),msg |
---|
909 | assert np.allclose(g,tg),msg |
---|
910 | tcell = Gmat2cell(g) |
---|
911 | assert np.allclose(cell,tcell),msg |
---|
912 | tcell = Gmat2cell(G) |
---|
913 | assert np.allclose(tcell,trcell),msg |
---|
914 | |
---|
915 | def test1(): |
---|
916 | if NeedTestData: TestData() |
---|
917 | msg = 'test cell2A and A2Gmat' |
---|
918 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
919 | G, g = A2Gmat(cell2A(cell)) |
---|
920 | assert np.allclose(G,tG),msg |
---|
921 | assert np.allclose(g,tg),msg |
---|
922 | |
---|
923 | def test2(): |
---|
924 | if NeedTestData: TestData() |
---|
925 | msg = 'test Gmat2A, A2cell, A2Gmat, Gmat2cell' |
---|
926 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
927 | G, g = cell2Gmat(cell) |
---|
928 | tcell = A2cell(Gmat2A(G)) |
---|
929 | assert np.allclose(cell,tcell),msg |
---|
930 | |
---|
931 | def test3(): |
---|
932 | if NeedTestData: TestData() |
---|
933 | msg = 'test invcell2Gmat' |
---|
934 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
935 | G, g = invcell2Gmat(trcell) |
---|
936 | assert np.allclose(G,tG),msg |
---|
937 | assert np.allclose(g,tg),msg |
---|
938 | |
---|
939 | def test4(): |
---|
940 | if NeedTestData: TestData() |
---|
941 | msg = 'test calc_rVsq, calc_rV, calc_V' |
---|
942 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
943 | assert np.allclose(calc_rV(cell2A(cell)),trV), msg |
---|
944 | assert np.allclose(calc_V(cell2A(cell)),tV), msg |
---|
945 | |
---|
946 | def test5(): |
---|
947 | if NeedTestData: TestData() |
---|
948 | msg = 'test A2invcell' |
---|
949 | for (cell, tg, tG, trcell, tV, trV) in CellTestData: |
---|
950 | rcell = A2invcell(cell2A(cell)) |
---|
951 | assert np.allclose(rcell,trcell),msg |
---|
952 | |
---|
953 | def test6(): |
---|
954 | if NeedTestData: TestData() |
---|
955 | msg = 'test cell2AB' |
---|
956 | for (cell,coordlist) in CoordTestData: |
---|
957 | A,B = cell2AB(cell) |
---|
958 | for (frac,ortho) in coordlist: |
---|
959 | to = np.inner(A,frac) |
---|
960 | tf = np.inner(B,to) |
---|
961 | assert np.allclose(ortho,to), msg |
---|
962 | assert np.allclose(frac,tf), msg |
---|
963 | to = np.sum(A*frac,axis=1) |
---|
964 | tf = np.sum(B*to,axis=1) |
---|
965 | assert np.allclose(ortho,to), msg |
---|
966 | assert np.allclose(frac,tf), msg |
---|
967 | |
---|
968 | # test GetBraviasNum(...) and GenHBravais(...) |
---|
969 | def test7(): |
---|
970 | import os.path |
---|
971 | import sys |
---|
972 | import GSASIIspc as spc |
---|
973 | testdir = os.path.join(os.path.split(os.path.abspath( __file__ ))[0],'testinp') |
---|
974 | if os.path.exists(testdir): |
---|
975 | if testdir not in sys.path: sys.path.insert(0,testdir) |
---|
976 | import sgtbxlattinp |
---|
977 | derror = 1e-4 |
---|
978 | def indexmatch(hklin, hkllist, system): |
---|
979 | for hklref in hkllist: |
---|
980 | hklref = list(hklref) |
---|
981 | # these permutations are far from complete, but are sufficient to |
---|
982 | # allow the test to complete |
---|
983 | if system == 'cubic': |
---|
984 | permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),] |
---|
985 | elif system == 'monoclinic': |
---|
986 | permlist = [(1,2,3),(-1,2,-3)] |
---|
987 | else: |
---|
988 | permlist = [(1,2,3)] |
---|
989 | |
---|
990 | for perm in permlist: |
---|
991 | hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm] |
---|
992 | if hkl == hklref: return True |
---|
993 | if [-i for i in hkl] == hklref: return True |
---|
994 | else: |
---|
995 | return False |
---|
996 | |
---|
997 | for key in sgtbxlattinp.sgtbx7: |
---|
998 | spdict = spc.SpcGroup(key) |
---|
999 | cell = sgtbxlattinp.sgtbx7[key][0] |
---|
1000 | system = spdict[1]['SGSys'] |
---|
1001 | center = spdict[1]['SGLatt'] |
---|
1002 | |
---|
1003 | bravcode = GetBraviasNum(center, system) |
---|
1004 | |
---|
1005 | g2list = GenHBravais(sgtbxlattinp.dmin, bravcode, cell2A(cell)) |
---|
1006 | |
---|
1007 | assert len(sgtbxlattinp.sgtbx7[key][1]) == len(g2list), 'Reflection lists differ for %s' % key |
---|
1008 | for h,k,l,d,num in g2list: |
---|
1009 | for hkllist,dref in sgtbxlattinp.sgtbx7[key][1]: |
---|
1010 | if abs(d-dref) < derror: |
---|
1011 | if indexmatch((h,k,l,), hkllist, system): |
---|
1012 | break |
---|
1013 | else: |
---|
1014 | assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key) |
---|
1015 | |
---|
1016 | def test8(): |
---|
1017 | import GSASIIspc as spc |
---|
1018 | import sgtbxlattinp |
---|
1019 | derror = 1e-4 |
---|
1020 | dmin = sgtbxlattinp.dmin |
---|
1021 | |
---|
1022 | def indexmatch(hklin, hklref, system, axis): |
---|
1023 | # these permutations are far from complete, but are sufficient to |
---|
1024 | # allow the test to complete |
---|
1025 | if system == 'cubic': |
---|
1026 | permlist = [(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1),] |
---|
1027 | elif system == 'monoclinic' and axis=='b': |
---|
1028 | permlist = [(1,2,3),(-1,2,-3)] |
---|
1029 | elif system == 'monoclinic' and axis=='a': |
---|
1030 | permlist = [(1,2,3),(1,-2,-3)] |
---|
1031 | elif system == 'monoclinic' and axis=='c': |
---|
1032 | permlist = [(1,2,3),(-1,-2,3)] |
---|
1033 | elif system == 'trigonal': |
---|
1034 | permlist = [(1,2,3),(2,1,3),(-1,-2,3),(-2,-1,3)] |
---|
1035 | elif system == 'rhombohedral': |
---|
1036 | permlist = [(1,2,3),(2,3,1),(3,1,2)] |
---|
1037 | else: |
---|
1038 | permlist = [(1,2,3)] |
---|
1039 | |
---|
1040 | hklref = list(hklref) |
---|
1041 | for perm in permlist: |
---|
1042 | hkl = [abs(i) * hklin[abs(i)-1] / i for i in perm] |
---|
1043 | if hkl == hklref: return True |
---|
1044 | if [-i for i in hkl] == hklref: return True |
---|
1045 | return False |
---|
1046 | |
---|
1047 | for key in sgtbxlattinp.sgtbx8: |
---|
1048 | spdict = spc.SpcGroup(key)[1] |
---|
1049 | cell = sgtbxlattinp.sgtbx8[key][0] |
---|
1050 | center = spdict['SGLatt'] |
---|
1051 | Laue = spdict['SGLaue'] |
---|
1052 | Axis = spdict['SGUniq'] |
---|
1053 | system = spdict['SGSys'] |
---|
1054 | |
---|
1055 | g2list = GenHLaue(dmin,Laue,center,Axis,cell2A(cell)) |
---|
1056 | #if len(g2list) != len(sgtbxlattinp.sgtbx8[key][1]): |
---|
1057 | # print 'failed',key,':' ,len(g2list),'vs',len(sgtbxlattinp.sgtbx8[key][1]) |
---|
1058 | # print 'GSAS-II:' |
---|
1059 | # for h,k,l,d in g2list: print ' ',(h,k,l),d |
---|
1060 | # print 'SGTBX:' |
---|
1061 | # for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref |
---|
1062 | assert len(g2list) == len(sgtbxlattinp.sgtbx8[key][1]), ( |
---|
1063 | 'Reflection lists differ for %s' % key |
---|
1064 | ) |
---|
1065 | #match = True |
---|
1066 | for h,k,l,d in g2list: |
---|
1067 | for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: |
---|
1068 | if abs(d-dref) < derror: |
---|
1069 | if indexmatch((h,k,l,), hkllist, system, Axis): break |
---|
1070 | else: |
---|
1071 | assert 0,'No match for %s at %s (%s)' % ((h,k,l),d,key) |
---|
1072 | #match = False |
---|
1073 | #if not match: |
---|
1074 | #for hkllist,dref in sgtbxlattinp.sgtbx8[key][1]: print ' ',hkllist,dref |
---|
1075 | #print center, Laue, Axis, system |
---|
1076 | |
---|
1077 | if __name__ == '__main__': |
---|
1078 | test0() |
---|
1079 | test1() |
---|
1080 | test2() |
---|
1081 | test3() |
---|
1082 | test4() |
---|
1083 | test5() |
---|
1084 | test6() |
---|
1085 | test7() |
---|
1086 | test8() |
---|
1087 | print "OK" |
---|