Changeset 3861

Timestamp:

Mar 25, 2019 1:32:58 PM (5 years ago)

Author:

vondreele

Message:

another try at incomm. mag. str. fctr.

Location:

trunk

Files:

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

trunk/GSASIImath.py

@@ -1356 +1356 @@
     return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt
 
-def MagMod(XYZ,modQ,MSSdata,SSGData):
+def MagMod(XYZ,modQ,MSSdata,SGData,SSGData):
+    
+    SGMT = np.array([ops[0] for ops in SGData['SGOps']])
+    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
+    if SGData['SGInv']:
+        SGMT = np.vstack((SGMT,-SGMT))
+        SGT = np.vstack((SGT,-SGT))
+    SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
+    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
+
+#works for DyMn6Ge6 but not MnWO4
     Bm = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
     Am = np.array(MSSdata[3:]).T   #...cos pos mods
@@ -1366 +1376 @@
     if nWaves:
         modind = np.arange(nWaves)+1.
-        phase = np.sum(twopi*XYZ[:,:,nxs,:]*modind[nxs,nxs,:,nxs]*modQ[nxs,nxs,nxs,:],axis=-1)
-        phase = np.swapaxes(phase,0,1)  #Nops,Natm,Nwave
+        phaseA = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ))).T         #= 2pimk.r Nops,Natm,Nwave
         if nCen > 0:
-            phshp = phase.shape
-            phase = np.reshape(phase,(nCen,mEqv,phshp[1],-1))
+            phshp = phaseA.shape
+            phaseA = np.reshape(phaseA,(nCen,mEqv,phshp[1],-1))
             for ic,cen in enumerate(SSGData['SSGCen']):
-                phase[ic] += twopi*cen[3]/2.
-            phase = np.reshape(phase,phshp)                
-        MmodA = Am[nxs,:,:,:]*np.cos(phase[:,:,:,nxs])-Bm[nxs,:,:,:]*np.sin(phase[:,:,:,nxs])
-        MmodB = Am[nxs,:,:,:]*np.sin(phase[:,:,:,nxs])+Bm[nxs,:,:,:]*np.cos(phase[:,:,:,nxs])
+                phaseA[ic] += twopi*cen[3]/2.
+            phaseA = np.reshape(phaseA,phshp)                
+        MmodAA = Am[nxs,:,:,:]*np.cos(phaseA[:,:,:,nxs])-Bm[nxs,:,:,:]*np.sin(phaseA[:,:,:,nxs])
+        MmodBA = Am[nxs,:,:,:]*np.sin(phaseA[:,:,:,nxs])+Bm[nxs,:,:,:]*np.cos(phaseA[:,:,:,nxs])    
+    
+#fails    
+    modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = 1.
+    poff = SGData['SpnFlp']*np.inner(modQp,SGT)
+    Bm = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
+    Am = np.array(MSSdata[3:]).T   #...cos pos mods
+    nWaves = Am.shape[1]
+    MmodA = 0; MmodB = 0
+    if nWaves:
+        modind = np.arange(nWaves)+1.
+        phase = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ)+poff)).T
+        psin = np.sin(phase)
+        pcos = np.cos(phase)
+        RAM = np.rollaxis(np.inner(Am,SGMT),2)
+        RBM = np.rollaxis(np.inner(Bm,SGMT),2)
+        MmodA = SGData['MagMom'][nxs,nxs,:,nxs]*(RAM*pcos[:,:,:,nxs]-RBM*psin[:,:,:,nxs])
+        MmodB = SGData['MagMom'][nxs,nxs,:,nxs]*(RAM*psin[:,:,:,nxs]+RBM*pcos[:,:,:,nxs])
+#        MmodA = RAM*pcos[:,:,:,nxs]-RBM*psin[:,:,:,nxs]
+#        MmodB = RAM*psin[:,:,:,nxs]+RBM*pcos[:,:,:,nxs]   
+
     return MmodA,MmodB    #cos & sin Nops,Natm,Nwaves,Mxyz
         

trunk/GSASIIstrMath.py

@@ -1518 +1518 @@
         
         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
-        MmodA,MmodB = G2mth.MagMod(mXYZ,modQ,MSSdata,SSGData)   #Re cos/Im sin,Nops,Natm,Nwaves,Mxyz
-        MmodA = np.inner(MmodA,uAmat.T)*(SGData['MagMom']/SGData['SpnFlp'])[:,nxs,nxs,nxs]   #make cartesian * apply det(ops)
-        MmodB = np.inner(MmodB,uAmat.T)*(SGData['MagMom']/SGData['SpnFlp'])[:,nxs,nxs,nxs]
+        MmodA,MmodB = G2mth.MagMod(mXYZ,modQ,MSSdata,SGData,SSGData)   #Re cos/Im sin,Nops,Natm,Nwaves,Mxyz
+        MmodA *= (SGData['MagMom']/SGData['SpnFlp'])[:,nxs,nxs,nxs]   #apply det(ops)
+        MmodB *= (SGData['MagMom']/SGData['SpnFlp'])[:,nxs,nxs,nxs]
+        MmodA = np.inner(MmodA,uAmat.T)   #make cartesian
+        MmodB = np.inner(MmodB,uAmat.T)
         
     FF = np.zeros(len(Tdata))
@@ -1590 +1592 @@
             phasem = twopi*np.inner(H.T[:,:3],mXYZ)
             phasem = np.swapaxes(phasem,1,2)
+            cosm = np.cos(phasem)
             sinm = np.sin(phasem)
-            cosm = np.cos(phasem)
             MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
             TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm