Changeset 1993


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

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

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1992 r1993  
    994994    Mast: array orthogonalization matrix for Uij
    995995    '''
    996    
    997     nxs = np.newaxis
    998     glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    999     Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    1000     dGdAx = np.zeros_like(Ax)
    1001     Bx = np.array(XSSdata[3:]).T   #...cos pos mods
    1002     dGdBx = np.zeros_like(Bx)
    1003     Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
    1004     dGdAf = np.zeros_like(Af)
    1005     Bf = np.array(FSSdata[1]).T    #cos frac mods...
    1006     dGdBf = np.zeros_like(Bf)
    1007     Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    1008     dGdAu = np.zeros_like(Au)
    1009     Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1010     dGdBu = np.zeros_like(Bu)
    1011    
    1012     if 'Fourier' in waveTypes:
    1013         nf = 0
    1014         nx = 0
    1015         XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    1016         FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    1017         if 'Crenel' in waveTypes:
    1018             nC = np.where('Crenel' in waveTypes)
    1019             nf = 1
    1020             #FmodZ = 0   replace
    1021     else:
    1022         nx = 1
    1023         if 'Sawtooth' in waveTypes:
    1024             nS = np.where('Sawtooth' in waveTypes)
    1025             #XmodZ = 0   replace
    1026     if Ax.shape[1]:
    1027         tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1028         StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #also dXmodA/dAx
    1029         CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #also dXmodB/dBx
    1030         XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32
    1031         XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
    1032         Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    1033         D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    1034         HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
    1035     if Af.shape[1]:
    1036         tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    1037         StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmodA/dAf
    1038         CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmodB/dBf
    1039         FmodA = Af[:,nf:,nxs]*StauF   #atoms X Fwaves X 32
    1040         FmodB = Bf[:,nf:,nxs]*CtauF
    1041         Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    1042     else:
    1043         Fmod = np.ones_like(HdotX)           
    1044     if Au.shape[1]:
    1045         tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    1046         StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodA/dAu
    1047         CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodB/dBu
    1048         UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32
    1049         UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
    1050         Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    1051         HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
    1052     else:
    1053         HbH = np.ones_like(HdotX)
    1054     HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]  #ops x atoms x waves x 32 x xyz
    1055     HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]
    1056     dHdXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:]    #ops x atoms x waves x 32 x xyz
    1057     dHdXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:]              #ditto
    1058     sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #ops x atoms x waves x xyz; sum for G-L integration
    1059     cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
    1060     dGdAx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXA*np.sin(twopi*HdotXA)+np.cos(twopi*HdotXA))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
    1061 # ops x atoms x waves x xyz
    1062     dGdBx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXB*np.cos(twopi*HdotXB)-np.sin(twopi*HdotXB))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
     996    import scipy.misc as sm
     997   
     998    cosHA,sinHA = Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)
     999    deltx = 0.00001
     1000    deltf = 0.0001
     1001    deltu = 0.000001
     1002    Mf = FSSdata.T.shape
     1003    dGdMfC = np.zeros(Mf)
     1004    dGdMfS = np.zeros(Mf)
     1005    Mx = XSSdata.T.shape
     1006    Mx = [H.shape[0],]+list(Mx)   #atoms x waves x sin+cos pos mods
     1007    dGdMxC = np.zeros(Mx)
     1008    dGdMxS = np.zeros(Mx)
     1009    for i in range(Mx[1]):  #atoms
     1010        for j in range(Mx[2]):  #waves
     1011            for k in range(Mx[3]):  #coeff
     1012                XSSdata[k,j,i] += deltx
     1013                Cap,Sap = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
     1014                XSSdata[k,j,i] -= 2*deltx
     1015                Cam,Sam = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast))
     1016                XSSdata[k,j,i] += deltx
     1017                dGdMxC[:,:,j,k] += (Cap-Cam)/(2*deltx)
     1018                dGdMxS[:,:,j,k] += (Sap-Sam)/(2*deltx)
     1019    Mu = USSdata.T.shape
     1020    dGdMuC = np.zeros(Mu)
     1021    dGdMuS = np.zeros(Mu)
    10631022#    GSASIIpath.IPyBreak()                 
    1064     return np.array([cosHA,sinHA]),np.concatenate((dGdAf,dGdBf),-1), \
    1065         np.concatenate((dGdAx,dGdBx),-1),np.concatenate((dGdAu,dGdBu),-1)      #ops X atoms
     1023               
     1024   
     1025#    nxs = np.newaxis
     1026#    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
     1027#    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
     1028#    dGdBx = np.zeros_like(Bx)
     1029#    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
     1030#    dGdAf = np.zeros_like(Af)
     1031#    Bf = np.array(FSSdata[1]).T    #cos frac mods...
     1032#    dGdBf = np.zeros_like(Bf)
     1033#    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
     1034#    dGdAu = np.zeros_like(Au)
     1035#    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
     1036#    dGdBu = np.zeros_like(Bu)
     1037#   
     1038#    GSASIIpath.IPyBreak()                 
     1039#    if 'Fourier' in waveTypes:
     1040#        nf = 0
     1041#        nx = 0
     1042#        XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     1043#        FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
     1044#        if 'Crenel' in waveTypes:
     1045#            nC = np.where('Crenel' in waveTypes)
     1046#            nf = 1
     1047#            #FmodZ = 0   replace
     1048#    else:
     1049#        nx = 1
     1050#        if 'Sawtooth' in waveTypes:
     1051#            nS = np.where('Sawtooth' in waveTypes)
     1052#            #XmodZ = 0   replace
     1053#    if Ax.shape[1]:
     1054#        tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     1055#        StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dAx
     1056#        CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #also dXmod/dBx
     1057#        XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32
     1058#        XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
     1059#        Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     1060#        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
     1061#        HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
     1062#    if Af.shape[1]:
     1063#        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
     1064#        StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
     1065#        CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
     1066#        FmodA = Af[:,nf:,nxs]*StauF   #atoms X Fwaves X 32
     1067#        FmodB = Bf[:,nf:,nxs]*CtauF
     1068#        Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
     1069#    else:
     1070#        Fmod = np.ones_like(HdotX)           
     1071#    if Au.shape[1]:
     1072#        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
     1073#        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dAu
     1074#        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmod/dBu
     1075#        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32
     1076#        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
     1077#        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
     1078#        HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
     1079#    else:
     1080#        HbH = np.ones_like(HdotX)
     1081#    HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:]  #ops x atoms x waves x 32 x xyz
     1082#    HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:]
     1083#    dHdXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:]    #ops x atoms x waves x 32 x xyz
     1084#    dHdXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:]              #ditto
     1085#    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #ops x atoms x waves x xyz; sum for G-L integration
     1086#    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
     1087#    dGdAx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXA*np.sin(twopi*HdotXA)+np.cos(twopi*HdotXA))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
     1088## ops x atoms x waves x xyz - imaginary part
     1089#    dGdBx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXB*np.cos(twopi*HdotXB)-np.sin(twopi*HdotXB))*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
     1090## ops x atoms x waves x xyz - real part
     1091    return np.array([cosHA,sinHA]),[dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS]
    10661092   
    10671093def posFourier(tau,psin,pcos,smul):
     
    22962322                    h,k,l,m = -hkl+Hmax
    22972323                    Fhkl[h,k,l,m] = dF*phasem
    2298     SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]                #4D map
     2324    SSrho = fft.fftn(fft.fftshift(Fhkl))/(10.*cell[6])          #4D map
    22992325    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
    23002326    map4DData['rho'] = np.real(SSrho)
  • trunk/GSASIIphsGUI.py

    r1988 r1993  
    13781378                    atomData[r][c] = atomData[r][c].replace(rbAtmDict.get(atomData[r][ci+8],''),'')
    13791379                if 'Atoms' in data['Drawing']:
     1380                    ci = colLabels.index('I/A')
    13801381                    DrawAtomsReplaceByID(data['Drawing'],ci+8,atomData[r],ID)
    13811382                wx.CallAfter(Paint)
     
    24832484                Map['MapType'] = mapType.GetValue()
    24842485               
    2485             Map['Resolution'] = 0.25
     2486            Map['Resolution'] = 0.5
    24862487            refsList = data['Histograms'].keys()
    24872488            mapSizer = wx.BoxSizer(wx.HORIZONTAL)
  • trunk/GSASIIplot.py

    r1991 r1993  
    44754475                G2frame.tau = 0.
    44764476            elif key in ['+','=']:
    4477                 G2frame.tau += 0.05
     4477                G2frame.tau += 0.1
    44784478            elif key == '-':
    4479                 G2frame.tau -= 0.05
     4479                G2frame.tau -= 0.1
    44804480            G2frame.tau %= 1.   #force 0-1 range
    44814481            G2frame.G2plotNB.status.SetStatusText('Modulation tau = %.2f'%(G2frame.tau),1)
     
    50565056# end of useful debug
    50575057        mapData = generalData['Map']
     5058        D4mapData = generalData.get('4DmapData',{})
    50585059        pageName = ''
    50595060        page = getSelection()
     
    50615062            pageName = G2frame.dataDisplay.GetPageText(page)
    50625063        rhoXYZ = []
    5063         if len(mapData['rho']):
     5064        if len(D4mapData.get('rho',[])):        #preferentially select 4D map if there
     5065            rho = D4mapData['rho'][:,:,:,int(G2frame.tau*10)]   #pick current tau 3D slice
     5066        elif len(mapData['rho']):               #ordinary 3D map
     5067            rho = mapData['rho']
     5068        if len(rho):
    50645069            VP = drawingData['viewPoint'][0]-np.array([.5,.5,.5])
    50655070            contLevel = drawingData['contourLevel']*mapData['rhoMax']
    50665071            if 'delt-F' in mapData['MapType'] or 'N' in mapData.get('Type',''):
    5067                 rho = ma.array(mapData['rho'],mask=(np.abs(mapData['rho'])<contLevel))
     5072                rho = ma.array(rho,mask=(np.abs(rho)<contLevel))
    50685073            else:
    5069                 rho = ma.array(mapData['rho'],mask=(mapData['rho']<contLevel))
     5074                rho = ma.array(rho,mask=(rho<contLevel))
    50705075            steps = 1./np.array(rho.shape)
    50715076            incre = np.where(VP>=0,VP%steps,VP%steps-steps)
     
    50775082            rcube = 2000.*Vol/(ForthirdPI*Nc)
    50785083            rmax = math.exp(math.log(rcube)/3.)**2
    5079             radius = min(drawingData['mapSize']**2,rmax)
     5084            radius = min(drawingData.get('mapSize',10.)**2,rmax)
    50805085            view = drawingData['viewPoint'][0]
    50815086            Rok = np.sum(np.inner(Amat,rhoXYZ-view).T**2,axis=1)>radius
  • 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.