Changeset 3415

Timestamp:

May 31, 2018 3:47:44 PM (5 years ago)

Author:

vondreele

Message:

change testDeriv.py so it can be run in a bat file.
mag moment derivatives now close to correct & are workable

Location:

trunk

Files:

• GSASIIstrMath.py (modified) (8 diffs)
• testDeriv.py (modified) (3 diffs)

trunk/GSASIIstrMath.py

@@ -860 +860 @@
         else:
             fotp = FPP*Tcorr     
-#            GSASIIpath.IPyBreak()
         if 'T' in calcControls[hfx+'histType']:
             fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
@@ -910 +909 @@
         dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \
                             fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
-#        GSASIIpath.IPyBreak()
         iBeg += blkSize
 #    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
@@ -982 +980 @@
     Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
     Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
-    Kdata = np.inner(Gdata.T,uAmat).T*np.sqrt(nl.det(Ginv))/Mag     #Cartesian unit vectors
+    VGi = np.sqrt(nl.det(Ginv))
+    Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag     #Cartesian unit vectors
     Uij = np.array(G2lat.U6toUij(Uijdata))
     bij = Mast*Uij.T
@@ -1075 +1074 @@
     :returns: dict dFdvDict: dictionary of derivatives
     '''
+    
     g = nl.inv(G)
     ast = np.sqrt(np.diag(G))
@@ -1097 +1097 @@
     Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
     dMdm = np.inner(Gdata.T,Ginv).T/Mag
+    Gones = np.ones_like(Gdata)
     Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
+    Gones = np.inner(Gones.T,SGMT).T
     if SGData['SGInv'] and not SGData['SGFixed']:
         Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
+        Gones = np.hstack((Gones,-Gones))       #inversion if any
     Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
+    Gones = np.hstack([Gones for icen in range(Ncen)])        #dup over cell centering
     Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
+    Gones = SGData['MagMom'][nxs,:,nxs]*Gones   #flip vectors according to spin flip
     Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
     VGi = np.sqrt(nl.det(Ginv))
     Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag       #make unit vectors in Cartesian space
-    dkdG = (np.inner(np.ones(3),uAmat)*VGi)[:,nxs,nxs]/Mag[nxs,:,:]
+#    Gones = np.ones_like(Gdata)
+    dkdG = (np.inner(Gones.T,uAmat).T*VGi)/Mag
     dkdm = dkdG-Kdata*dMdm[:,nxs,:]/Mag[nxs,:,:]
     dFdMx = np.zeros((nRef,mSize,3))
@@ -1159 +1165 @@
         Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
         dqdk = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T     #Mxyz**2,Nref
-#        NQ = np.where(np.abs(Q)>0.,1./np.abs(Q),0.)     #this sort of works esp for 1 axis moments
         NQ = np.sqrt(np.sum(Q*Q,axis=0))
         NQ = np.where(NQ > 0.,1./NQ,0.)
         dqdm = dqdk[:,:,:,nxs,nxs]*dkdm[:,nxs,nxs,:,:]   #Mxyz**2,Nref,Nops,Natms
-        dmx = NQ*Q*dMdm[:,nxs,nxs,:]
+        dmx = Q*dMdm[:,nxs,nxs,:]/2.
         dmx = dmx[nxs,:,:,:,:]+dqdm*Mag[nxs,nxs,nxs,:,:]
         
@@ -1177 +1182 @@
         dfadmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*cosm[nxs,nxs,:,:,:],axis=3)
         dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
-        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)    #deriv OK? not U12 & U23 in sarc
+        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
         # imaginary part; array(3,refBlk,nAtom,3) & array(3,refBlk,nAtom,6)
         dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem 
@@ -1187 +1192 @@
         dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
         dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
-        dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)/len(SGT)          #problems
+        dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)          #problems
         dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
         dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok

trunk/testDeriv.py

@@ -22 +22 @@
 
 import sys
+import os
 import platform
 if '2' in platform.python_version_tuple()[0]:
@@ -69 +70 @@
         self.testDerivPanel = wx.ScrolledWindow(self)        
         self.plotNB = plot.PlotNotebook()
+        self.testFile = ''
+        arg = sys.argv
+        if len(arg) > 1 and arg[1]:
+            try:
+                self.testFile = os.path.splitext(arg[1])[0]+u'.testDeriv'
+            except:
+                self.testFile = os.path.splitext(arg[1])[0]+'.testDeriv'
+            self.TestRead()
+            self.UpdateControls(None)
         
     def __init__(self, parent):
@@ -94 +104 @@
             if dlg.ShowModal() == wx.ID_OK:
                 self.dirname = dlg.GetDirectory()
-                testFile = dlg.GetPath()
-                file = open(testFile,'rb')
-                if '2' in platform.python_version_tuple()[0]:
-                    self.values = cPickle.load(file)
-                    self.HistoPhases = cPickle.load(file)
-                    (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
-                    self.parmDict = cPickle.load(file)
-                    self.varylist = cPickle.load(file)
-                    self.calcControls = cPickle.load(file)
-                    self.pawleyLookup = cPickle.load(file)
-                else:
-                    self.values = cPickle.load(file,encoding='Latin-1')
-                    self.HistoPhases = cPickle.load(file,encoding='Latin-1')
-                    (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
-                    self.parmDict = cPickle.load(file,encoding='Latin-1')
-                    self.varylist = cPickle.load(file,encoding='Latin-1')
-                    self.calcControls = cPickle.load(file,encoding='Latin-1')
-                    self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
-                self.use = [False for i in range(len(self.varylist+self.depVarList))]
-                self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
-                file.close()
-                groups,parmlist = G2mv.GroupConstraints(self.constrDict)
-                G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
+                self.testFile = dlg.GetPath()
+                self.TestRead()
                 self.UpdateControls(event)
-                print(G2mv.VarRemapShow(self.varylist))
-                print('Dependent Vary List:',self.depVarList)
         finally:
             dlg.Destroy()
+            
+    def TestRead(self):
+        file = open(self.testFile,'rb')
+        if '2' in platform.python_version_tuple()[0]:
+            self.values = cPickle.load(file)
+            self.HistoPhases = cPickle.load(file)
+            (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file)
+            self.parmDict = cPickle.load(file)
+            self.varylist = cPickle.load(file)
+            self.calcControls = cPickle.load(file)
+            self.pawleyLookup = cPickle.load(file)
+        else:
+            self.values = cPickle.load(file,encoding='Latin-1')
+            self.HistoPhases = cPickle.load(file,encoding='Latin-1')
+            (self.constrDict,self.fixedList,self.depVarList) = cPickle.load(file,encoding='Latin-1')
+            self.parmDict = cPickle.load(file,encoding='Latin-1')
+            self.varylist = cPickle.load(file,encoding='Latin-1')
+            self.calcControls = cPickle.load(file,encoding='Latin-1')
+            self.pawleyLookup = cPickle.load(file,encoding='Latin-1')
+        self.use = [False for i in range(len(self.varylist+self.depVarList))]
+        self.delt = [max(abs(self.parmDict[name])*0.0001,1e-6) for name in self.varylist+self.depVarList]
+        file.close()
+        groups,parmlist = G2mv.GroupConstraints(self.constrDict)
+        G2mv.GenerateConstraints(groups,parmlist,self.varylist,self.constrDict,self.fixedList,self.parmDict)
+        print(G2mv.VarRemapShow(self.varylist))
+        print('Dependent Vary List:',self.depVarList)
             
     def UpdateControls(self,event):