Changeset 1990


Ignore:
Timestamp:
Oct 7, 2015 1:06:38 PM (7 years ago)
Author:
vondreele
Message:

Structure factor derivs for SS continued...

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1988 r1990  
    928928def Modulation(waveTypes,H,FSSdata,XSSdata,USSdata,Mast):
    929929    '''
    930     H: array ops X hklt
    931     Ax: array atoms X waves X xyz
    932     Bx: array atoms X waves X xyz
     930    H: array nRefBlk x ops X hklt
     931    FSSdata: array 2 x atoms x waves    (sin,cos terms)
     932    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
     933    USSdata: array 2x6 x atoms X waves (sin,cos terms)
     934    Mast: array orthogonalization matrix for Uij
    933935    '''
    934936   
    … …  
    987989    '''
    988990    H: array ops X hklt
    989     Ax: array atoms X waves X xyz
    990     Bx: array atoms X waves X xyz
     991    FSSdata: array 2 x atoms x waves    (sin,cos terms)
     992    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
     993    USSdata: array 2x6 x atoms X waves (sin,cos terms)
     994    Mast: array orthogonalization matrix for Uij
    991995    '''
    992996   
    … …  
    10201024            nS = np.where('Sawtooth' in waveTypes)
    10211025            #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
    10221035    if Af.shape[1]:
    10231036        tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    … …  
    10281041        Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    10291042    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
     1043        Fmod = np.ones_like(HdotX)           
    10381044    if Au.shape[1]:
    10391045        tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    … …  
    10451051        HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.?
    10461052    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
     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
    10511059    cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1)       #ditto
    1052 #   GSASIIpath.IPyBreak()                 
     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)
     1063#    GSASIIpath.IPyBreak()                 
    10531064    return np.array([cosHA,sinHA]),dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu      #ops X atoms
    10541065   

  • trunk/GSASIIstrMath.py

    r1988 r1990  
    12571257        Uniq = np.inner(H,SSGMT)
    12581258        Phi = np.inner(H,SSGT)
     1259        if SGInv:   #if centro - expand HKL sets
     1260            Uniq = np.vstack((Uniq,-Uniq))
     1261            Phi = np.hstack((Phi,-Phi))
    12591262        phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs])
    12601263        sinp = np.sin(phase)
    12611264        cosp = np.cos(phase)
    1262         occ = Mdata*Fdata/len(Uniq)
     1265        occ = Mdata*Fdata/Uniq.shape[0]
    12631266        biso = -SQfactor*Uisodata[:,nxs]
    12641267        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T    #ops x atoms
    … …  
    12681271        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
    12691272        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)     #ops x atoms
    1270         Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #ops x atoms
     1273        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0]  #ops x atoms
    12711274        fot = (FF+FP-Bab)*occ*Tcorr     #ops x atoms
    12721275        fotp = FPP*occ*Tcorr            #ops x atoms
    … …  
    12861289        #sum below is over Uniq
    12871290        dfadfr = np.sum(fag/occ,axis=1)        #Fdata != 0 ever avoids /0. problem
     1291        dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
    12881292        dfadba = np.sum(-cosp*(occ*Tcorr)[:,nxs],axis=1)
     1293        dfbdba = np.sum(-sinp*(occ*Tcorr)[:,nxs],axis=1)
    12891294        dfadui = np.sum(-SQfactor*fag,axis=1)
     1295        dfbdui = np.sum(-SQfactor*fbg,axis=1)
    12901296        if nTwin > 1:
    12911297            dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     1298            dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])           
    12921299            dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
     1300            dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)])
    12931301            # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
    12941302        else:
    12951303            dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1304            dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=-2)           
    12961305            dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,np.newaxis],axis=-2)
     1306            dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=-2)
    12971307            # array(2,nAtom,3) & array(2,nAtom,6)
    12981308        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!   
    1299         if not SGData['SGInv']:
    1300             dfbdfr = np.sum(fbg/occ,axis=1)        #Fdata != 0 avoids /0. problem
    1301             dfbdba = np.sum(-sinp*Tcorr[:,np.newaxis],axis=1)
    1302             dfbdui = np.sum(-SQfactor*fb,axis=1)
    1303             if len(TwinLaw) > 1:
    1304                 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])           
    1305                 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,np.newaxis],axis=2) for it in range(nTwin)])
    1306             else:
    1307                 dfadfl = np.sum(-FPP*Tcorr*sinp)
    1308                 dfbdfl = np.sum(FPP*Tcorr*cosp)
    1309                 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)           
    1310                 dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,np.newaxis],axis=2)
     1309        if not SGData['SGInv'] and len(TwinLaw) == 1:   #Flack derivative
     1310            dfadfl = np.sum(-FPP*Tcorr*sinp)
     1311            dfbdfl = np.sum(FPP*Tcorr*cosp)
    13111312        else:
    1312             dfbdfr = np.zeros_like(dfadfr)
    1313             dfbdx = np.zeros_like(dfadx)
    1314             dfbdui = np.zeros_like(dfadui)
    1315             dfbdua = np.zeros_like(dfadua)
    1316             dfbdba = np.zeros_like(dfadba)
    1317             dfadfl = 0.0
    1318             dfbdfl = 0.0
     1313            dfadfl = 1.0
     1314            dfbdfl = 1.0
    13191315        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!   
    13201316        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
Note: See TracChangeset for help on using the changeset viewer.