Changeset 4054

Timestamp:

Jul 6, 2019 7:00:41 PM (4 years ago)

Author:

vondreele

Message:

changes to incommensurate magnetism math

Location:

trunk

Files:

• GSASIImath.py (modified) (2 diffs)
• GSASIIstrMath.py (modified) (2 diffs)

trunk/GSASIImath.py

@@ -1402 +1402 @@
     return Mmod,MmodA,MmodB    #Ntau,Nops,Natm,,Mxyz; sum,sin & cos parts
         
+def MagMod2(XYZ,modQ,MSSdata,SGData,SSGData):
+    '''
+    this needs to make magnetic moment modulations & magnitudes
+    '''
+    Am = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
+    Bm = np.array(MSSdata[3:]).T  #...cos pos mods
+    nWaves = Am.shape[1]
+    if not nWaves:
+        return 0.0,0.0
+    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
+    Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
+    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
+    if SGData['SGInv']:
+        SGMT = np.vstack((SGMT,-SGMT))
+        Sinv =np.vstack((Sinv,-Sinv))
+        SGT = np.vstack((SGT,-SGT))
+    SGMT = np.vstack([SGMT for cen in SGData['SGCen']])     #Nops,3,3
+    Sinv = np.vstack([Sinv for cen in SGData['SGCen']])     #Nops,4,4
+    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
+    detSM = nl.det(SGMT)            #Nops
+    epsinv = Sinv[:,3,3]            #Nops
+    kdr = np.inner(XYZ,modQ).T      #Nops,Natm
+    phase = kdr+(epsinv*(np.inner(modQ,SGT[:,:3])-SGT[:,3]))[:,nxs]     #Nops,Natm
+    
+    psin = np.sin(twopi*phase)      #Nops,Natm
+    pcos = np.cos(twopi*phase)
+    MmodB = np.sum(Bm[nxs,:,:,:]*pcos[:,:,nxs,nxs],axis=2)      #Nops,Natm,3
+    MmodA = np.sum(Am[nxs,:,:,:]*psin[:,:,nxs,nxs],axis=2)
+#    if SGData['SGGray']:
+#        MmodA = -np.sum(SGMT[:,nxs,:,:]*MmodA[:,:,nxs,:],axis=-1)*detSM[:,nxs,nxs]*SGData['SpnFlp'][:,nxs,nxs]
+#        MmodB = -np.sum(SGMT[:,nxs,:,:]*MmodB[:,:,nxs,:],axis=-1)*detSM[:,nxs,nxs]*SGData['SpnFlp'][:,nxs,nxs]
+#    else:
+    MmodA = np.sum(SGMT[:,nxs,:,:]*MmodA[:,:,nxs,:],axis=-1)*SGData['MagMom'][:,nxs,nxs]
+    MmodB = np.sum(SGMT[:,nxs,:,:]*MmodB[:,:,nxs,:],axis=-1)*SGData['MagMom'][:,nxs,nxs]
+    Mmod = MmodA+MmodB
+    return Mmod,MmodA,MmodB    #Nops,Natm,,Mxyz; sum,sin & cos parts
+        
 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
     '''
@@ -1431 +1468 @@
     return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
     
-def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
-    '''
-    H: array nRefBlk x ops X hklt
-    HP: array nRefBlk x ops X hklt proj to hkl
-    nWaves: list number of waves for frac, pos, uij & mag
-    Fmod: array 2 x atoms x waves    (sin,cos terms)
-    Xmod: array atoms X 3 X ngl
-    Umod: array atoms x 3x3 x ngl
-    glTau,glWt: arrays Gauss-Lorentzian pos & wts
-    '''
-    
-    if nWaves[2]:       #uij (adp) waves
-        if len(HP.shape) > 2:
-            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
-        else:
-            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
-    else:
-        HbH = 1.0
-    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
-    if len(H.shape) > 2:
-        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
-        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
-    else:
-        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
-        HdotXD = twopi*(HdotX+D[:,nxs,:])
-    M = np.swapaxes(Mmod,1,2)
-    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
-    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
-    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
-
+#def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
+#    '''
+#    H: array nRefBlk x ops X hklt
+#    HP: array nRefBlk x ops X hklt proj to hkl
+#    nWaves: list number of waves for frac, pos, uij & mag
+#    Fmod: array 2 x atoms x waves    (sin,cos terms)
+#    Xmod: array atoms X 3 X ngl
+#    Umod: array atoms x 3x3 x ngl
+#    glTau,glWt: arrays Gauss-Lorentzian pos & wts
+#    '''
+#    
+#    if nWaves[2]:       #uij (adp) waves
+#        if len(HP.shape) > 2:
+#            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
+#        else:
+#            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
+#    else:
+#        HbH = 1.0
+#    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
+#    if len(H.shape) > 2:
+#        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
+#        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
+#    else:
+#        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
+#        HdotXD = twopi*(HdotX+D[:,nxs,:])
+#    M = np.swapaxes(Mmod,1,2)
+#    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
+#    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
+#    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
+#
 def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
     '''

trunk/GSASIIstrMath.py

@@ -1510 +1510 @@
         
         mXYZ = np.array([[xyz[0] for xyz in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T])%1. #Natn,Nop,xyz
-        Tmag,TmagA,TmagB = G2mth.MagMod(ngl,mXYZ,modQ,MSSdata,SGData,SSGData)   #Ntau,Nops,Natm,Mxyz-Tmag matches drawing moments
+        Tmag,TmagA,TmagB = G2mth.MagMod2(mXYZ,modQ,MSSdata,SGData,SSGData)   #Nops,Natm,Mxyz-Tmag matches drawing moments
         
         if not SGData['SGGray']:    #for fixed Mx,My,Mz
@@ -1612 +1612 @@
             M = np.array(np.abs(H[3]),dtype=np.int)-1
         
-            fam = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*cosm[i,nxs,:,:,nxs]-    \
-                np.sign(H[3,i])*TmagA*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])
-            fbm = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*sinm[i,nxs,:,:,nxs]+    \
-                np.sign(H[3,i])*TmagA*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])
+            fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*cosm[i,:,:,nxs]-    \
+                np.sign(H[3,i])*TmagA*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
+            fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*sinm[i,:,:,nxs]+    \
+                np.sign(H[3,i])*TmagA*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
                        
             famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
             fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
             
-            fa = np.sum(famq,axis=-1)**2-np.sum(eM.T[:,nxs,:]*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
-            fb = np.sum(fbmq,axis=-1)**2-np.sum(eM.T[:,nxs,:]*fbmq,axis=-1)**2
-            
-            fas = np.sum(fa*glWt,axis=1)
-            fbs = np.sum(fb*glWt,axis=1)
+            fas = np.sum(famq,axis=-1)**2-np.sum(eM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
+            fbs = np.sum(fbmq,axis=-1)**2-np.sum(eM.T*fbmq,axis=-1)**2
             
             refl.T[10] = np.where(H[3],fas+fbs,fas0+fbs0)