Changeset 5515

Timestamp:

Mar 16, 2023 5:56:09 AM (2 years ago)

Author:

vondreele

Message:

fixes to spin RB derivatives

File:

• trunk/GSASIIstrMath.py (modified) (12 diffs)

trunk/GSASIIstrMath.py ¶

@@ -364 +364 @@
         
     dFFdS = {}
+    atFlg = []
     for iAt,Atype in enumerate(Tdata):
         if 'Q' in Atype:
+            atFlg.append(1.0)
             SHdat = SHCdict[iAt]
             symAxis = np.array(SHdat['symAxis'])
@@ -387 +389 @@
             FF[:,iAt] = 0.
             ishl = 0
+#            dBSdR = np.zeros(HKL.shape[0]*HKL.shape[1])
+            dSHdO = np.zeros(HKL.shape[0]*HKL.shape[1])
+            dSHdOi = np.zeros(HKL.shape[0]*HKL.shape[1])
+            dSHdOj = np.zeros(HKL.shape[0]*HKL.shape[1])
+            dSHdOk = np.zeros(HKL.shape[0]*HKL.shape[1])
             while True:
                 shl = '%d'%ishl
@@ -407 +414 @@
                 R0M = sp.spherical_jn(0,QR*(R-0.01))/(4.*np.pi)
                 dBSdR = Nat*SFF*(R0P-R0M)/0.02
-                dSHdO = np.zeros(HKL.shape[0]*HKL.shape[1])
-                dSHdOi = np.zeros(HKL.shape[0]*HKL.shape[1])
-                dSHdOj = np.zeros(HKL.shape[0]*HKL.shape[1])
-                dSHdOk = np.zeros(HKL.shape[0]*HKL.shape[1])
                 FF[:,iAt] += Nat*SFF*R0    #Bessel function; L=0 term
                 for item in Shell:
@@ -438 +441 @@
                         dFFdS[name] = Nat*SFF*BS*SH
                         #fill derivatives here wrt iAt,ishl,item or l,m
+                ishl += 1
                 dFFdS[Rname] = dBSdR
-                dFFdS[Oname] = dSHdO
-                dFFdS[Oiname] = dSHdOi
-                dFFdS[Ojname] = dSHdOj
-                dFFdS[Okname] = dSHdOk
-                ishl += 1
+            dFFdS[Oname] = dSHdO
+            dFFdS[Oiname] = dSHdOi
+            dFFdS[Ojname] = dSHdOj
+            dFFdS[Okname] = dSHdOk
+        else:
+            atFlg.append(0.)
     if ifDeriv:
-        return dFFdS    
+        return dFFdS,atFlg    
             
 def GetSHC(pfx,parmDict):
@@ -1065 +1070 @@
         return {}
     mSize = len(Mdata)
+    nOps = len(SGMT)
     FF = np.zeros(len(Tdata))
     if 'NC' in calcControls[hfx+'histType'] or 'NB' in calcControls[hfx+'histType']:
@@ -1078 +1084 @@
     dFdvDict = {}
     dFdfr = np.zeros((nRef,mSize))
-    dFdff = np.zeros((nRef,mSize))
+    dFdff = np.zeros((nRef,nOps,mSize))
     dFdx = np.zeros((nRef,mSize,3))
     dFdui = np.zeros((nRef,mSize))
@@ -1123 +1129 @@
         Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
         if len(SHCdict):
-            dffdsh = MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,BLtables,FF,SQ,True)
+            dffdsh,atFlg = MakeSpHarmFF(Uniq,Bmat,SHCdict,Tdata,hType,FFtables,BLtables,FF,SQ,True)
             if len(dffdSH):
                 for item in dffdSH:
@@ -1146 +1152 @@
         #sum below is over Uniq
         dfadfr = np.sum(fa/occ,axis=-2)        #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
-        dfadff = np.sum(fa[0]*Tcorr/fot,axis=-2)           #ignores resonant scattering? #TODO needs only for Q atoms!
-        for i,aname in enumerate(Tdata):
-            if 'Q' not in aname:
-                dfadff[:,i] = 0.0
+#        dfadff = np.sum(fa[0]*Tcorr*atFlg/fot,axis=-2)           #ignores resonant scattering? #TODO needs only for Q atoms!
+        dfadff = fa[0]*Tcorr*atFlg/fot           #ignores resonant scattering? no sum on Uniq
         dfadba = np.sum(-cosp*Tcorr,axis=-2)  #array(refBlk,nAtom)
         dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
@@ -1157 +1161 @@
         if not SGData['SGInv']:
             dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
-            dfbdff = np.sum(fb[0]*Tcorr/fot,axis=-2)           #ignores resonant scattering?
-            for i,aname in enumerate(Tdata):
-                if 'Q' not in aname:
-                    dfbdff[:,i] = 0.0
+#            dfbdff = np.sum(fb[0]*Tcorr*atFlg/fot,axis=-2)           #ignores resonant scattering?
+            dfbdff = fb[0]*Tcorr*atFlg/fot           #ignores resonant scattering? no sum on Uniq
             dfbdba = np.sum(-sinp*Tcorr,axis=-2)
             dfadfl = np.sum(np.sum(-fotp*sinp,axis=-1),axis=-1)
@@ -1181 +1183 @@
         if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
             dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
-            dFdff[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadff+fbs[:,:,nxs]*dfbdff,axis=0)/len(SGMT)
+#            dFdff[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadff+fbs[:,:,nxs]*dfbdff,axis=0)/len(SGMT)
+            dFdff[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadff[nxs,:,:,:]+fbs[:,:,nxs,nxs]*dfbdff[nxs,:,:,:],axis=0) #not summed on Uniq yet
             dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
             dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
@@ -1187 +1190 @@
         else:
             dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
-            dFdff[iBeg:iFin] = (2.*SA[:,nxs]*dfadff+2.*SB[:,nxs]*dfbdff)/len(SGMT)
+#            dFdff[iBeg:iFin] = (2.*SA[:,nxs]*dfadff+2.*SB[:,nxs]*dfbdff)/len(SGMT)
+            dFdff[iBeg:iFin] = (2.*SA[:,nxs,nxs]*dfadff+2.*SB[:,nxs,nxs]*dfbdff)      #not summed on Uniq yet
             dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
             dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
@@ -1210 +1214 @@
         dFdvDict[pfx+'AU23:'+str(i)] = dFdua.T[5][i]
         for item in dffdSH:
-            dFdvDict[pfx+'RBS'+item] = dFdff.T[i]*np.sum(np.reshape(dffdSH[item],(nRef,-1)),axis=-1)
+            i = int(item.split(':')[1])
+            dFdvDict[pfx+'RBS'+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSH[item],(nRef,-1)),axis=-1)
     dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
     dFdvDict[phfx+'BabA'] = dFdbab.T[0]