Changeset 2089


Ignore:
Timestamp:
Dec 11, 2015 1:34:43 PM (6 years ago)
Author:
vondreele
Message:

fix mod vector derivs
work on ZigZag/Block? derivs - seem ok on value & powder data, not pos or single xtal.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r2084 r2089  
    10341034    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
    10351035    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
    1036     ZtauXt = np.zeros((Ax.shape[0],2,ngl))               #atoms x Tminmax x ngl
     1036    ZtauXt = np.zeros((Ax.shape[0],2,3,ngl))               #atoms x Tminmax x 3 x ngl
    10371037    ZtauXx = np.zeros((Ax.shape[0],3,ngl))               #atoms x XYZmax x ngl
    10381038    for iatm in range(Ax.shape[0]):
    … …  
    10421042            Tmm = Ax[iatm][0][:2]                       
    10431043            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
    1044             for i in range(2):
    1045                 Tmm[i] += dT
    1046                 Zp = posZigZag(glTau,Tmm,XYZmax).T
    1047                 Tmm[i] -= 2.*dT
    1048                 Zm = posZigZag(glTau,Tmm,XYZmax).T
    1049                 Tmm[i] += dT
    1050                 ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
    1051 #            for i in range(3):
    1052 #                XYZmax[i] += dX
    1053 #                Zp = posZigZag(glTau,Tmm,XYZmax).T
    1054 #                XYZmax[i] -= 2.*dX
    1055 #                Zm = posZigZag(glTau,Tmm,XYZmax).T
    1056 #                XYZmax[i] += dX
    1057 #                ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX)
    1058             dAdX = posZigZagDerv(glTau,Tmm,XYZmax)
    1059             ZtauXx[iatm] = dAdX
     1044            ZtauXt[iatm],ZtauXx[iatm] = posZigZagDerv(glTau,Tmm,XYZmax)
    10601045        elif 'Block' in waveTypes[iatm]:
    10611046            nx = 1
    10621047            Tmm = Ax[iatm][0][:2]                       
    10631048            XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]])           
    1064             for i in range(2):
    1065                 Tmm[i] += dT
    1066                 Zp = posBlock(glTau,Tmm,XYZmax).T
    1067                 Tmm[i] -= 2.*dT
    1068                 Zm = posBlock(glTau,Tmm,XYZmax).T
    1069                 Tmm[i] += dT
    1070                 ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT)
    1071 #            for i in range(3):
    1072 #                XYZmax[i] += dX
    1073 #                Zp = posBlock(glTau,Tmm,XYZmax).T
    1074 #                XYZmax[i] -= 2.*dX
    1075 #                Zm = posBlock(glTau,Tmm,XYZmax).T
    1076 #                XYZmax[i] += dX
    1077             ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax)
     1049            ZtauXt[iatm],ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax)
    10781050        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
    10791051        if nx:   
    … …  
    11451117    dHdXA = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
    11461118    dHdXB = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:]    #ditto - cos waves
     1119# ops x atoms x waves x 2xyz - real part - good
    11471120    dGdMxCa = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11481121    dGdMxCb = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11491122    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
    1150 # ops x atoms x waves x 2xyz - real part - good
     1123# ops x atoms x waves x 2xyz - imag part - good
    11511124    dGdMxSa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
    11521125    dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
    11531126    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
    1154 # ZigZag/Block waves
    1155     dHdXZt = twopi*np.ones(HP.shape[0])[:,nxs,nxs,nxs]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:]          #??ops x atoms x ngl x 3(ZigZag/Block Tminmax)
     1127# ZigZag/Block waves - problems here?
     1128    dHdXZt = -twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:,:]          #ops x atoms x ngl x 2(ZigZag/Block Tminmax)
    11561129    dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:]          #ops x atoms x ngl x 3(ZigZag/Block XYZmax)
    1157     dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1130    dGdMzCt = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11581131    dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
    1159     dGdMzC = np.concatenate((dGdMzCt,dGdMzCx),axis=-1)
    1160     dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
     1132    dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1)
     1133    dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    11611134    dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)
    1162     dGdMzS = np.concatenate((dGdMzSt,dGdMzSx),axis=-1)
     1135    dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1)
    11631136#    GSASIIpath.IPyBreak()
    1164 # ops x atoms x waves x 2xyz - imaginary part - good
    11651137    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS]
    11661138   
    … …  
    11701142    return np.sum(A,axis=0)+np.sum(B,axis=0)
    11711143   
    1172 def posZigZag(tau,Tmm,XYZmax):
     1144def posZigZag(T,Tmm,Xmax):
    11731145    DT = Tmm[1]-Tmm[0]
    1174     slopeUp = 2.*XYZmax/DT
    1175     slopeDn = 2.*XYZmax/(1.-DT)
    1176     A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in tau])
     1146    Su = 2.*Xmax/DT
     1147    Sd = 2.*Xmax/(1.-DT)
     1148    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax+Su*((t-Tmm[0])%1.),Xmax-Sd*((t-Tmm[1])%1.)) for t in T])
    11771149    return A
    11781150   
    1179 def posZigZagDerv(tau,Tmm,XYZmax):
     1151def posZigZagDerv(T,Tmm,Xmax):
    11801152    DT = Tmm[1]-Tmm[0]
    1181     slopeUp = 2.*XYZmax/DT
    1182     slopeDn = 2.*XYZmax/(1.-DT)
    1183     dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1+2*(t-Tmm[0])%1./DT,1-2*(t-Tmm[1])%1./DT) for t in tau])
    1184     return dAdX
    1185    
    1186 
    1187 def posBlock(tau,Tmm,XYZmax):
    1188     A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in tau])
     1153    Su = 2.*Xmax/DT
     1154    Sd = 2.*Xmax/(1.-DT)
     1155    dAdT = np.zeros((2,3,len(T)))
     1156    dAdT[0] = np.array([np.where(Tmm[0] < t <= Tmm[1],Su*(t-Tmm[0]-1)/DT,-Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
     1157    dAdT[1] = np.array([np.where(Tmm[0] < t <= Tmm[1],-Su*(t-Tmm[0])/DT,Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
     1158    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.+2.*(t-Tmm[0])/DT,1.-2.*(t-Tmm[1])%1./DT) for t in T])
     1159    return dAdT,dAdX
     1160
     1161def posBlock(T,Tmm,Xmax):
     1162    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T])
    11891163    return A
    11901164   
    1191 def posBlockDerv(tau,Tmm,XYZmax):
    1192     dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.,1.) for t in tau])
    1193     return dAdX
     1165def posBlockDerv(T,Tmm,Xmax):
     1166    dAdT = np.zeros((2,3,len(T)))
     1167    ind = np.searchsorted(T,Tmm)
     1168    dAdT[0,:,ind[0]] = -Xmax/len(T)
     1169    dAdT[1,:,ind[1]] = Xmax/len(T)
     1170    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T])  #OK
     1171    return dAdT,dAdX
    11941172   
    11951173def fracCrenel(tau,Toff,Twid):

  • trunk/GSASIIstrMath.py

    r2078 r2089  
    888888            dfbdfl = 0.0
    889889        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
     890        SA = fas[0]+fas[1]
     891        SB = fbs[0]+fbs[1]
    890892        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    891893            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
    … …  
    898900                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
    899901        else:
    900             SA = fas[0]+fas[1]
    901             SB = fbs[0]+fbs[1]
    902902            if nTwin > 1:
    903903                dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0][it]+SA[it]*dfadfr[1][it]+SB[it]*dfbdfr[0][it]+SB[it]*dfbdfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
    … …  
    12581258            dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)
    12591259            dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1)   
     1260#        GSASIIpath.IPyBreak()
    12601261        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
    12611262            dfadfl = np.sum(-FPP*Tcorr*sinp)
    … …  
    12681269        SB = fbs[0]+fbs[1]      #float = B+B' (might be array[nTwin])
    12691270        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    1270             dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)
    1271             dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1])    #array(nRef,nAtom,3)
    1272             dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1])   #array(nRef,nAtom)
    1273             dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1])    #array(nRef,nAtom,6)
    12741271            dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
    1275                            
    1276             dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1])      #array(nRef,natom,nwave,2)
    1277             dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1])      #array(nRef,natom,nwave,6)
    1278             dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
    1279             dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
    1280 #            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
    1281 #                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
    1282 #            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
    1283 #                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
    1284 #            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
    1285 #                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
    1286 #            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
    1287 #                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
    1288 #            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
    1289 #                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
    1290 #            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
    1291 #                2.*(fbs[0]*dfbdGx[0]+fbs[1]*dfbdGx[1])
    1292 #            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
    1293 #                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
     1272            dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+   \
     1273                2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     1274            dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+  \
     1275                2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     1276            dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+   \
     1277                2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     1278            dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+   \
     1279                2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     1280            dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+  \
     1281                2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1])
     1282            dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+  \
     1283                2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1])
     1284            dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+  \
     1285                2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1])
     1286            dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+  \
     1287                2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1])
    12941288        else:
    12951289            if nTwin > 1:
    … …  
    13151309                dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1])      #array(nRef,natom,5)
    13161310                dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1])      #array(nRef,natom,nwave,12)
     1311#            GSASIIpath.IPyBreak()
    13171312        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    13181313            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    … …  
    13641359            dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i]
    13651360           
    1366 #    GSASIIpath.IPyBreak()
     1361#        GSASIIpath.IPyBreak()
    13671362    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    13681363    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
    … …  
    19931988        vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    19941989        dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec)
     1990        h,k,l = [h+m*vec[0],k+m*vec[1],l+m*vec[2]]          #do proj of hklm to hkl so dPdA & dPdV come out right
    19951991    else:
    19961992        m = 0
    … …  
    20232019        dpdDB = 1./dsp
    20242020        dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5],
    2025             2*l*A[2]+h*A[4]+k*A[5]])*m**parmDict[hfx+'difC']*dsp**3/2.
     2021            2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2.
    20262022        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV
    20272023           
Note: See TracChangeset for help on using the changeset viewer.