Changeset 3808


Ignore:
Timestamp:
Feb 3, 2019 7:35:59 AM (3 years ago)
Author:
vondreele
Message:

latest attempt at incommensurate mag str fctr & some cleanup of dead ends

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3777 r3808  
    13491349        MmodA = Am[:,:,:,nxs]*np.sin(twopi*tauM[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
    13501350        MmodB = Bm[:,:,:,nxs]*np.cos(twopi*tauM[nxs,:,nxs,:]) #ditto
    1351         Mmod = np.sum(MmodA+MmodB,axis=1)                #atoms X 3 X ngl; sum waves
    1352         Mmod = np.swapaxes(Mmod,1,2)
     1351        Mmod = np.array([MmodB,MmodA])
     1352        Mmod = np.sum(Mmod,axis=2)                #atoms X 3 X ngl; sum waves
     1353        Mmod = np.swapaxes(Mmod,1,-1)       #ReIm,mxyz,Ntau,Natm
    13531354    else:
    13541355        Mmod = 1.0
    13551356    return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt
    13561357       
    1357 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
     1358def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    13581359    '''
    13591360    H: array nRefBlk x ops X hklt
     
    13631364    Xmod: array atoms X 3 X ngl
    13641365    Umod: array atoms x 3x3 x ngl
    1365     Mmod: array atoms X 3 X ngl         #TODO: add mag moment modulation math
    13661366    glTau,glWt: arrays Gauss-Lorentzian pos & wts
    13671367    '''
     
    14131413    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    14141414   
    1415 def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast):
     1415def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast):
    14161416    '''
    14171417    Only for Fourier waves for fraction, position & adp (probably not used for magnetism)
     
    14291429    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    14301430    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1431     Am = np.array(MSSdata[:3]).T   #...cos pos mods x waves x atoms
    1432     Bm = np.array(MSSdata[3:]).T   #...cos pos mods
    1433     nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]]
     1431    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]]
    14341432    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
    14351433    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
     
    14701468        UmodA = 0.
    14711469        UmodB = 0.
    1472     if nWaves[3]:
    1473         tauM = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
    1474         StauM = np.ones_like(Am)[:,:,:,nxs]*np.sin(twopi*tauM)[nxs,:,nxs,:]   #also dMmodA/dAm
    1475         CtauM = np.ones_like(Bm)[:,:,:,nxs]*np.cos(twopi*tauM)[nxs,:,nxs,:]   #also dMmodB/dBm
    1476         MmodA = Am[:,:,:,nxs]*StauM #atoms x waves x 3x3 x ngl
    1477         MmodB = Bm[:,:,:,nxs]*CtauM #ditto
    1478     else:
    1479         StauM = 1.0
    1480         CtauM = 1.0
    1481         MmodA = 0.
    1482         MmodB = 0.
    1483     return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB,[StauM,CtauM],MmodA+MmodB
     1470    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
    14841471   
    14851472def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
  • trunk/GSASIIstrMath.py

    r3806 r3808  
    11081108        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    11091109        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    1110         fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     1110        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #Mxyz,Nref
    11111111        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    11121112        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     
    15081508        return
    15091509    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    1510     ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
     1510    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)      #NB: Mmod is ReIm,Mxyz,Ntau,Natm
    15111511    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    15121512
    15131513    if parmDict[pfx+'isMag']:       #TODO: fix the math
    1514         GSdata = Gdata[:,nxs,:]+Mmod.T                  #Mxyz,Ntau,nAtm
    1515         GSdata = np.inner(GSdata.T,np.swapaxes(SGMT,1,2)).T            #apply sym. ops.
     1514        GSdata = Gdata[nxs,nxs,:,:]+Mmod    #ReIm,ntau,Mxyz,Natm
     1515        GSdata = np.inner(np.swapaxes(GSdata,2,3),np.swapaxes(SGMT,1,2)).T  #apply sym. ops.--> Mxyz,Nops,Natm,Ntau,ReIm
    15161516        if SGData['SGInv'] and not SGData['SGFixed']:
    15171517            GSdata = np.hstack((GSdata,-GSdata))       #inversion if any
    15181518        GSdata = np.hstack([GSdata for icen in range(Ncen)])        #dup over cell centering
    1519         GSdata = SGData['MagMom'][nxs,:,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
    1520         SMag = np.sqrt(np.sum((np.inner(GSdata.T,Ginv)*GSdata.T),axis=-1)).T
     1519        GSdata = SGData['MagMom'][nxs,:,nxs,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
     1520        SMag = np.sqrt(np.sum((np.inner(GSdata.T,Ginv)*GSdata.T),axis=-1)).T    #Nops,Natm,Ntau,ReIm
    15211521        Kdata = np.inner(GSdata.T,uAmat).T     #Cartesian unit vectors
    1522         Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)    #normalization --> unit vectors
    1523         Kdata /= Kmean      #mxyz,nops,ntau,natm
     1522        Kmean = np.mean(np.sqrt(np.sum(np.sum(Kdata**2,axis=0),axis=-1)),axis=0)    #normalization --> unit vectors
     1523        Kdata /= Kmean.T      #mxyz,nops,natm,ntau,ReIm
    15241524
    15251525    FF = np.zeros(len(Tdata))
     
    15851585        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    15861586        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1]  #refBlk x ops x atoms
    1587         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x refBlk x sym X atoms
     1587        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
     1588        if 'T' in calcControls[hfx+'histType']:
     1589            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     1590            fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1591        else:
     1592            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     1593            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    15881594        fams = 0.
    15891595        fbms = 0.
     
    15921598            MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    15931599            TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
    1594             mphase = np.array([phase+twopi*np.inner(cen,HP.T)[:,nxs,nxs] for cen in SGData['SGCen']])
    1595             mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
    1596             sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    1597             cosm = np.cos(mphase)                               #ditto
     1600            D = twopi*H.T[:,3:]*glTau[nxs,:]
     1601            mphase = phase[:,:,nxs,:]+D[:,nxs,:,nxs]
     1602            mphase = np.array([mphase+twopi*np.inner(cen,HP.T)[:,nxs,nxs,nxs] for cen in SGData['SGCen']])
     1603            mphase = np.concatenate(mphase,axis=1)    #remove extra axis; Nref,Nop,Natm
     1604            sinm = np.swapaxes(np.sin(mphase),2,3)    #--> Nref,Nop,Natm,Ntau
     1605            cosm = np.swapaxes(np.cos(mphase),2,3)                               #ditto
    15981606           
    15991607            HM = np.inner(Bmat,HP.T)                             #put into cartesian space
    16001608            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    1601             eDotK = np.sum(HM[:,:,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:],axis=0)
    1602             Q = HM[:,:,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:]-Kdata[:,nxs,:,:,:] #Mxyz,Nref,Nop,Ntau,Natm
     1609            eDotK = np.sum(HM[:,:,nxs,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:,:],axis=0)
     1610            Q = HM[:,:,nxs,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:,:]-Kdata[:,nxs,:,:,:,:] #Mxyz,Nref,Nop,Natm,Ntau,ReIm
     1611
     1612#flawed but "close"           
     1613            fam = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*cosm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,:]).T   #Mxyz,Nref,Nop,Natm,Ntau,ReIm
     1614            fbm = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*sinm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,:]).T  #--> ReIm,Ntau,Natm,Nop,Nref,Mxyzditto
    16031615           
    1604             fam = Q*TMcorr[nxs,:,nxs,nxs,:]*cosm[nxs,:,:,nxs,:]*SMag[nxs,nxs,:,:,:]*glWt[nxs,nxs,nxs,:,nxs]    #ditto
    1605             fbm = Q*TMcorr[nxs,:,nxs,nxs,:]*sinm[nxs,:,:,nxs,:]*SMag[nxs,nxs,:,:,:]*glWt[nxs,nxs,nxs,:,nxs]    #ditto
     1616            famg = (fam[0]*fbm[0]-fam[1]*fbm[1]).T          #--> Mxyz,Nref,Nop,Natm,Ntau
     1617            fbmg = (fam[0]*fbm[1]+fam[1]*fbm[1]).T
    16061618           
    1607             fams = np.sum(np.sum(fam,axis=2),axis=-1)                          #xyz,Nref,ntau
    1608             fbms = np.sum(np.sum(fbm,axis=2),axis=-1)                          #ditto
     1619            fams = np.sum(np.sum(famg,axis=2),axis=2)      #xyz,Nref,ntau; sum ops & atoms
     1620            fbms = np.sum(np.sum(fbmg,axis=2),axis=2)      #ditto
    16091621           
    1610             fams = np.sum(np.sum(fams,axis=0)**2,axis=-1)
    1611             fbms = np.sum(np.sum(fbms,axis=0)**2,axis=-1)
     1622            fams = np.sum(np.sum(fams**2,axis=0)*glWt,axis=-1)
     1623            fbms = np.sum(np.sum(fbms**2,axis=0)*glWt,axis=-1)
    16121624
    1613         if 'T' in calcControls[hfx+'histType']:
    1614             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    1615             fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1616         else:
    1617             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    1618             fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
     1625#why doesn't this work?
     1626           
     1627#            fam = Q*TMcorr[nxs,:,nxs,nxs,:]*cosm[nxs,:,:,nxs,:]*SMag[nxs,nxs,:,:,:]    #RI,Mxyz,Nref,Nop,Ntau,Natm
     1628#            fbm = Q*TMcorr[nxs,:,nxs,nxs,:]*sinm[nxs,:,:,nxs,:]*SMag[nxs,nxs,:,:,:]    #ditto
     1629#           
     1630#            famg = fa[:,nxs,:,:,nxs,:]*fam[nxs,:,:,:,:]-fb[:,nxs,:,:,nxs,:]*fbm[nxs,:,:,:,:]
     1631#            fbmg = fb[:,nxs,:,:,nxs,:]*fam[nxs,:,:,:,:]+fa[:,nxs,:,:,nxs,:]*fbm[nxs,:,:,:,:]
     1632#           
     1633#            fams = np.sum(np.sum(famg,axis=3),axis=-1)     #RI,xyz,Nref,ntau
     1634#            fbms = np.sum(np.sum(fbmg,axis=3),axis=-1)     #ditto
     1635#           
     1636#            fams = np.sum(np.sum(np.sum(fams**2,axis=0),axis=0)*glWt[nxs,nxs,nxs,:],axis=-1)
     1637#            fbms = np.sum(np.sum(np.sum(fbms**2,axis=0),axis=0)*glWt[nxs,nxs,nxs,:],axis=-1)
     1638
     1639#or this
     1640
     1641#            fam = TMcorr[:,nxs,:]*cosm    #ditto
     1642#            fbm = TMcorr[:,nxs,:]*sinm    #ditto
     1643#
     1644#            cosQ = np.sum(np.cos(twopi*Q)*SMag[nxs,nxs,:,:,:]*glWt[nxs,nxs,nxs,:,nxs],axis=3)           
     1645#            sinQ = np.sum(np.sin(twopi*Q)*SMag[nxs,nxs,:,:,:]*glWt[nxs,nxs,nxs,:,nxs],axis=3)
     1646#           
     1647#            fagm = fam[nxs,:,:,:]*cosQ-fbm[nxs,:,:,:]*sinQ
     1648#            fbgm = fbm[nxs,:,:,:]*cosQ+fam[nxs,:,:,:]*sinQ
     1649#           
     1650#            fams = np.sum(np.sum(fagm,axis=-1),axis=-1)                          #xyz,Nref,ntau
     1651#            fbms = np.sum(np.sum(fbgm,axis=-1),axis=-1)                          #ditto
     1652#           
     1653#            fams = np.sum(fams**2,axis=0)
     1654#            fbms = np.sum(fbms**2,axis=0)
     1655
    16191656        fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    16201657        fbg = fb*GfpuA[0]+fa*GfpuA[1]
     
    18051842    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    18061843    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
    1807     waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5]
     1844    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
    18081845    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    18091846    FF = np.zeros(len(Tdata))
     
    18701907        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
    18711908        fotp = FPP*Tcorr            #ops x atoms
    1872         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms
     1909        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
    18731910        dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    18741911        # GfpuA is 2 x ops x atoms
     
    20492086    mSize = len(Mdata)  #no. atoms
    20502087    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    2051     ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
    2052     waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5]
     2088    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)     #NB: Mmod is ReIm,Mxyz,Ntau,Natm
     2089    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast)
    20532090    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    20542091    FF = np.zeros(len(Tdata))
     
    21222159        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
    21232160        fotp = FPP*Tcorr            #ops x atoms
    2124         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms
     2161        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
    21252162        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    21262163        # GfpuA is 2 x ops x atoms
Note: See TracChangeset for help on using the changeset viewer.