Changeset 303 for trunk/GSASIIindex.py


Ignore:
Timestamp:
Jun 20, 2011 2:37:07 PM (11 years ago)
Author:
vondreele
Message:

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIindex.py

    r302 r303  
    1616import pypowder as pyp              #assumes path has been amended to include correctr bin directory
    1717import GSASIIlattice as G2lat
     18import scipy.optimize as so
    1819
    1920# trig functions in degrees
     
    194195def IndexPeaks(peaks,HKL):
    195196    import bisect
    196     hklds = [1000.0]                                    #make bounded list of available d-spacings
    197197    N = len(HKL)
    198198    if N == 0: return False
    199     for hkl in HKL:
    200         hklds.append(hkl[3])
    201     hklds.append(0.0)
     199    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
    202200    hklds.sort()                                        # ascending sort - upper bound at end
    203201    hklmax = [0,0,0]
     
    215213                dnew = min(dm,dp)
    216214                if dold > dnew:                             # new better - zero out old
    217                     for j in range(3):
    218                         opeak[j+4] = 0
     215                    opeak[4:7] = [0,0,0]
    219216                    opeak[8] = 0.
    220217                else:                                       # old better - do nothing
    221218                    continue               
    222219            hkl[4] = ipk
    223             for j in range(3):
    224                 peak[j+4] = hkl[j]
     220            peak[4:7] = hkl[:3]
    225221            peak[8] = hkl[3]                                # fill in d-calc
    226222    for peak in peaks:
     
    235231    else:
    236232        return False
    237    
    238 def FitHKL(ibrav,peaks,A,wtP):
    239     def ShiftTest(a,b):
    240         if b < -0.1*a:
    241             b = -0.0001*a
    242         return b
    243     smin = 0.
    244     first = True
    245     for peak in peaks:
    246         if peak[2] and peak[3]:
    247             h,k,l = H = peak[4:7]
    248             Qo = 1./peak[7]**2
    249             Qc = G2lat.calc_rDsq(H,A)
    250             try:
    251                 peak[8] = 1./math.sqrt(Qc)
    252             except:
    253                 print G2lat.A2invcell(A)
    254             delt = Qo-Qc
    255             smin += delt**2
    256             dp = []
    257             if ibrav in [0,1,2]:            #m3m
    258                 dp.append(h*h+k*k+l*l)
    259             elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    260                 dp.append(h*h+k*k+h*k)
    261                 dp.append(l*l)
    262             elif ibrav in [5,6]:            #4/mmm
    263                 dp.append(h*h+k*k)
    264                 dp.append(l*l)
    265             elif ibrav in [7,8,9,10]:       #mmm
    266                 dp.append(h*h)
    267                 dp.append(k*k)
    268                 dp.append(l*l)
    269             elif ibrav in [11,12]:          #2/m
    270                 dp.append(h*h)
    271                 dp.append(k*k)
    272                 dp.append(l*l)
    273                 dp.append(h*l)
    274             else:                           #1
    275 #    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
    276                 dp.append(h*h)
    277                 dp.append(k*k)
    278                 dp.append(l*l)
    279                 dp.append(h*k)
    280                 dp.append(h*l)
    281                 dp.append(k*l)
    282             if first:
    283                 first = False
    284                 M = len(dp)
    285                 B = np.zeros(shape=(M,M))
    286                 V = np.zeros(shape=(M))
    287             dwt = peak[7]**wtP
    288             B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
    289     if nl.det(B) > 0.0:
    290         try:
    291             b = nl.solve(B,V)
    292             B = nl.inv(B)
    293             sig = np.diag(B)
    294         except SingularMatrix:
    295             return False,0
    296         if ibrav in [0,1,2]:            #m3m
    297             A[0] += ShiftTest(A[0],b[0])
    298             A[1] = A[2] = A[0]
    299         elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    300             A[0] += ShiftTest(A[0],b[0])
    301             A[1] = A[3] = A[0]
    302             A[2] += ShiftTest(A[2],b[1])
    303         elif ibrav in [5,6]:            #4/mmm
    304             A[0] += ShiftTest(A[0],b[0])
    305             A[1] = A[0]
    306             A[2] += ShiftTest(A[2],b[1])
    307         elif ibrav in [7,8,9,10]:       #mmm
    308             for i in range(3):
    309                 A[i] += ShiftTest(A[i],b[i])
    310         elif ibrav in [11,12]:          #2/m
    311             for i in range(3):
    312                 A[i] += ShiftTest(A[i],b[i])
    313             A[4] += ShiftTest(A[4],b[3])
    314             A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
    315         else:                           #1
    316             oldV = math.sqrt(1./G2lat.calc_rVsq(A))
    317             oldA = A[:]
    318             for i in range(6):
    319                 A[i] += b[i]*0.2
    320             A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
    321             A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
    322             A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
    323             ratio = math.sqrt(1./G2lat.calc_rVsq(A))/oldV
    324             if 0.9 > ratio or ratio > 1.1:
    325                 A = oldA
    326 #    else:
    327 #        print B
    328 #        print V
    329 #        for peak in peaks:
    330 #            print peak
    331     return True,smin
    332        
     233
     234def FitHKL(ibrav,peaks,A,Pwr):
     235   
     236    def Values2A(ibrav,values):
     237        if ibrav in [0,1,2]:
     238            return [values[0],values[0],values[0],0,0,0]
     239        elif ibrav in [3,4]:
     240            return [values[0],values[0],values[1],values[0],0,0]
     241        elif ibrav in [5,6]:
     242            return [values[0],values[0],values[1],0,0,0]
     243        elif ibrav in [7,8,9,10]:
     244            return [values[0],values[1],values[2],0,0,0]
     245        elif ibrav in [11,12]:
     246            return [values[0],values[1],values[2],0,values[3],0]
     247        else:
     248            return values
     249           
     250    def A2values(ibrav,A):
     251        if ibrav in [0,1,2]:
     252            return [A[0],]
     253        elif ibrav in [3,4,5,6]:
     254            return [A[0],A[2]]
     255        elif ibrav in [7,8,9,10]:
     256            return [A[0],A[1],A[2]]
     257        elif ibrav in [11,12]:
     258            return [A[0],A[1],A[2],A[4]]
     259        else:
     260            return A
     261   
     262    def errFit(values,ibrav,d,H,Pwr):
     263        A = Values2A(ibrav,values)
     264        Qo = 1./d**2
     265        Qc = G2lat.calc_rDsq(H,A)
     266        return (Qo-Qc)*d**Pwr
     267   
     268    Peaks = np.array(peaks).T
     269    values = A2values(ibrav,A)   
     270    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),full_output=True)
     271    A = Values2A(ibrav,result[0])
     272    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A
     273           
    333274def FitHKLZ(ibrav,peaks,Z,A):
    334275    return A,Z
     
    377318    dmin = getDmin(peaks)
    378319    smin = 1.0e10
    379     pwr = 6
     320    pwr = 3
    380321    maxTries = 10
    381322    if ibrav == 13:
     
    386327    HKL = G2lat.GenHBravais(dmin,ibrav,A)
    387328    while IndexPeaks(peaks,HKL):
    388         Pwr = pwr - 2*(tries % 2)
     329        Pwr = pwr - (tries % 2)
    389330        HKL = []
    390331        tries += 1
    391332        osmin = smin
    392333        oldA = A
    393         OK,smin = FitHKL(ibrav,peaks,A,Pwr)
    394         for a in A[:3]:
    395             if a < 0:
    396                 A = oldA
    397                 OK = False
    398                 break
     334        OK,smin,A = FitHKL(ibrav,peaks,A,Pwr)
     335        if min(A[:3]) <= 0:
     336            A = oldA
     337            OK = False
     338            break
    399339        if OK:
    400340            HKL = G2lat.GenHBravais(dmin,ibrav,A)
     
    404344        if tries > maxTries: break
    405345    if OK:
    406         OK,smin = FitHKL(ibrav,peaks,A,4)
     346        OK,smin,A = FitHKL(ibrav,peaks,A,2)
     347        Peaks = np.array(peaks).T
     348        H = Peaks[4:7]
     349        Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
     350        peaks = Peaks.T
     351       
    407352    M20,X20 = calc_M20(peaks,HKL)
    408     return len(HKL),M20,X20
     353    return len(HKL),M20,X20,A
    409354       
    410355def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
     
    433378            Anew = ranAbyV(ibrav,amin,amax,V1)
    434379        HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    435         if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
    436             Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     380       
     381        if IndexPeaks(peaks,HKL) and len(HKL) > mHKL[ibrav]:
     382            Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    437383            Asave.append([calc_M20(peaks,HKL),Anew[:]])
    438384            if ibrav == 9:                          #C-centered orthorhombic
    439385                for i in range(2):
    440386                    Anew = rotOrthoA(Anew[:])
    441                     Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     387                    Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    442388                    HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    443389                    IndexPeaks(peaks,HKL)
     
    445391            elif ibrav == 11:                      #C-centered monoclinic
    446392                Anew = swapMonoA(Anew[:])
    447                 Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     393                Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    448394                HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    449395                IndexPeaks(peaks,HKL)
     
    461407    X = sortM20(Asave)
    462408    if X:
    463         Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
     409        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
    464410        return GoOn,Lhkl,M20,X20,X[0][1]
    465411    else:
    466         return GoOn,0,0,0,0
     412        return GoOn,0,0,0,Anew
    467413       
    468414def monoCellReduce(ibrav,A):
     
    489435def DoIndexPeaks(peaks,inst,controls,bravais):
    490436   
    491     def peakDspace(peaks,A):
    492         for peak in peaks:
    493             if peak[3]:
    494                 dsq = calc_rDsq(peak[4:7],A)
    495                 if dsq > 0:
    496                     peak[8] = 1./math.sqrt(dsq)
    497         return
    498437    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
    499438    amin = 2.5
     
    541480                            if ibrav > 2:
    542481                                if not N2:
    543                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
     482                                    A = []
     483                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
    544484                                if A:
    545485                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
Note: See TracChangeset for help on using the changeset viewer.