[71] | 1 | '''Perform lattice-related computations''' |
---|
| 2 | |
---|
[155] | 3 | import math |
---|
[70] | 4 | import numpy as np |
---|
| 5 | import numpy.linalg as nl |
---|
| 6 | |
---|
| 7 | # trig functions in degrees |
---|
| 8 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
| 9 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
| 10 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
| 11 | atand = lambda x: 180.*np.arctan(x)/np.pi |
---|
| 12 | atan2d = lambda y,x: 180.*np.atan2(y,x)/np.pi |
---|
| 13 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
| 14 | acosd = lambda x: 180.*np.arccos(x)/np.pi |
---|
| 15 | rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p) |
---|
| 16 | |
---|
[78] | 17 | def 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) |
---|
[158] | 22 | |
---|
| 23 | def 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)]]) |
---|
[164] | 33 | |
---|
| 34 | def 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 | |
---|
[70] | 41 | def 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 |
---|
[73] | 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]]) |
---|
[70] | 51 | return g |
---|
| 52 | |
---|
| 53 | def 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 | |
---|
| 62 | def A2Gmat(A): |
---|
[73] | 63 | '''Fill reciprocal metric tensor (G) from A |
---|
[70] | 64 | returns reciprocal (G) & real (g) metric tensors (list of two 3x3 arrays) |
---|
| 65 | ''' |
---|
| 66 | G = np.zeros(shape=(3,3)) |
---|
[73] | 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]]] |
---|
[70] | 71 | g = nl.inv(G) |
---|
| 72 | return G,g |
---|
| 73 | |
---|
| 74 | def Gmat2A(G): |
---|
[73] | 75 | 'Extract A from reciprocal metric tensor (G)' |
---|
[70] | 76 | return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]] |
---|
| 77 | |
---|
| 78 | def cell2A(cell): |
---|
| 79 | G,g = cell2Gmat(cell) |
---|
| 80 | return Gmat2A(G) |
---|
| 81 | |
---|
| 82 | def 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 | |
---|
| 89 | def Gmat2cell(g): |
---|
| 90 | '''Compute lattice parameters from real metric tensor (g) |
---|
| 91 | returns tuple with a,b,c,alpha, beta, gamma (degrees) |
---|
[73] | 92 | Alternatively,compute reciprocal lattice parameters from inverse metric tensor (G) |
---|
| 93 | returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) |
---|
[70] | 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 | |
---|
| 103 | def 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 | |
---|
| 113 | def calc_rVsq(A): |
---|
[73] | 114 | 'Compute the square of the reciprocal lattice volume (V* **2) from A' |
---|
| 115 | G,g = A2Gmat(A) |
---|
| 116 | rVsq = nl.det(G) |
---|
[70] | 117 | if rVsq < 0: |
---|
| 118 | return 1 |
---|
| 119 | return rVsq |
---|
| 120 | |
---|
| 121 | def calc_rV(A): |
---|
[76] | 122 | 'Compute the reciprocal lattice volume (V*) from A' |
---|
[70] | 123 | return np.sqrt(calc_rVsq(A)) |
---|
| 124 | |
---|
| 125 | def calc_V(A): |
---|
[76] | 126 | 'Compute the real lattice volume (V) from A' |
---|
[70] | 127 | return 1./calc_rV(A) |
---|
| 128 | |
---|
| 129 | def A2invcell(A): |
---|
[76] | 130 | '''Compute reciprocal unit cell constants from A |
---|
[70] | 131 | returns tuple with a*,b*,c*,alpha*, beta*, gamma* (degrees) |
---|
| 132 | ''' |
---|
| 133 | G,g = A2Gmat(A) |
---|
| 134 | return Gmat2cell(G) |
---|
| 135 | |
---|
| 136 | def cell2AB(cell): |
---|
| 137 | '''Computes orthogonalization matrix from unit cell constants |
---|
| 138 | cell is tuple with a,b,c,alpha, beta, gamma (degrees) |
---|
[155] | 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 |
---|
[70] | 142 | ''' |
---|
[71] | 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 |
---|
[158] | 155 | |
---|
[155] | 156 | def 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 | |
---|
| 167 | def 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 | |
---|
| 178 | def 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] |
---|
[70] | 195 | |
---|
[155] | 196 | #Permutations and Combinations |
---|
| 197 | # Four routines: combinations,uniqueCombinations, selections & permutations |
---|
| 198 | #These taken from Python Cookbook, 2nd Edition. 19.15 p724-726 |
---|
| 199 | # |
---|
| 200 | def _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 |
---|
| 209 | def 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) |
---|
| 214 | def 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) |
---|
| 219 | def 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) |
---|
| 224 | def permutations(items): |
---|
| 225 | ''' take all items, order matters ''' |
---|
| 226 | return combinations(items, len(items)) |
---|
| 227 | |
---|
[78] | 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 | |
---|
| 232 | def 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 | |
---|
[166] | 236 | def calc_rDsq2(H,G): |
---|
| 237 | return np.inner(H,np.inner(G,H)) |
---|
| 238 | |
---|
[78] | 239 | def 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 | |
---|
| 244 | def 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 | |
---|
| 254 | def 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 | |
---|
| 274 | def 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 | |
---|
[108] | 282 | def 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 | |
---|
| 289 | def 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 | |
---|
[78] | 303 | def 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 | |
---|
[93] | 320 | def 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 | |
---|
[78] | 358 | def GenHBravais(dmin,Bravais,A): |
---|
| 359 | '''Generate the positionally unique powder diffraction reflections |
---|
[93] | 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 |
---|
[78] | 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]) |
---|
[70] | 458 | |
---|
[78] | 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 | |
---|
| 471 | def GenHLaue(dmin,Laue,Cent,Axis,A): |
---|
| 472 | '''Generate the crystallographically unique powder diffraction reflections |
---|
| 473 | for a lattice and Bravais type |
---|
| 474 | ''' |
---|
| 475 | # dmin - minimum d-spacing |
---|
[103] | 476 | # Laue - Laue group symbol = '-1','2/m','mmm','4/m','6/m','4/mmm','6/mmm', |
---|
[78] | 477 | # '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m' |
---|
| 478 | # Cent - lattice centering = 'P','A','B','C','I','F' |
---|
| 479 | # Axis - 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 |
---|
[108] | 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)) |
---|
[109] | 492 | #print Hmax,aHx,cHx |
---|
[108] | 493 | else: # all others |
---|
| 494 | Hmax = MaxIndex(dmin,A) |
---|
| 495 | |
---|
[78] | 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(Cent,[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 |
---|
[102] | 510 | axisnum = 1 + ['a','b','c'].index(Axis) |
---|
| 511 | Hmax = SwapIndx(axisnum,Hmax) |
---|
[78] | 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): |
---|
[102] | 517 | [h,k,l] = SwapIndx(-axisnum,[h,k,l]) |
---|
[78] | 518 | H = [] |
---|
| 519 | if CentCheck(Cent,[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)]) |
---|
[102] | 524 | [h,k,l] = SwapIndx(axisnum,[h,k,l]) |
---|
[78] | 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 |
---|
[102] | 530 | for k in range(kmin,Hmax[1]+1): |
---|
[78] | 531 | H = [] |
---|
| 532 | if CentCheck(Cent,[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(Cent,[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 |
---|
[109] | 554 | kmin = -int((h-1.)/2.) |
---|
[78] | 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(Cent,[h,k,l]): H=[h,k,l] |
---|
[108] | 563 | if Laue in ['3R','3mR']: |
---|
| 564 | H = Hx2Rh(H) |
---|
[78] | 565 | if H: |
---|
| 566 | rdsq = calc_rDsq(H,A) |
---|
| 567 | if 0 < rdsq <= dminsq: |
---|
[109] | 568 | HKL.append([H[0],H[1],H[2],1/math.sqrt(rdsq)]) |
---|
[78] | 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(Cent,[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 | |
---|
[70] | 586 | # output from uctbx computed on platform darwin on 2010-05-28 |
---|
[93] | 587 | NeedTestData = True |
---|
| 588 | def TestData(): |
---|
| 589 | array = np.array |
---|
| 590 | global NeedTestData |
---|
| 591 | NeedTestData = False |
---|
| 592 | global CellTestData |
---|
| 593 | CellTestData = [ |
---|
[70] | 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], |
---|
[93] | 615 | ] |
---|
| 616 | global CoordTestData |
---|
| 617 | CoordTestData = [ |
---|
[71] | 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 | ] |
---|
[70] | 649 | |
---|
| 650 | def test0(): |
---|
[93] | 651 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
| 662 | def test1(): |
---|
[93] | 663 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
| 670 | def test2(): |
---|
[93] | 671 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
| 678 | def test3(): |
---|
[93] | 679 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
| 686 | def test4(): |
---|
[93] | 687 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
| 693 | def test5(): |
---|
[93] | 694 | if NeedTestData: TestData() |
---|
[70] | 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 | |
---|
[71] | 700 | def test6(): |
---|
[93] | 701 | if NeedTestData: TestData() |
---|
[71] | 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 |
---|
[108] | 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 |
---|
[71] | 714 | |
---|
[102] | 715 | # test GetBraviasNum(...) and GenHBravais(...) |
---|
| 716 | def 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 | |
---|
| 763 | def 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)] |
---|
[108] | 782 | elif system == 'rhombohedral': |
---|
[162] | 783 | permlist = [(1,2,3),(2,3,1),(3,1,2)] |
---|
[102] | 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 | |
---|
[70] | 824 | if __name__ == '__main__': |
---|
| 825 | test0() |
---|
| 826 | test1() |
---|
| 827 | test2() |
---|
| 828 | test3() |
---|
| 829 | test4() |
---|
| 830 | test5() |
---|
[71] | 831 | test6() |
---|
[102] | 832 | test7() |
---|
| 833 | test8() |
---|
[70] | 834 | print "OK" |
---|