Changeset 1993 for trunk/GSASIImath.py


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/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)
Note: See TracChangeset for help on using the changeset viewer.