Changeset 1988 for trunk/GSASIImath.py


Ignore:
Timestamp:
Oct 6, 2015 2:53:02 PM (8 years ago)
Author:
vondreele
Message:

fix atom id problem in AddRestraints?
refactor Modulation & ModulationDerv? to remove expModInt
fix AtomRefine?
fix Uij derivative in SS problems

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1987 r1988  
    926926################################################################################
    927927   
    928 def Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast):
     928def Modulation(waveTypes,H,FSSdata,XSSdata,USSdata,Mast):
     929    '''
     930    H: array ops X hklt
     931    Ax: array atoms X waves X xyz
     932    Bx: array atoms X waves X xyz
     933    '''
    929934   
    930935    nxs = np.newaxis
    931936    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    932    
    933     def expModInt(H,Af,Bf,Ax,Bx,Au,Bu):
    934         '''
    935         H: array ops X hklt
    936         Ax: array atoms X waves X xyz
    937         Bx: array atoms X waves X xyz
    938         S: array ops
    939         '''
    940         if 'Fourier' in waveTypes:
    941             nf = 0
    942             nx = 0
    943             XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    944             FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    945             if 'Crenel' in waveTypes:
    946                 nC = np.where('Crenel' in waveTypes)
    947                 nf = 1
    948                 #FmodZ = 0   replace
    949         else:
    950             nx = 1
    951             if 'Sawtooth' in waveTypes:
    952                 nS = np.where('Sawtooth' in waveTypes)
    953                 #XmodZ = 0   replace
    954         if Af.shape[1]:
    955             tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    956             FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X 32
    957             FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:])
    958             Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    959         else:
    960             Fmod = 1.0           
    961         if Ax.shape[1]:
    962             tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    963             XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
    964             XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    965             Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    966         if Au.shape[1]:
    967             tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    968             UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32
    969             UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
    970             Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    971             HbH = np.exp(-np.sum(H[:,:,nxs,nxs,:3]*np.inner(H[:,:,:3],Umod),axis=-1)) # refBlk x ops x atoms x 32 add Overhauser corr.?
    972         else:
    973             HbH = 1.0
    974         D = H[:,:,3][:,:,nxs]*glTau[nxs,nxs,:]              #m*tau;refBlk x ops X 32
    975         HdotX = np.inner(H[:,:,:3],np.swapaxes(Xmod,1,2))+D[:,:,nxs,:]     #refBlk X ops x atoms X 32
    976         sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #refBlk X ops x atoms; sum for G-L integration
    977         cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
    978 #        GSASIIpath.IPyBreak()                 
    979         return cosHA,sinHA     #each refBlk x ops X atoms
    980 
    981937    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    982938    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
    … …  
    985941    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    986942    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    987     GpA = np.array(expModInt(Uniq,Af,Bf,Ax,Bx,Au,Bu))
    988     return GpA             # 2 x refBlk x SGops x atoms
    989    
    990 def ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):
     943   
     944    if 'Fourier' in waveTypes:
     945        nf = 0
     946        nx = 0
     947        XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     948        FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
     949        if 'Crenel' in waveTypes:
     950            nC = np.where('Crenel' in waveTypes)
     951            nf = 1
     952            #FmodZ = 0   replace
     953    else:
     954        nx = 1
     955        if 'Sawtooth' in waveTypes:
     956            nS = np.where('Sawtooth' in waveTypes)
     957            #XmodZ = 0   replace
     958    if Af.shape[1]:
     959        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
     960        FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X 32
     961        FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:])
     962        Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
     963    else:
     964        Fmod = 1.0           
     965    if Ax.shape[1]:
     966        tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     967        XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
     968        XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     969        Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     970    if Au.shape[1]:
     971        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
     972        UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32
     973        UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
     974        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
     975        HbH = np.exp(-np.sum(H[:,:,nxs,nxs,:3]*np.inner(H[:,:,:3],Umod),axis=-1)) # refBlk x ops x atoms x 32 add Overhauser corr.?
     976    else:
     977        HbH = 1.0
     978    D = H[:,:,3][:,:,nxs]*glTau[nxs,nxs,:]              #m*tau;refBlk x ops X 32
     979    HdotX = np.inner(H[:,:,:3],np.swapaxes(Xmod,1,2))+D[:,:,nxs,:]     #refBlk X ops x atoms X 32
     980    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #refBlk X ops x atoms; sum for G-L integration
     981    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
     982#        GSASIIpath.IPyBreak()                 
     983
     984    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
     985   
     986def ModulationDerv(waveTypes,H,FSSdata,XSSdata,USSdata,Mast):
     987    '''
     988    H: array ops X hklt
     989    Ax: array atoms X waves X xyz
     990    Bx: array atoms X waves X xyz
     991    '''
    991992   
    992993    nxs = np.newaxis
    993994    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    994    
    995     def expModInt(H,Af,Bf,Ax,Bx,Au,Bu):
    996         '''
    997         H: array ops X hklt
    998         Ax: array atoms X waves X xyz
    999         Bx: array atoms X waves X xyz
    1000         S: array ops
    1001         '''
    1002         if 'Fourier' in waveTypes:
    1003             nf = 0
    1004             nx = 0
    1005             XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    1006             FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    1007             if 'Crenel' in waveTypes:
    1008                 nC = np.where('Crenel' in waveTypes)
    1009                 nf = 1
    1010                 #FmodZ = 0   replace
    1011         else:
    1012             nx = 1
    1013             if 'Sawtooth' in waveTypes:
    1014                 nS = np.where('Sawtooth' in waveTypes)
    1015                 #XmodZ = 0   replace
    1016         if Af.shape[1]:
    1017             tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    1018             FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X 32
    1019             FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:])
    1020             Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    1021         else:
    1022             Fmod = 1.0           
    1023         if Ax.shape[1]:
    1024             tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1025             XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
    1026             XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    1027             Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    1028         if Au.shape[1]:
    1029             tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    1030             UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32
    1031             UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
    1032             Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    1033             HbH = np.exp(np.swapaxes(np.array([-np.inner(h,np.inner(Umod,h)) for h in H[:,:3]]),1,2)).T #Overhauser corr.?
    1034         else:
    1035             HbH = 1.0
    1036         D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    1037         HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:]     #atoms X 32 X ops
    1038         sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #atoms X ops; sum for G-L integration
    1039         cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #ditto
    1040 #        GSASIIpath.IPyBreak()                 
    1041         return cosHA.T,sinHA.T      #ops X atoms
    1042 
    1043995    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    1044996    dGdAx = np.zeros_like(Ax)
    … …  
    10531005    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    10541006    dGdBu = np.zeros_like(Bu)
    1055     GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu))
    1056     return GpA,dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu
     1007   
     1008    if 'Fourier' in waveTypes:
     1009        nf = 0
     1010        nx = 0
     1011        XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     1012        FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
     1013        if 'Crenel' in waveTypes:
     1014            nC = np.where('Crenel' in waveTypes)
     1015            nf = 1
     1016            #FmodZ = 0   replace
     1017    else:
     1018        nx = 1
     1019        if 'Sawtooth' in waveTypes:
     1020            nS = np.where('Sawtooth' in waveTypes)
     1021            #XmodZ = 0   replace
     1022    if Af.shape[1]:
     1023        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
     1024        StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmodA/dAf
     1025        CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmodB/dBf
     1026        FmodA = Af[:,nf:,nxs]*StauF   #atoms X Fwaves X 32
     1027        FmodB = Bf[:,nf:,nxs]*CtauF
     1028        Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
     1029    else:
     1030        Fmod = 1.0           
     1031    if Ax.shape[1]:
     1032        tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     1033        StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #also dXmodA/dAx
     1034        CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #also dXmodB/dBx
     1035        XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32
     1036        XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto
     1037        Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     1038    if Au.shape[1]:
     1039        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
     1040        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodA/dAu
     1041        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodB/dBu
     1042        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32
     1043        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
     1044        Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
     1045        HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
     1046    else:
     1047        HbH = 1.0
     1048    D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
     1049    HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:]     #ops x atoms X 32
     1050    sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1)       #ops x atoms; sum for G-L integration
     1051    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
     1052#   GSASIIpath.IPyBreak()                 
     1053    return np.array([cosHA,sinHA]),dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu      #ops X atoms
    10571054   
    10581055def posFourier(tau,psin,pcos,smul):
Note: See TracChangeset for help on using the changeset viewer.