Ignore:
Timestamp:
Oct 8, 2015 4:09:47 PM (7 years ago)
Author:
vondreele
Message:

numerical derivatives for x modulations - work but slow
show density map slices as tau changed for SS structures

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1992 r1993  
    11341134        dFdfl = np.zeros((nRef,nTwin))
    11351135        dFdtw = np.zeros((nRef,nTwin))
    1136         dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1],2))
    1137         dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],6))
    1138         dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],12))
     1136        dFdGfA = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
     1137        dFdGfB = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1]))
     1138        dFdGxA = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
     1139        dFdGxB = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3))
     1140        dFdGuA = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
     1141        dFdGuB = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6))
    11391142    else:
    11401143        dFdfr = np.zeros((nRef,mSize))
     
    11901193        fotp = FPP*Tcorr            #ops x atoms
    11911194        GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast)
     1195        # derivs are: ops x atoms x waves x 1,3,or 6 parms as [real,imag] parts
    11921196        fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms)
    11931197        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     
    12081212        dfadui = np.sum(-SQfactor*fag,axis=1)
    12091213        dfbdui = np.sum(-SQfactor*fbg,axis=1)
     1214#        GSASIIpath.IPyBreak()                 
    12101215        if nTwin > 1:
    1211             dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    1212             dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])           
    1213             dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    1214             dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    1215             # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
     1216            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1217            dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])           
     1218            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1219            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
     1220            # array(nTwin,2,nAtom,3) & array(2,nTwin,nAtom,6)
     1221            dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1222            dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1)
     1223            # array (2,nAtom,wave,3)
    12161224        else:
    1217             dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
    1218             dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=-2)           
    1219             dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,np.newaxis],axis=-2)
    1220             dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1225            dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)
     1226            dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2)           
     1227            dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2)
     1228            dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2)
    12211229            # array(2,nAtom,3) & array(2,nAtom,6)
     1230            dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1231            dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1)
     1232            # array (2,nAtom,wave,3)
    12221233        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
    12231234#        GSASIIpath.IPyBreak()
     
    12551266                dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1])   #array(nRef,nAtom)
    12561267                dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1])    #array(nRef,nAtom,6)
    1257                 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl  #array(nRef,)
     1268                dFdfl[iref] = -SA*dfadfl-SB*dfbdfl                  #array(nRef,)
     1269                dFdGx[iref] = 2.*SA*(dfadGx[0]+dfbdGx[1])+2.*SB*(dfadGx[1]+dfbdGx[0])      #array(nRef,natom,nwave,6)
    12581270        dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    12591271            2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    12601272        #loop over atoms - each dict entry is list of derivatives for all the reflections
    1261     for i in range(len(Mdata)):     
     1273    for i in range(len(Mdata)):     #loop over atoms 
    12621274        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    12631275        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     
    12711283        dFdvDict[pfx+'AU13:'+str(i)] = .5*dFdua.T[4][i]
    12721284        dFdvDict[pfx+'AU23:'+str(i)] = .5*dFdua.T[5][i]
    1273         #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations...
     1285        for j in range(XSSdata.shape[1]):
     1286            dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i]
     1287            dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i]
     1288            dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j)] = dFdGx.T[2][j][i]
     1289            dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j)] = dFdGx.T[3][j][i]
     1290            dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i]
     1291            dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i]
    12741292    dFdvDict[phfx+'BabA'] = dFdbab.T[0]
    12751293    dFdvDict[phfx+'BabU'] = dFdbab.T[1]
Note: See TracChangeset for help on using the changeset viewer.