Changeset 380 for trunk/GSASIIspc.py


Ignore:
Timestamp:
Oct 3, 2011 3:43:49 PM (11 years ago)
Author:
vondreele
Message:

completed Rietveld refinement version

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIspc.py

    r372 r380  
    322322    h,k,l,f = Uniq
    323323    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
     324    phi = f[:Nuniq]
    324325    Uniq = np.array(Uniq)
    325326   
    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(HPT-HPT[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(hpt-hpt[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                                 
    430329def GetOprPtrName(key):           
    431330    OprPtrName = {
Note: See TracChangeset for help on using the changeset viewer.