Changeset 380 for trunk/GSASIIspc.py
 Timestamp:
 Oct 3, 2011 3:43:49 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIspc.py
r372 r380 322 322 h,k,l,f = Uniq 323 323 Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq])) 324 phi = f[:Nuniq] 324 325 Uniq = np.array(Uniq) 325 326 326 return iabsnt,2*mulp,Uniq #include Friedel pairs in powder mulp 327 328 329 def GenHKL(HKL,SGData,Friedel=False): #gives wrong answers! 330 ''' 331 Generate set of equivalent reflections 332 input: 333 HKL  [h,k,l] 334 SGData  space group data obtained from SpcGroup 335 Friedel = True to retain Friedel pairs for centrosymmetric case 336 returns: 337 iabsnt = True is reflection is forbidden by symmetry 338 mulp = reflection multiplicity including Fridel pairs 339 Uniq = numpy array of equivalent hkl in descending order of h,k,l 340 Phs = numpy array of corresponding phase factors (multiples of 2pi) 341 342 NB: only works on one HKL at a time  343 it can be made to do an array of HKLs  not any faster! 344 ''' 345 iabsnt = False 346 H = np.array(HKL) 347 Inv = SGData['SGInv'] 348 Cen = SGData['SGCen'][1:] #skip P 349 # do cell centering extinction check 350 for cen in Cen: 351 iabsnt = bool(int(round(np.sum(cen*H*12)))%12) 352 # get operators & expand if centrosymmetric 353 Ops = SGData['SGOps'] 354 opM = np.array([op[0] for op in Ops]) 355 opT = np.array([op[1] for op in Ops]) 356 if Inv: 357 opM = np.concatenate((opM,opM)) 358 opT = np.concatenate((opT,opT)) 359 # generate full hkl table and phase factors; find duplicates for multiplicity 360 eqH = np.inner(opM,H) 361 HT = np.dot(opT,H)%1. 362 dup = len(ma.nonzero([ma.allclose(H,hkl,atol=0.001) for hkl in eqH])[0]) 363 mulp = len(eqH)/dup 364 # generate unique reflection set (with or without Friedel pairs) 365 HPT = [tuple(hp) for hp in np.column_stack((eqH,HT))] 366 HPT.sort() 367 UPT = HPT[::dup] 368 UPT.reverse() 369 if Inv and not Friedel: 370 UPT = UPT[:len(UPT)/2] 371 Uniq = np.split(UPT,[3,4],axis=1) 372 Phs = Uniq[1].T[0] 373 Uniq = Uniq[0] 374 # determine if extinct 375 HPT = np.split(HPT,[3,4],axis=1)[1].T[0] 376 HPT = np.reshape(HPT,(mulp,1)).T 377 HPT = np.around(HPTHPT[0],4)%1. 378 iabsnt = np.any(HPT) 379 return iabsnt,mulp,Uniq 380 381 #def GenHKLs(H,SGData,Friedel=False): 382 # ''' 383 # Generate set of equivalent reflections 384 # input: 385 # H  np.array of [h,k,l] assumed to exclude lattice centering extinctions 386 # SGData  space group data obtained from SpcGroup 387 # Friedel = True to retain Friedel pairs for centrosymmetric case 388 # returns: 389 # iabsnt = True is reflection is forbidden by symmetry 390 # mulp = reflection multiplicity including Fridel pairs 391 # Uniq = numpy array of equivalent hkl in descending order of h,k,l 392 # Phs = numpy array of corresponding phase factors (multiples of 2pi) 393 # this is not any faster than GenHKL above  oh well 394 # ''' 395 # H = np.array(H,dtype=float) 396 # Inv = SGData['SGInv'] 397 # Ops = SGData['SGOps'] 398 # opM = np.array([op[0] for op in Ops]) 399 # opT = np.array([op[1] for op in Ops]) 400 # if Inv: 401 # opM = np.concatenate((opM,opM)) 402 # opT = np.concatenate((opT,opT)) 403 # # generate full hkl table and phase factors; find duplicates for multiplicity 404 # eqH = np.swapaxes(np.swapaxes(np.inner(opM,H),0,2),1,2) 405 # HeqH = zip(H,eqH) 406 # dup = np.array([sum([ma.allclose(h,hkl,atol=0.001) for hkl in HKL]) for h,HKL in HeqH]) 407 # mulp = len(eqH[0])/dup 408 # # generate unique reflection set (with or without Friedel pairs) 409 # HT = np.dot(H,opT.T)%1. 410 # HPT = zip(mulp,dup,[[tuple(hp) for hp in HP] for HP in np.dstack((eqH,HT))]) 411 # Uniq = [] 412 # Phs = [] 413 # iabsnt = [] 414 # for mp,dp,hpt in HPT: 415 # hpt.sort() 416 # upt = hpt[::dp] 417 # upt.reverse() 418 # if Inv and not Friedel: 419 # upt = upt[:len(upt)/2] 420 # uniq = np.split(upt,[3,4],axis=1) 421 # Phs.append(uniq[1].T[0]) 422 # Uniq.append(uniq[0]) 423 # # determine if extinct 424 # hpt = np.split(hpt,[3,4],axis=1)[1].T[0] 425 # hpt = np.reshape(hpt,(mp,1)).T 426 # hpt = np.around(hpthpt[0],4)%1. 427 # iabsnt.append(np.any(hpt)) 428 # return iabsnt,mulp,Uniq,Phs 429 327 return iabsnt,2*mulp,Uniq,phi #include Friedel pairs in powder mulp 328 430 329 def GetOprPtrName(key): 431 330 OprPtrName = {
Note: See TracChangeset
for help on using the changeset viewer.