Changeset 885


Ignore:
Timestamp:
Apr 13, 2013 1:38:32 PM (11 years ago)
Author:
vondreele
Message:

add 'Global' to constraints
remove a stray debug print
fix U6toUij so it returns a numpy array
fixes to Q2AV & Q2AVdeg
RB changes - get derivatives in

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r871 r885  
    387387            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
    388388        IFD[Tag] = [Type,nVal,Value]
    389         print Tag,IFD[Tag]
     389#        print Tag,IFD[Tag]
    390390    sizexy = [IFD[256][2][0],IFD[257][2][0]]
    391391    [nx,ny] = sizexy

  • trunk/GSASIIconstrGUI.py

    r883 r885  
    9797    '''
    9898    if not data:
    99         data.update({'Hist':[],'HAP':[],'Phase':[]})       #empty dict - fill it
     99        data.update({'Hist':[],'HAP':[],'Phase':[],'Global':[]})       #empty dict - fill it
     100    if 'Global' not in data:                                            #patch
     101        data['Global'] = []
    100102    Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree()
    101103    rigidbodyDict = G2frame.PatternTree.GetItemPyData(   
    … …  
    103105    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    104106    rbVary,rbDict = G2str.GetRigidBodyModels(rigidbodyDict,Print=False)
     107    globalList = rbDict.keys()
     108    globalList.sort()
    105109    AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
    106110    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,rbIds=rbIds,Print=False)
    … …  
    142146        hlegend = '\n In :h:name'
    143147        phlegend = '\n In p:h:name'
     148        glegend = '\n In ::name'
    144149        for phase in Phases:
    145150            plegend += '\n p:: = '+str(Phases[phase]['pId'])+':: for '+phase
    … …  
    160165                break
    161166            count += 1
    162         return plegend,hlegend,phlegend
     167        return plegend,hlegend,phlegend,glegend
    163168       
    164169    def FindEquivVarb(name,nameList):
    … …  
    171176        elif 'AU' in name:
    172177            namelist = ['AUiso','AU11','AU22','AU33','AU12','AU13','AU23']
     178        elif 'RB' in name:
     179            rbfx = 'RB'+items[2][2]
     180            if 'T' in name and 'Tr' not in name:
     181                namelist = [rbfx+'T11',rbfx+'T22',rbfx+'T33',rbfx+'T12',rbfx+'T13',rbfx+'T23']
     182            if 'L' in name:
     183                namelist = [rbfx+'L11',rbfx+'L22',rbfx+'L33',rbfx+'L12',rbfx+'L13',rbfx+'L23']
     184            if 'S' in name:
     185                namelist = [rbfx+'S12',rbfx+'S13',rbfx+'S21',rbfx+'S23',rbfx+'S31',rbfx+'S32',rbfx+'SAA',rbfx+'SBB']
     186            if 'U' in name:
     187                namelist = [rbfx+'U',]
    173188        for item in nameList:
    174189            keys = item.split(':')
    … …  
    558573        Size[1] = min(Size[1],250)
    559574        G2frame.dataFrame.setSizePosLeft(Size)
     575
     576    def UpdateGlobalConstr():
     577        '''Responds to press on Global Constraint tab,
     578        shows constraints in data window'''
     579        GlobalConstr.DestroyChildren()
     580        GlobalDisplay = wx.Panel(GlobalConstr)
     581        GlobalSizer = wx.BoxSizer(wx.VERTICAL)
     582        GlobalSizer.Add((5,5),0)       
     583        GlobalSizer.Add(ConstSizer('Global',GlobalDisplay))
     584        GlobalDisplay.SetSizer(GlobalSizer,True)
     585        Size = GlobalSizer.GetMinSize()
     586        Size[0] += 40
     587        Size[1] = max(Size[1],250) + 20
     588        GlobalDisplay.SetSize(Size)
     589        GlobalConstr.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
     590        Size[1] = min(Size[1],250)
     591        G2frame.dataFrame.setSizePosLeft(Size)
    560592   
    561593    def OnPageChanged(event):
    … …  
    575607            G2frame.Page = [page,'phs']
    576608            UpdatePhaseConstr()
     609        elif text == 'Global constraints':
     610            G2frame.Page = [page,'glb']
     611            UpdateGlobalConstr()
     612           
    577613
    578614    def SetStatusLine(text):
    579615        Status.SetStatusText(text)                                     
    580616       
    581     plegend,hlegend,phlegend = GetPHlegends(Phases,Histograms)
     617    plegend,hlegend,phlegend,glegend = GetPHlegends(Phases,Histograms)
    582618    scope = {'hst':['Histogram contraints:',hlegend,histList,'Hist',UpdateHistConstr],
    583619        'hap':['Histogram * Phase contraints:',phlegend,hapList,'HAP',UpdateHAPConstr],
    584         'phs':['Phase contraints:',plegend,phaseList,'Phase',UpdatePhaseConstr]}
     620        'phs':['Phase contraints:',plegend,phaseList,'Phase',UpdatePhaseConstr],
     621        'glb':['Global constraints:',glegend,globalList,'Global',UpdateGlobalConstr]}
    585622    if G2frame.dataDisplay:
    586623        G2frame.dataDisplay.Destroy()
    … …  
    604641    HistConstr = wx.ScrolledWindow(G2frame.dataDisplay)
    605642    G2frame.dataDisplay.AddPage(HistConstr,'Histogram constraints')
     643    GlobalConstr = wx.ScrolledWindow(G2frame.dataDisplay)
     644    G2frame.dataDisplay.AddPage(GlobalConstr,'Global constraints')   
    606645    UpdatePhaseConstr()
    607646

  • trunk/GSASIIlattice.py

    r825 r885  
    240240        Uij - numpy [3][3] array of uij
    241241    """
    242     U = np.zeros(shape=(3,3))
    243     U = [
     242    U = np.array([
    244243        [U6[0],  U6[3],  U6[4]],
    245244        [U6[3],  U6[1],  U6[5]],
    246         [U6[4],  U6[5],  U6[2]]]
     245        [U6[4],  U6[5],  U6[2]]])
    247246    return U
    248247
    … …  
    251250    NB: there is a non numpy version in GSASIIspc: Uij2U
    252251    """
    253     U6 = np.zeros(6)
    254     U6 = [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
     252    U6 = np.array([U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]])
    255253    return U6
    256254

  • trunk/GSASIImath.py

    r882 r885  
    268268    '''
    269269    TLStype,TLS = RBObj['ThermalMotion'][:2]
    270     Tmat = np.zeros((3,3))
    271     Lmat = np.zeros((3,3))
    272     Smat = np.zeros((3,3))
     270    T = np.zeros(6)
     271    L = np.zeros(6)
     272    S = np.zeros(8)
    273273    if 'T' in TLStype:
    274274        T = TLS[:6]
    275 #        Tmat = G2lat.U6toUij(T)
    276275    if 'L' in TLStype:
    277276        L = np.array(TLS[6:12])*(np.pi/180.)**2
    278 #        Lmat = np.array(G2lat.U6toUij(L))
    279277    if 'S' in TLStype:
    280278        S = np.array(TLS[12:])*(np.pi/180.)
    281 #        Smat = np.array([[S[6],S[0],S[1]],[S[2],S[7],S[3]],[S[4],S[5],0]])
    282279    g = np.inner(Bmat,Bmat.T)
    283280    gvec = 1./np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
    … …  
    287284    for X in Cart:
    288285        X = prodQVQ(Q,X)
    289 #        Axyz = np.array([[0,X[2],-X[1]], [-X[2],0,X[0]], [X[1],-X[0],0]])
    290286        if 'U' in TLStype:
    291287            Uout.append(['I',TLS[0],0,0,0,0,0,0])
    … …  
    302298                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
    303299            Umat = G2lat.U6toUij(U)
    304 #            print 'Umat: ',Umat
    305 #wrong?      Umat1 = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz.T,Lmat),Axyz)
    306 #            print 'Umat1: ',Umat1
    307300            beta = np.inner(np.inner(Bmat,Umat),Bmat.T)
    308301            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
     302        else:
     303            Uout.append(['N',])
    309304    return Uout
    310305   
    … …  
    14501445    '''
    14511446    A = 2.*acosd(Q[0])
    1452     V = np.array([0.,0.,1.])
     1447    V = np.array(Q[1:])
    14531448    if nl.norm(Q[1:]):
    14541449        V = Q[1:]/nl.norm(Q[1:])
    … …  
    14601455    '''
    14611456    A = 2.*np.arccos(Q[0])
    1462     V = np.array([0.,0.,1.])
     1457    V = np.array(Q[1:])
    14631458    if nl.norm(Q[1:]):
    14641459        V = Q[1:]/nl.norm(Q[1:])

  • trunk/GSASIIphsGUI.py

    r882 r885  
    30903090            RBObj = Indx[Obj.GetId()]
    30913091            val = Obj.GetValue()
    3092             RBObj['ThermalMotion'] = ['None',[],[]]
    30933092            if val == 'Uiso':
    3094                 RBObj['ThermalMotion'] = ['Uiso',[0.01,],[False,]]
     3093                RBObj['ThermalMotion'][0] = 'Uiso'
    30953094            elif val == 'T':
    3096                 RBObj['ThermalMotion'] = ['T',[0.0 for i in range(6)],[False for i in range(6)]]
     3095                RBObj['ThermalMotion'][0] = 'T'
    30973096            elif val == 'TL':
    3098                 RBObj['ThermalMotion'] = ['TL',[0.0 for i in range(12)],[False for i in range(12)]]
     3097                RBObj['ThermalMotion'][0] = 'TL'
    30993098            elif val == 'TLS':
    3100                 RBObj['ThermalMotion'] = ['TLS',[0.0 for i in range(20)],[False for i in range(20)]] #SAA = S33-S11 & SBB = S22-S33
     3099                RBObj['ThermalMotion'][0] = 'TLS'
    31013100            wx.CallAfter(FillRigidBodyGrid,True)
    31023101           
    … …  
    34503449                    return
    34513450                rbObj['Ids'] = Ids
    3452                 rbObj['ThermalMotion'] = ['None',[],[]] #type,values,flags
     3451                rbObj['ThermalMotion'] = ['None',[0. for i in range(21)],[False for i in range(21)]] #type,values,flags
    34533452                rbObj['RBname'] += ':'+str(RBData[rbType][rbId]['useCount'])
    34543453                RBObjs.append(rbObj)
    … …  
    37663765                QuatC = G2mth.prodQQ(QuatB,QuatA)
    37673766                rbObj['Orient'] = [QuatC,' ']
    3768                 rbObj['ThermalMotion'] = ['None',[],[]] #type,values,flags
     3767                rbObj['ThermalMotion'] = ['None',[0. for i in range(21)],[False for i in range(21)]] #type,values,flags
    37693768                SXYZ = []
    37703769                TXYZ = []

  • trunk/GSASIIplot.py

    r867 r885  
    34593459                    RenderLines(x,y,z,mapBonds[ind],Wt)
    34603460        if len(testRBObj) and pageName == 'RB Models':
    3461             XYZ = G2mth.UpdateRBAtoms(Bmat,testRBObj['rbObj'],testRBObj['rbData'],testRBObj['rbType'])
     3461            XYZ = G2mth.UpdateRBXYZ(Bmat,testRBObj['rbObj'],testRBObj['rbData'],testRBObj['rbType'])[0]
    34623462            rbBonds = FindPeaksBonds(XYZ)
    34633463            for ind,[x,y,z] in enumerate(XYZ):

  • trunk/GSASIIstruct.py

    r883 r885  
    3737   
    3838ateln2 = 8.0*math.log(2.0)
    39 DEBUG = False
     39DEBUG = False       #only for powders!
    4040
    4141def GetControls(GPXfile):
    … …  
    580580                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
    581581                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
    582                     pid = '::RBV;'+str(irb)+':'+str(i)
     582                    pid = '::RBV;'+str(i)+':'+str(irb)
    583583                    rbDict[pid] = mag
    584584                    if ref:
    … …  
    621621            RBmags = rigidbodyDict['Vector'][item]['VectMag']
    622622            for i,mag in enumerate(RBmags):
    623                 name = '::RBV;'+str(irb)+':'+str(i)
     623                name = '::RBV;'+str(i)+':'+str(irb)
    624624                mag = parmDict[name]
    625625                if name in sigDict:
    … …  
    645645        VRBData = RBData['Vector']
    646646        for i,rbId in enumerate(VRBIds):
    647             pfxRB = '::RBV;'+str(i)+':'
    648647            for j in range(len(VRBData[rbId]['VectMag'])):
    649                 VRBData[rbId]['VectMag'][j] = parmDict[pfxRB+str(j)]
     648                name = '::RBV;'+str(j)+':'+str(i)
     649                VRBData[rbId]['VectMag'][j] = parmDict[name]
    650650    for phase in Phases:
    651651        Phase = Phases[phase]
    … …  
    669669            if 'T' in TLS[0]:
    670670                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
    671                     RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
     671                    TLS[1][i] = parmDict[pfx+pt+rbsx]
    672672            if 'L' in TLS[0]:
    673673                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
    674                     RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
     674                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
    675675            if 'S' in TLS[0]:
    676676                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
    677                     RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
     677                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
    678678            if 'U' in TLS[0]:
    679                 RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBVU:'+rbsx]
     679                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
    680680            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
    681681            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
    … …  
    684684                for j in [0,1,2]:
    685685                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
    686                 if UIJ[0] == 'A':
     686                if UIJ[i][0] == 'A':
    687687                    for j in range(6):
    688                         parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = Uij[j+2]
    689                 elif UIJ[0] == 'I':
    690                     parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = Uij[j+1]
     688                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
     689                elif UIJ[i][0] == 'I':
     690                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
     691           
    691692        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
    692693            jrb = RRBIds.index(RBObj['RBId'])
    … …  
    716717                for j in [0,1,2]:
    717718                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
    718                 if UIJ[0] == 'A':
     719                if UIJ[i][0] == 'A':
    719720                    for j in range(6):
    720                         parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = Uij[j+2]
    721                 elif UIJ[0] == 'I':
    722                     parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = Uij[j+1]
     721                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
     722                elif UIJ[i][0] == 'I':
     723                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
     724                   
     725def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
     726    atxIds = ['dAx:','dAy:','dAz:']
     727    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
     728    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
     729    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
     730    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
     731    PIds = ['Px:','Py:','Pz:']
     732    OIds = ['Oa:','Oi:','Oj:','Ok:']
     733    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
     734    if not RBIds['Vector'] and not RBIds['Residue']:
     735        return
     736    VRBIds = RBIds['Vector']
     737    RRBIds = RBIds['Residue']
     738    RBData = rigidbodyDict
     739    for item in parmDict:
     740        if 'RB' in item:
     741            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
     742    General = Phase['General']
     743    cell = General['Cell'][1:7]
     744    Amat,Bmat = G2lat.cell2AB(cell)
     745    rpd = np.pi/180.
     746    rpd2 = rpd**2
     747    g = np.inner(Bmat,Bmat.T)
     748    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
     749        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
     750    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
     751    pfx = str(Phase['pId'])+'::'
     752    RBModels =  Phase['RBModels']
     753    dx = 0.0001
     754    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
     755        VModel = RBData['Vector'][RBObj['RBId']]
     756        Q = RBObj['Orient'][0]
     757        Pos = RBObj['Orig'][0]
     758        jrb = VRBIds.index(RBObj['RBId'])
     759        rbsx = str(irb)+':'+str(jrb)
     760        dXdv = []
     761        for iv in range(len(VModel['VectMag'])):
     762            dXdv.append(np.inner(Bmat,VModel['rbVect'][iv]).T)
     763        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
     764        for ia,atId in enumerate(RBObj['Ids']):
     765            atNum = AtLookup[atId]
     766            for iv in range(len(VModel['VectMag'])):
     767                for ix in [0,1,2]:
     768                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
     769            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
     770                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
     771            for iv in range(4):
     772                Q[iv] -= dx
     773                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
     774                Q[iv] += 2.*dx
     775                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
     776                Q[iv] -= dx
     777                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
     778                for ix in [0,1,2]:
     779                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
     780            X = G2mth.prodQVQ(Q,Cart[ia])
     781            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T*gvec
     782            dFdu = G2lat.U6toUij(dFdu.T)
     783            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
     784            dFdu = G2lat.UijtoU6(dFdu)
     785            atNum = AtLookup[atId]
     786            if 'T' in RBObj['ThermalMotion'][0]:
     787                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
     788                    dFdvDict[pfx+name+rbsx] += dFdu[i]
     789            if 'L' in RBObj['ThermalMotion'][0]:
     790                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
     791                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
     792                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
     793                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
     794                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
     795                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
     796                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
     797                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[0]*X[1]*X[2]+
     798                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
     799            if 'S' in RBObj['ThermalMotion'][0]:
     800                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
     801                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
     802                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
     803                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
     804                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
     805                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
     806                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
     807                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
     808            if 'U' in RBObj['ThermalMotion'][0]:
     809                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
     810
     811    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
     812        Q = RBObj['Orient'][0]
     813        Pos = RBObj['Orig'][0]
     814        jrb = RRBIds.index(RBObj['RBId'])
     815        rbsx = str(irb)+':'+str(jrb)
     816        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     817        for ia,atId in enumerate(RBObj['Ids']):
     818            atNum = AtLookup[atId]
     819            dx = 0.0001
     820            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
     821                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
     822            for iv in range(4):
     823                Q[iv] -= dx
     824                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     825                Q[iv] += 2.*dx
     826                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     827                Q[iv] -= dx
     828                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
     829                for ix in [0,1,2]:
     830                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
     831            X = G2mth.prodQVQ(Q,Cart[ia])
     832            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T*gvec
     833            dFdu = G2lat.U6toUij(dFdu.T)
     834            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
     835            dFdu = G2lat.UijtoU6(dFdu)
     836            atNum = AtLookup[atId]
     837            if 'T' in RBObj['ThermalMotion'][0]:
     838                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
     839                    dFdvDict[pfx+name+rbsx] += dFdu[i]
     840            if 'L' in RBObj['ThermalMotion'][0]:
     841                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
     842                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
     843                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
     844                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
     845                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
     846                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
     847                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
     848                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[0]*X[1]*X[2]+
     849                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
     850            if 'S' in RBObj['ThermalMotion'][0]:
     851                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
     852                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
     853                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
     854                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
     855                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
     856                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
     857                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
     858                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
     859            if 'U' in RBObj['ThermalMotion'][0]:
     860                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
     861   
    723862       
    724863################################################################################
    … …  
    13251464        print >>pFile,sigstr
    13261465       
    1327     def PrintRBObjTLSAndSig(rbfx,rbsx):
     1466    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
    13281467        namstr = '  names :'
    13291468        valstr = '  values:'
    13301469        sigstr = '  esds  :'
    1331         for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
    1332             name = pfx+rbfx+pt+rbsx
    1333             namstr += '%12s'%(pt[:3])
    1334             valstr += '%12.5f'%(parmDict[name])
     1470        if 'T' in TLS:
     1471            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
     1472                name = pfx+rbfx+pt+rbsx
     1473                namstr += '%12s'%(pt[:3])
     1474                valstr += '%12.5f'%(parmDict[name])
     1475                if name in sigDict:
     1476                    sigstr += '%12.5f'%(sigDict[name])
     1477                else:
     1478                    sigstr += 12*' '
     1479            print >>pFile,namstr
     1480            print >>pFile,valstr
     1481            print >>pFile,sigstr
     1482        if 'L' in TLS:
     1483            namstr = '  names :'
     1484            valstr = '  values:'
     1485            sigstr = '  esds  :'
     1486            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
     1487                name = pfx+rbfx+pt+rbsx
     1488                namstr += '%12s'%(pt[:3])
     1489                valstr += '%12.3f'%(parmDict[name])
     1490                if name in sigDict:
     1491                    sigstr += '%12.3f'%(sigDict[name])
     1492                else:
     1493                    sigstr += 12*' '
     1494            print >>pFile,namstr
     1495            print >>pFile,valstr
     1496            print >>pFile,sigstr
     1497        if 'S' in TLS:
     1498            namstr = '  names :'
     1499            valstr = '  values:'
     1500            sigstr = '  esds  :'
     1501            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
     1502                name = pfx+rbfx+pt+rbsx
     1503                namstr += '%12s'%(pt[:3])
     1504                valstr += '%12.3f'%(parmDict[name])
     1505                if name in sigDict:
     1506                    sigstr += '%12.3f'%(sigDict[name])
     1507                else:
     1508                    sigstr += 12*' '
     1509            print >>pFile,namstr
     1510            print >>pFile,valstr
     1511            print >>pFile,sigstr
     1512        if 'U' in TLS:
     1513            name = pfx+rbfx+'U:'+rbsx
     1514            namstr = '  names :'+'%12s'%('U')
     1515            valstr = '  values:'+'%12.3f'%(parmDict[name])
    13351516            if name in sigDict:
    1336                 sigstr += '%12.5f'%(sigDict[name])
     1517                sigstr = '  esds  :'+'%12.3f'%(sigDict[name])
    13371518            else:
    1338                 sigstr += 12*' '
    1339         print >>pFile,namstr
    1340         print >>pFile,valstr
    1341         print >>pFile,sigstr
    1342         namstr = '  names :'
    1343         valstr = '  values:'
    1344         sigstr = '  esds  :'
    1345         for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
    1346             name = pfx+rbfx+pt+rbsx
    1347             namstr += '%12s'%(pt[:3])
    1348             valstr += '%12.3f'%(parmDict[name])
    1349             if name in sigDict:
    1350                 sigstr += '%12.3f'%(sigDict[name])
    1351             else:
    1352                 sigstr += 12*' '
    1353         print >>pFile,namstr
    1354         print >>pFile,valstr
    1355         print >>pFile,sigstr
    1356         namstr = '  names :'
    1357         valstr = '  values:'
    1358         sigstr = '  esds  :'
    1359         for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
    1360             name = pfx+rbfx+pt+rbsx
    1361             namstr += '%12s'%(pt[:3])
    1362             valstr += '%12.3f'%(parmDict[name])
    1363             if name in sigDict:
    1364                 sigstr += '%12.3f'%(sigDict[name])
    1365             else:
    1366                 sigstr += 12*' '
    1367         print >>pFile,namstr
    1368         print >>pFile,valstr
    1369         print >>pFile,sigstr
     1519                sigstr = '  esds  :'+12*' '
     1520            print >>pFile,namstr
     1521            print >>pFile,valstr
     1522            print >>pFile,sigstr
     1523           
    13701524       
    13711525    def PrintRBObjTorAndSig(rbsx):
    … …  
    14741628                print >>pFile,' Vector rigid body '
    14751629                PrintRBObjPOAndSig('RBV',rbsx)
    1476                 PrintRBObjTLSAndSig('RBV',rbsx)
     1630                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
    14771631            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
    1478                 jrb = VRBIds.index(RBObj['RBId'])
     1632                jrb = RRBIds.index(RBObj['RBId'])
    14791633                rbsx = str(irb)+':'+str(jrb)
    14801634                print >>pFile,' Residue rigid body '
    14811635                PrintRBObjPOAndSig('RBR',rbsx)
    1482                 PrintRBObjTLSAndSig('RBR',rbsx)
     1636                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
    14831637                PrintRBObjTorAndSig(rbsx)
    14841638            atomsSig = {}
    … …  
    28713025        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
    28723026        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
    2873 # or else do RB modification of dFdvDict here? Perhaps not...
    28743027    return dFdvDict
    28753028   
    … …  
    35293682    return yc,yb
    35303683   
    3531 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     3684def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
    35323685   
    35333686    def cellVaryDerv(pfx,SGData,dpdA):
    … …  
    36083761        if not Phase['General'].get('doPawley'):
    36093762            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
     3763            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    36103764        for iref,refl in enumerate(refList):
    36113765            if 'C' in calcControls[hfx+'histType']:        #CW powder
    … …  
    37473901                print 'TOF Undefined at present'
    37483902                raise Exception    #no TOF yet
    3749             #do atom derivatives -  for F,X & U so far             
     3903            #do atom derivatives -  for RB,F,X & U so far             
    37503904            corr = dervDict['int']/refl[9]
    37513905            if Ka2:
    … …  
    37543908                try:
    37553909                    aname = name.split(pfx)[1][:2]
    3756                     if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
     3910                    if aname not in ['Af','dA','AU','RB']: continue # skip anything not an atom or rigid body param
    37573911                except IndexError:
    37583912                    continue
    … …  
    37693923    return dMdv
    37703924
    3771 def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
    3772     parmdict.update(zip(varylist,values))
    3773     G2mv.Dict2Map(parmdict,varylist)
     3925def dervRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
     3926    parmDict.update(zip(varylist,values))
     3927    G2mv.Dict2Map(parmDict,varylist)
    37743928    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
    37753929    nvar = len(varylist)
    … …  
    37883942            xB = np.searchsorted(x,Limits[0])
    37893943            xF = np.searchsorted(x,Limits[1])
    3790             dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
    3791                 varylist,Histogram,Phases,calcControls,pawleyLookup)
     3944            dMdvh = np.sqrt(W[xB:xF])*getPowderProfileDerv(parmDict,x[xB:xF],
     3945                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
    37923946        elif 'HKLF' in histogram[:4]:
    37933947            Histogram = Histograms[histogram]
    … …  
    38013955            phfx = '%d:%d:'%(Phase['pId'],hId)
    38023956            SGData = Phase['General']['SGData']
    3803             A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
     3957            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    38043958            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    38053959            refList = Histogram['Data']
    3806             dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
    3807 # do RB modification of dFdvDict here; varylist will contain RB variables so dFdvDict needs corresponding entries
     3960            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
     3961            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    38083962            dMdvh = np.zeros((len(varylist),len(refList)))
    38093963            for iref,ref in enumerate(refList):
    38103964                if ref[6] > 0:
    3811                     dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
     3965                    dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
    38123966                    if calcControls['F**2']:
    38133967                        if ref[5]/ref[6] >= calcControls['minF/sig']:
    … …  
    38393993            dMdv = dMdvh
    38403994           
    3841     pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     3995    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
    38423996    if np.any(pVals):
    3843         dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
     3997        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
    38443998        dMdv = np.concatenate((dMdv.T,dpdv.T)).T
    38453999       
    38464000    return dMdv
    38474001
    3848 def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
    3849     parmdict.update(zip(varylist,values))
    3850     G2mv.Dict2Map(parmdict,varylist)
     4002def HessRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):
     4003    parmDict.update(zip(varylist,values))
     4004    G2mv.Dict2Map(parmDict,varylist)
    38514005    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
     4006    ApplyRBModels(parmDict,Phases,rigidbodyDict)        #,Update=True??
    38524007    nvar = len(varylist)
    38534008    Hess = np.empty(0)
    … …  
    38664021            xB = np.searchsorted(x,Limits[0])
    38674022            xF = np.searchsorted(x,Limits[1])
    3868             dMdvh = getPowderProfileDerv(parmdict,x[xB:xF],
    3869                 varylist,Histogram,Phases,calcControls,pawleyLookup)
     4023            dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
     4024                varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup)
    38704025            Wt = np.sqrt(W[xB:xF])[np.newaxis,:]
    38714026            Dy = dy[xB:xF][np.newaxis,:]
    … …  
    38924047            phfx = '%d:%d:'%(Phase['pId'],hId)
    38934048            SGData = Phase['General']['SGData']
    3894             A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
     4049            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    38954050            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    38964051            refList = Histogram['Data']
    3897             dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmdict)
    3898 # do RB modification of dFdvDict here; varylist will contain RB variables so dFdvDict needs corresponding entries
     4052            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
     4053            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    38994054            dMdvh = np.zeros((len(varylist),len(refList)))
    39004055            wdf = np.zeros(len(refList))
    39014056            for iref,ref in enumerate(refList):
    39024057                if ref[6] > 0:
    3903                     dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
     4058                    dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
    39044059                    if calcControls['F**2']:
    39054060                        if ref[5]/ref[6] >= calcControls['minF/sig']:
    … …  
    39374092        else:
    39384093            continue        #skip non-histogram entries
    3939     pNames,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     4094    pNames,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
    39404095    if np.any(pVals):
    3941         dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmdict,varylist)
     4096        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
    39424097        Vec += np.sum(dpdv*pWt*pVals,axis=1)
    39434098        Hess += np.inner(dpdv*pWt,dpdv)
    39444099    return Vec,Hess
    39454100
    3946 def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
    3947     parmdict.update(zip(varylist,values))
    3948     Values2Dict(parmdict, varylist, values)
    3949     G2mv.Dict2Map(parmdict,varylist)
     4101def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg):       
     4102    parmDict.update(zip(varylist,values))
     4103    Values2Dict(parmDict, varylist, values)
     4104    G2mv.Dict2Map(parmDict,varylist)
    39504105    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
    39514106    M = np.empty(0)
    39524107    SumwYo = 0
    39534108    Nobs = 0
    3954     ApplyRBModels(parmdict,Phases,rigidbodyDict)
     4109    ApplyRBModels(parmDict,Phases,rigidbodyDict)
    39554110    histoList = Histograms.keys()
    39564111    histoList.sort()
    … …  
    39734128            Histogram['sumwYo'] = np.sum(W[xB:xF]*y[xB:xF]**2)
    39744129            SumwYo += Histogram['sumwYo']
    3975             yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
     4130            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
    39764131                varylist,Histogram,Phases,calcControls,pawleyLookup)
    39774132            yc[xB:xF] += yb[xB:xF]
    … …  
    39944149            phfx = '%d:%d:'%(Phase['pId'],hId)
    39954150            SGData = Phase['General']['SGData']
    3996             A = [parmdict[pfx+'A%d'%(i)] for i in range(6)]
     4151            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    39974152            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    39984153            refList = Histogram['Data']
    3999             refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmdict)
     4154            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
    40004155            df = np.zeros(len(refList))
    40014156            sumwYo = 0
    … …  
    40074162            for i,ref in enumerate(refList):
    40084163                if ref[6] > 0:
    4009                     SCExtinction(ref,phfx,hfx,pfx,calcControls,parmdict,varylist) #puts correction in refl[13]
    4010                     ref[7] = parmdict[phfx+'Scale']*ref[9]
     4164                    SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13]
     4165                    ref[7] = parmDict[phfx+'Scale']*ref[9]
    40114166                    ref[7] *= ref[13]
    4012                     ref[8] = ref[5]/parmdict[phfx+'Scale']
     4167                    ref[8] = ref[5]/parmDict[phfx+'Scale']
    40134168                    if calcControls['F**2']:
    40144169                        if ref[5]/ref[6] >= calcControls['minF/sig']:
    … …  
    40524207        GoOn = dlg.Update(Rw,newmsg='%s%8.3f%s'%('All data Rw =',Rw,'%'))[0]
    40534208        if not GoOn:
    4054             parmdict['saved values'] = values
     4209            parmDict['saved values'] = values
    40554210            dlg.Destroy()
    40564211            raise Exception         #Abort!!
    4057     pDict,pVals,pWt = penaltyFxn(HistoPhases,parmdict,varylist)
     4212    pDict,pVals,pWt = penaltyFxn(HistoPhases,parmDict,varylist)
    40584213    if np.any(pVals):
    40594214        pSum = np.sum(pWt*pVals**2)
    … …  
    42224377        cPickle.dump(Histogram,fl,1)
    42234378        cPickle.dump(Phases,fl,1)
     4379        cPickle.dump(rigidbodyDict,fl,1)
    42244380        cPickle.dump(calcControls,fl,1)
    42254381        cPickle.dump(pawleyLookup,fl,1)

  • trunk/testGSASIIstruct.py

    r866 r885  
    2020    global Phases
    2121    Phases = cPickle.load(file)
     22    global RBData
     23    RBData = cPickle.load(file)
    2224    global calcControls
    2325    calcControls = cPickle.load(file)
    … …  
    3840    fplot = plotter.add('function test').gca()
    3941    yc,yb = G2st.getPowderProfile(parmDict,xdata[xB:xF],varylist,Histogram,Phases,calcControls,pawleyLookup)
    40     fplot.plot(xdata[xB:xF],yc+yb,'r')
     42    fplot.plot(xdata[xB:xF],yc+yb,'r',label='calc+bkg')
     43    fplot.legend()
    4144
    4245def test2(name,delt):
    … …  
    4952    xF = np.searchsorted(xdata,limits[1])
    5053    hplot = plotter.add('derivatives test for '+name).gca()
    51     ya = G2st.getPowderProfileDerv(parmDict,xdata[xB:xF],varyList,Histogram,Phases,calcControls,pawleyLookup)[0]
    52     hplot.plot(xdata[xB:xF],ya,'b')
     54    ya = G2st.getPowderProfileDerv(parmDict,xdata[xB:xF],varyList,Histogram,Phases,RBData,calcControls,pawleyLookup)[0]
     55    hplot.plot(xdata[xB:xF],ya,'b',label='analytic deriv')
    5356    if 'dA' in name:
    5457        name = ''.join(name.split('d'))
    … …  
    6164    y1 += yb
    6265    yn = (y1-y0)/(2.*delt)
    63     hplot.plot(xdata[xB:xF],yn,'r+')
    64     hplot.plot(xdata[xB:xF],ya-yn,'g')
     66    hplot.plot(xdata[xB:xF],yn,'r+',label='numeric deriv')
     67    hplot.plot(xdata[xB:xF],ya-yn,'g',label='diff')
     68    hplot.legend()
    6569   
    6670if __name__ == '__main__':
    … …  
    7579        print name,parmDict[name]
    7680    names = [
    77         ['0:0:Size;i',0.01],
    78         ['0:0:Mustrain:0',0.001],
    79         ['0:0:Mustrain:1',0.001],
     81        ['0::AUiso:0',0.001],
     82        ['::RBV;0:0',0.001],
     83        ['0::RBVT11:0:0',0.1],
     84        ['0::RBVL11:0:0',0.1],
     85        ['0::RBVPz:0:0',0.0001],
     86        ['0::RBVOa:0:0',0.001],
    8087        ]
    8188    for [name,delt] in names:
Note: See TracChangeset for help on using the changeset viewer.