Skip to content

Commit 3bca4db

Browse files
vondreelevondreele
committed
repair to Babinet fxns - now ok
defo derivs
1 parent d7cb7c3 commit 3bca4db

1 file changed

Lines changed: 22 additions & 26 deletions

File tree

‎GSASII/GSASIIstrMath.py‎

Lines changed: 22 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -605,14 +605,14 @@ def MakePolar(Orient,QB):
605605
dFFdSR[name] = SHR*ffkp #ok
606606
dFFdSR["Akappa'1:%d"%iAt] += -2.0*SQkp*SHR*orbs[term]*dffdkp/kappap #ok
607607
dFFdSI[name] = SHI*ffkp #ok
608-
dFFdSI["Akappa'1:%d"%iAt] += 2.0*SQkp*SHI*orbs[term]*dffdkp/kappap #ok
608+
dFFdSI["Akappa'1:%d"%iAt] += -2.0*SQkp*SHI*orbs[term]*dffdkp/kappap #ok
609609
#end test
610610
FFSH += SH*orbs[term]*ffkp
611611
dFFdS[name] = SH*ffkp #ok
612612
dFFdS["Akappa'1:%d"%iAt] += -2.0*SQkp*SH*orbs[term]*dffdkp/kappap #ok
613613
FF[:,iAt] = FFcore+FFval+FFSH
614614
FFR[:,iAt] = FFcore+FFval+FFSHR
615-
FFI[:,iAt] = FFSHI
615+
FFI[:,iAt] = np.round(FFSHI,10)
616616
else:
617617
FFR[:,iAt] = FF[:,iAt]
618618
atFlg.append(0.)
@@ -1232,12 +1232,8 @@ def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
12321232
Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),nOps*len(TwinLaw))
12331233
fp = np.reshape(((FFR+FP).T-Bab).T,cosp.shape)*Tcorr
12341234
fpp = np.reshape(Flack*(FFI+FPP),sinp.shape)*Tcorr
1235-
if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1236-
fa = np.array([fp*cosp,-fpp*sinp])
1237-
fb = np.array([fp*sinp,fpp*cosp])
1238-
else:
1239-
fa = np.array([fp*cosp,-fpp*sinp])
1240-
fb = np.array([fp*sinp,fpp*cosp])
1235+
fa = np.array([fp*cosp,-fpp*sinp])
1236+
fb = np.array([fp*sinp,fpp*cosp])
12411237
fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
12421238
fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag
12431239
if 'P' in hType: #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
@@ -1308,7 +1304,7 @@ def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
13081304
bij = Mast*Uij.T
13091305
dFdvDict = {}
13101306
dFdfr = np.zeros((nRef,mSize))
1311-
dFdff = np.zeros((nRef,nOps,mSize))
1307+
dFdff = np.zeros((2,nRef,nOps,mSize))
13121308
dFdx = np.zeros((nRef,mSize,3))
13131309
dFdui = np.zeros((nRef,mSize))
13141310
dFdua = np.zeros((nRef,mSize,6))
@@ -1360,33 +1356,32 @@ def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
13601356
if len(dffdSHR):
13611357
for item in dFFdSR:
13621358
dffdSHR[item] = np.hstack((dffdSHR[item],dFFdSR[item]))
1363-
dffdSHI[item] = np.hstack((dffdSHI[item],dFFdSI[item]))
13641359
else:
13651360
dffdSHR.update(dFFdSR)
1361+
if len(dffdSHI):
1362+
for item in dFFdSI:
1363+
dffdSHI[item] = np.hstack((dffdSHI[item],dFFdSI[item]))
1364+
else:
13661365
dffdSHI.update(dFFdSI)
13671366
Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),nOps)
1368-
dBabdA = np.repeat(np.exp(-parmDict[phfx+'BabU']*SQfactor),nOps)
1367+
dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
13691368
fotr = np.reshape(((FFR+FP).T-Bab).T,cosp.shape)*Tcorr
13701369
foti = np.reshape(Flack*(FFI+FPP),sinp.shape)*Tcorr
1371-
if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
1372-
fa = np.array([fotr*cosp,-foti*sinp])
1373-
fb = np.array([fotr*sinp,foti*cosp])
1374-
else:
1375-
fa = np.array([fotr*cosp,-foti*sinp])
1376-
fb = np.array([fotr*sinp,foti*cosp])
1370+
fa = np.array([fotr*cosp,-foti*sinp])
1371+
fb = np.array([fotr*sinp,foti*cosp])
13771372
fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,refBlk,nTwins)
13781373
fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl
13791374
fax = np.array([-fotr*sinp,-foti*cosp]) #positions array(2,refBlk,nEqv,nAtoms)
13801375
fbx = np.array([fotr*cosp,-foti*sinp])
13811376
#sum below is over Uniq
13821377
dfadfr = np.sum(fa/occ,axis=-2) #array(2,refBlk,nAtom) Fdata != 0 avoids /0. problem
1383-
dfadff = fa[0]*Tcorr*atFlg/fotr #ignores resonant scattering? no sum on Uniq; array(refBlk,nEqv,nAtom)
1378+
dfadff = np.array([cosp*Tcorr*atFlg,-Flack*sinp*Tcorr*atFlg]) # no sum on Uniq; array(2,refBlk,nEqv,nAtom)
13841379
dfadba = np.sum(-cosp*Tcorr,axis=-2) #array(refBlk,nAtom)
13851380
dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2)
13861381
dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nAtoms)
13871382
dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2)
13881383
dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem
1389-
dfbdff = fb[0]*Tcorr*atFlg/fotr #ignores resonant scattering? no sum on Uniq; array(refBlk,nEqv,nAtom)
1384+
dfbdff = np.array([sinp*Tcorr*atFlg,Flack*cosp*Tcorr*atFlg])
13901385
dfbdba = np.sum(-sinp*Tcorr,axis=-2)
13911386
dfadfl = np.sum(np.sum(-foti*sinp,axis=-1),axis=-1)
13921387
dfbdfl = np.sum(np.sum(foti*cosp,axis=-1),axis=-1)
@@ -1398,19 +1393,20 @@ def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
13981393
SB = fbs[0]+fbs[1]
13991394
if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
14001395
dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/nOps
1401-
dFdff[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadff[nxs,:,:,:]+fbs[:,:,nxs,nxs]*dfbdff[nxs,:,:,:],axis=0) #not summed on Uniq yet
1396+
dFdff[:,iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadff+fbs[:,:,nxs,nxs]*dfbdff,axis=0) #not summed on Uniq yet
14021397
dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
14031398
dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
14041399
dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
14051400
else:
14061401
dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/nOps
1407-
dFdff[iBeg:iFin] = (2.*SA[:,nxs,nxs]*dfadff+2.*SB[:,nxs,nxs]*dfbdff) #not summed on Uniq yet array(Nref,nEqv,nAtom)
1402+
dFdff[:,iBeg:iFin] = [2.*(fas[0,:,nxs,nxs]*dfadff[0]+fbs[0,:,nxs,nxs]*dfbdff[0]),
1403+
2.*(fas[1,:,nxs,nxs]*dfadff[1]+fbs[1,:,nxs,nxs]*dfbdff[1])] #not summed on Uniq yet array(Nref,nEqv,nAtom)
14081404
dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])
14091405
dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1])
14101406
dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])
14111407
dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl #array(nRef,)
1412-
# 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)])+ \
1413-
# fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
1408+
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)])+ \
1409+
fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T
14141410
iBeg += blkSize
14151411
# print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
14161412
#loop over atoms - each dict entry is list of derivatives for all the reflections
@@ -1429,11 +1425,11 @@ def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
14291425
for item in dffdSHR:
14301426
if 'Sh' in item or 'O' in item:
14311427
if i == int(item.split(':')[1]):
1432-
dFdvDict[pfx+'RBS'+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSHR[item],(nRef,-1)),axis=1)
1428+
dFdvDict[pfx+'RBS'+item] = np.sum(dFdff[0,:,:,i]*np.reshape(dffdSHR[item],(nRef,-1)),axis=1)
14331429
else:
14341430
if i == int(item.split(':')[1]):
1435-
dFdvDict[pfx+item] = np.sum(dFdff[:,:,i]*np.reshape(dffdSHR[item],(nRef,-1)),axis=1)- \
1436-
np.sum(dFdff[:,:,i]*np.reshape(dffdSHI[item],(nRef,-1)),axis=1)
1431+
dFdvDict[pfx+item] = np.sum(dFdff[0,:,:,i]*np.reshape(dffdSHR[item],(nRef,-1)),axis=1)- \
1432+
np.sum(dFdff[1,:,:,i]*np.reshape(dffdSHI[item],(nRef,-1)),axis=1)
14371433
dFdvDict[phfx+'Flack'] = 4.*dFdfl.T
14381434
dFdvDict[phfx+'BabA'] = dFdbab.T[0]
14391435
dFdvDict[phfx+'BabU'] = dFdbab.T[1]

0 commit comments

Comments
 (0)