Changeset 2123
- Timestamp:
- Jan 13, 2016 1:19:43 PM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIstrMath.py
r2122 r2123 1018 1018 dFdbab = np.zeros((nRef,2)) 1019 1019 dFdfl = np.zeros((nRef)) 1020 dFdtw = np.zeros((nRef))1021 1020 Flack = 1.0 1022 1021 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1100 1099 SB = fbs[0]+fbs[1] 1101 1100 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 1102 dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+ \ 1103 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT) 1104 dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+ \ 1105 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1]) 1106 dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+ \ 1107 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1]) 1108 dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+ \ 1109 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1]) 1101 dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT) 1102 # dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+ \ 1103 # 2.*(fbs[0,:,nxs]*dfbdfr[0]+fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT) 1104 dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0) 1105 # dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+ \ 1106 # 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1]) 1107 dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0) 1108 # dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+ \ 1109 # 2.*(fbs[0,:,nxs]*dfbdui[0]+fbs[1,:,nxs]*dfbdui[1]) 1110 dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0) 1111 # dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+ \ 1112 # 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1]) 1110 1113 else: 1111 1114 dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT) … … 1191 1194 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1192 1195 # H = np.array(refl[:3]) 1193 # H = np. squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)1196 # H = np.inner(H.T,TwinLaw) #maybe array(3,nTwins) or (3) 1194 1197 # TwMask = np.any(H,axis=-1) 1195 1198 # if iref in TwDict: … … 1214 1217 # HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1215 1218 # Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1216 # Hij = np. squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))1219 # Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)) 1217 1220 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1218 1221 # Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ … … 1332 1335 dFdtw = np.zeros((nRef,nTwin)) 1333 1336 time0 = time.time() 1334 nref = len(refDict['RefList'])/100 1335 for iref,refl in enumerate(refDict['RefList']): 1337 #reflection processing begins here - big arrays! 1338 iBeg = 0 1339 blkSize = 16 #no. of reflections in a block - optimized for speed 1340 while iBeg < nRef: 1341 iFin = min(iBeg+blkSize,nRef) 1342 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1343 H = refl.T[:3] 1344 H = np.inner(H.T,TwinLaw) #array(3,nTwins) 1345 TwMask = np.any(H,axis=-1) 1346 for ir in range(blkSize): 1347 iref = ir+iBeg 1348 if iref in TwDict: 1349 for i in TwDict[iref]: 1350 for n in range(NTL): 1351 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1352 TwMask = np.any(H,axis=-1) 1353 SQ = 1./(2.*refl.T[4])**2 # or (sin(theta)/lambda)**2 1354 SQfactor = 8.0*SQ*np.pi**2 1336 1355 if 'T' in calcControls[hfx+'histType']: 1337 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1338 H = np.array(refl[:3]) 1339 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 1340 TwMask = np.any(H,axis=-1) 1341 if iref in TwDict: 1342 for i in TwDict[iref]: 1343 for n in range(NTL): 1344 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1345 TwMask = np.any(H,axis=-1) 1346 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 1347 SQfactor = 8.0*SQ*np.pi**2 1356 if 'P' in calcControls[hfx+'histType']: 1357 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14]) 1358 else: 1359 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 1360 FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0) 1361 FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0) 1348 1362 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1349 Bab = parmDict[phfx+'BabA']*dBabdA1363 Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin) 1350 1364 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1351 FF = refDict['FF']['FF'][iref].T[Tindx].T1352 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or(nTwin,nSGOp,3)1365 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0) 1366 Uniq = np.inner(H,SGMT) # (nTwin,nSGOp,3) 1353 1367 Phi = np.inner(H,SGT) 1354 1368 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T … … 1360 1374 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1361 1375 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1362 Hij = np. squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))1376 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6)) 1363 1377 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1364 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T* occ1365 fot = (FF+FP-Bab)*Tcorr1378 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT) 1379 fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr 1366 1380 fotp = FPP*Tcorr 1367 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr]) 1368 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr]) 1369 # GSASIIpath.IPyBreak() 1381 if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms 1382 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr]) 1383 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr]) 1384 else: 1385 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr]) 1386 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr]) 1370 1387 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 1371 1388 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl … … 1373 1390 fbs[0] *= 0. 1374 1391 fas[1] *= 0. 1375 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,n twi,nEqv,nAtoms)1392 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,nRef,ntwi,nEqv,nAtoms) 1376 1393 fbx = np.array([fot*cosp,-fotp*sinp]) 1377 1394 #sum below is over Uniq 1378 dfadfr = np.sum( fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem1395 dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0) #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem 1379 1396 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1380 dfadui = np.sum(-SQfactor*fa,axis=-2) 1381 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1382 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1383 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 1397 dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0) 1398 dfadx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fax[:,:,:,:,:,nxs],axis=-3),axis=0) # nRef x nTwin x nAtoms x xyz; sum on ops & A,A' 1399 dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0) 1384 1400 if not SGData['SGInv']: 1385 dfbdfr = np.sum( fb/occ,axis=-2) #Fdata != 0 avoids /0. problem1401 dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0) #Fdata != 0 avoids /0. problem 1386 1402 dfadba /= 2. 1387 1403 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1388 dfbdui = np.sum( -SQfactor*fb,axis=-2)1389 dfbdx = np. array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])1390 dfbdua = np. array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])1404 dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0) 1405 dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0) 1406 dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0) 1391 1407 else: 1392 1408 dfbdfr = np.zeros_like(dfadfr) … … 1397 1413 SA = fas[0]+fas[1] 1398 1414 SB = fbs[0]+fbs[1] 1399 dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)] 1400 dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)] 1401 dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)] 1402 dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)] 1415 # GSASIIpath.IPyBreak() 1416 dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT) 1417 dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx 1418 dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui 1419 dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua 1403 1420 if SGData['SGInv']: #centrosymmetric; B=0 1404 dFdtw[i ref] = np.sum(TwMask[nxs,:]*fas,axis=0)**21421 dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2 1405 1422 else: 1406 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1407 dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1408 fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1423 dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1424 # dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1425 # fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1426 iBeg += blkSize 1409 1427 # GSASIIpath.IPyBreak() 1410 1428 print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 1411 1429 #loop over atoms - each dict entry is list of derivatives for all the reflections 1412 for i in range(len(Mdata)): #these all OK ?1430 for i in range(len(Mdata)): #these all OK 1413 1431 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1414 1432 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) … … 3667 3685 if im: # split to nontwin/twin versions 3668 3686 if len(TwinLaw) > 1: 3669 dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3687 dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) #? 3670 3688 else: 3671 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3689 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) #OK 3672 3690 else: 3673 3691 if len(TwinLaw) > 1:
Note: See TracChangeset
for help on using the changeset viewer.