Changeset 311 for trunk/GSASIIindex.py


Ignore:
Timestamp:
Jun 27, 2011 10:27:29 AM (10 years ago)
Author:
vondreele
Message:

finish new indexing refinement stuff
more on texture plotting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIindex.py

    r307 r311  
    245245            return [values[0],values[1],values[2],0,values[3],0]
    246246        else:
    247             return values
     247            return list(values)
    248248           
    249249    def A2values(ibrav,A):
     
    258258        else:
    259259            return A
     260           
     261    def Jacobian(values,ibrav,d,H,Pwr):
     262# derivatives of 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]
     263        if ibrav in [0,1,2]:            #m3m
     264            return [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]*d**Pwr
     265        elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
     266            return [H[0]*H[0]+H[1]*H[0]+H[1]*H[1],H[2]*H[2]]*d**Pwr
     267        elif ibrav in [5,6]:            #4/mmm
     268            return [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]*d**Pwr
     269        elif ibrav in [7,8,9,10]:       #mmm
     270            return [H[0]*H[0],H[1]*H[1],H[2]*H[2]]*d**Pwr
     271        elif ibrav in [11,12]:          #2/m
     272            return [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]*d**Pwr
     273        else:                           #-1
     274            return [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[1]*H[0],H[0]*H[2],H[1]*H[2]]*d**Pwr
    260275   
    261276    def errFit(values,ibrav,d,H,Pwr):
     
    266281   
    267282    Peaks = np.array(peaks).T
    268     values = A2values(ibrav,A)   
    269     result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),full_output=True)
     283   
     284    values = A2values(ibrav,A)
     285#    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),
     286#        Dfun=Jacobian,col_deriv=True,full_output=True)
     287    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),
     288        full_output=True,factor=0.1)
    270289    A = Values2A(ibrav,result[0])
    271     return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A
     290    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
    272291           
    273292def FitHKLZ(ibrav,peaks,Z,A):
     
    314333    return peaks[0][7]
    315334   
    316 def refinePeaks(peaks,ibrav,A):
     335def refinePeaks(peaks,ibrav,A,Zero):
     336    #Zero is list(zero value, flag)
    317337    dmin = getDmin(peaks)
    318338    smin = 1.0e10
    319     pwr = 3
     339    pwr = 4
    320340    maxTries = 10
    321     if ibrav == 13:
     341    if ibrav == 13: #-1 - triclinic
    322342        pwr = 4
    323343        maxTries = 10
     
    325345    tries = 0
    326346    HKL = G2lat.GenHBravais(dmin,ibrav,A)
    327     while IndexPeaks(peaks,HKL):
     347    while len(HKL) > 2 and IndexPeaks(peaks,HKL):
    328348        Pwr = pwr - (tries % 2)
    329349        HKL = []
    330350        tries += 1
    331351        osmin = smin
    332         oldA = A
    333         OK,smin,A = FitHKL(ibrav,peaks,A,Pwr)
    334         if min(A[:3]) <= 0:
     352        oldA = A[:]
     353        Vold = G2lat.calc_V(oldA)
     354        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
     355        Vnew = G2lat.calc_V(A)
     356        if Vnew > 2.0*Vold or Vnew < 2.:
     357            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
     358            OK = False
     359            continue
     360        try:
     361            HKL = G2lat.GenHBravais(dmin,ibrav,A)
     362        except FloatingPointError:
    335363            A = oldA
    336364            OK = False
    337365            break
    338         if OK:
    339             HKL = G2lat.GenHBravais(dmin,ibrav,A)
    340366        if len(HKL) == 0: break                         #absurd cell obtained!
    341367        rat = (osmin-smin)/smin
     
    343369        if tries > maxTries: break
    344370    if OK:
    345         OK,smin,A = FitHKL(ibrav,peaks,A,2)
     371        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
    346372        Peaks = np.array(peaks).T
    347373        H = Peaks[4:7]
     
    350376       
    351377    M20,X20 = calc_M20(peaks,HKL)
    352     return len(HKL),M20,X20,A
     378    return len(HKL),M20,X20,A,Zero
    353379       
    354380def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
     
    371397    while tries < Ntries:
    372398        if A:
    373             Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
    374             if ibrav in [11,12,13]:
    375                 Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
     399            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
     400            if ibrav in [11,12,13]:         #monoclinic & triclinic
     401                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
    376402        else:
    377             Anew = ranAbyV(ibrav,amin,amax,V1)
    378         HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
     403            Abeg = ranAbyV(ibrav,amin,amax,V1)
     404        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
    379405       
    380406        if IndexPeaks(peaks,HKL) and len(HKL) > mHKL[ibrav]:
    381             Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    382             Asave.append([calc_M20(peaks,HKL),Anew[:]])
     407            Lhkl,M20,X20,Aref,Zero = refinePeaks(peaks,ibrav,Abeg,[0,False])
     408            Asave.append([calc_M20(peaks,HKL),Aref[:]])
    383409            if ibrav == 9:                          #C-centered orthorhombic
    384410                for i in range(2):
    385                     Anew = rotOrthoA(Anew[:])
    386                     Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    387                     HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
     411                    Abeg = rotOrthoA(Abeg[:])
     412                    Lhkl,M20,X20,Aref,Zero = refinePeaks(peaks,ibrav,Abeg,[0,False])
     413                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
    388414                    IndexPeaks(peaks,HKL)
    389                     Asave.append([calc_M20(peaks,HKL),Anew[:]])
     415                    Asave.append([calc_M20(peaks,HKL),Aref[:]])
    390416            elif ibrav == 11:                      #C-centered monoclinic
    391                 Anew = swapMonoA(Anew[:])
    392                 Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    393                 HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
     417                Abeg = swapMonoA(Abeg[:])
     418                Lhkl,M20,X20,Aref,Zero = refinePeaks(peaks,ibrav,Abeg,[0,False])
     419                HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
    394420                IndexPeaks(peaks,HKL)
    395                 Asave.append([calc_M20(peaks,HKL),Anew[:]])
     421                Asave.append([calc_M20(peaks,HKL),Aref[:]])
    396422        else:
    397423            break
     
    406432    X = sortM20(Asave)
    407433    if X:
    408         Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
    409         return GoOn,Lhkl,M20,X20,X[0][1]
    410     else:
    411         return GoOn,0,0,0,Anew
     434        Lhkl,M20,X20,A,Zero = refinePeaks(peaks,ibrav,X[0][1],[0,False])
     435        return GoOn,Lhkl,M20,X20,A
     436       
     437    else:
     438        return GoOn,0,0,0,0
    412439       
    413440def monoCellReduce(ibrav,A):
     
    537564        return False,0,0
    538565       
     566       
     567NeedTestData = True
     568def TestData():
     569    array = np.array
     570    global NeedTestData
     571    NeedTestData = False
     572    global TestData
     573    TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
     574        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436],
     575        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777],
     576        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718],
     577        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106],
     578        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946],
     579        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333],
     580        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
     581        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592],
     582        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413],
     583        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829],
     584        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108],
     585        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944],
     586        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193],
     587        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
     588        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389],
     589        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764],
     590        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602],
     591        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214],
     592        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812],
     593        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784],
     594        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369],
     595        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648],
     596        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553],
     597        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181],
     598        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037],
     599        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951],
     600        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487],
     601        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
     602        ]
     603    global TestData2
     604    TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
     605        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801],
     606        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881],
     607        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0],
     608        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192],
     609        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0],
     610        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479],
     611        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754],
     612        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0],
     613        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0],
     614        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534],
     615        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
     616        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0],
     617        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0],
     618        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0],
     619        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441],
     620        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0],
     621        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986],
     622        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005],
     623        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166],
     624        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428],
     625        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503],
     626        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0],
     627        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509],
     628        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121],
     629        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0],
     630        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0],
     631        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0],
     632        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
     633        ]
     634#
     635def test0():
     636    if NeedTestData: TestData()
     637    msg = 'test FitHKL'
     638    ibrav,cell,bestcell,Pwr,peaks = TestData
     639    print 'best cell:',bestcell
     640    print 'old cell:',cell
     641    Peaks = np.array(peaks)
     642    HKL = Peaks[4:7]
     643    print calc_M20(peaks,HKL)
     644    A = G2lat.cell2A(cell)
     645    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
     646    print 'new cell:',G2lat.A2cell(A)   
     647    print 'x:',result[0]
     648    print 'cov_x:',result[1]
     649    print 'infodict:'
     650    for item in result[2]:
     651        print item,result[2][item]
     652    print 'msg:',result[3]
     653    print 'ier:',result[4]
     654    result = refinePeaks(peaks,ibrav,A,[0,False])
     655    N,M20,X20,A,Zero = result
     656    print 'refinePeaks:',N,M20,X20,G2lat.A2cell(A)
     657    print 'compare bestcell:',bestcell
     658#
     659def test1():
     660    if NeedTestData: TestData()
     661    msg = 'test FitHKL'
     662    ibrav,A,Pwr,peaks = TestData2
     663    print 'bad cell:',G2lat.A2cell(A)
     664    print 'FitHKL'
     665    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
     666    result = refinePeaks(peaks,ibrav,A,[0,False])
     667    N,M20,X20,A,Zero = result
     668    print 'refinePeaks:',N,M20,X20,A
     669#    Peaks = np.array(peaks)
     670#    HKL = Peaks[4:7]
     671#    print calc_M20(peaks,HKL)
     672#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
     673#    print 'new cell:',G2lat.A2cell(A)   
     674#    print 'x:',result[0]
     675#    print 'cov_x:',result[1]
     676#    print 'infodict:'
     677#    for item in result[2]:
     678#        print item,result[2][item]
     679#    print 'msg:',result[3]
     680#    print 'ier:',result[4]
     681   
     682#
     683if __name__ == '__main__':
     684    test0()
     685    test1()
     686#    test2()
     687#    test3()
     688#    test4()
     689#    test5()
     690#    test6()
     691#    test7()
     692#    test8()
     693    print "OK"
Note: See TracChangeset for help on using the changeset viewer.