Changeset 2115
- Timestamp:
- Jan 7, 2016 9:48:49 AM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r2110 r2115 1009 1009 else: 1010 1010 D = H[:,3:]*glTau[nxs,:] #m*e*tau; refBlk x ops X ngl 1011 HdotXD = twopi*(HdotX+D[:,nxs,:]) 1012 cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1) #real part; refBlk X ops x atoms; sum for G-L integration 1013 sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1) #imag part; ditto 1014 return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms 1015 1016 def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt): 1017 ''' 1018 H: array nRefBlk x tw x ops X hklt 1019 HP: array nRefBlk x tw x ops X hklt proj to hkl 1020 Fmod: array 2 x atoms x waves (sin,cos terms) 1021 Xmod: array atoms X ngl X 3 1022 Umod: array atoms x ngl x 3x3 1023 glTau,glWt: arrays Gauss-Lorentzian pos & wts 1024 ''' 1025 1026 if nWaves[2]: 1027 if len(HP.shape) > 3: #Blocks of reflections 1028 HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.? 1029 else: #single reflections 1030 HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.? 1031 else: 1032 HbH = 1.0 1033 HdotX = np.inner(HP,Xmod) #refBlk x tw x ops x atoms X ngl 1034 if len(H.shape) > 3: 1035 D = glTau*H[:,:,:,3:,nxs] #m*e*tau; refBlk x tw x ops X ngl 1036 HdotXD = twopi*(HdotX+D[:,:,:,nxs,:]) 1037 else: 1038 D = H*glTau[nxs,:] #m*e*tau; refBlk x ops X ngl 1011 1039 HdotXD = twopi*(HdotX+D[:,nxs,:]) 1012 1040 cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1) #real part; refBlk X ops x atoms; sum for G-L integration -
trunk/GSASIIstrMath.py
r2111 r2115 1425 1425 return dFdvDict 1426 1426 1427 def SStructureFactor(refDict, im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):1427 def SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1428 1428 ''' 1429 Compute super structure factors for all h,k,l,m for phase 1429 Compute super structure factors for all h,k,l,m for phase - no twins 1430 puts the result, F^2, in each ref[9] in refList 1431 works on blocks of 32 reflections for speed 1432 input: 1433 1434 :param dict refDict: where 1435 'RefList' list where each ref = h,k,l,m,it,d,... 1436 'FF' dict of form factors - filed in below 1437 :param np.array G: reciprocal metric tensor 1438 :param str pfx: phase id string 1439 :param dict SGData: space group info. dictionary output from SpcGroup 1440 :param dict calcControls: 1441 :param dict ParmDict: 1442 1443 ''' 1444 phfx = pfx.split(':')[0]+hfx 1445 ast = np.sqrt(np.diag(G)) 1446 Mast = twopisq*np.multiply.outer(ast,ast) 1447 SGInv = SGData['SGInv'] 1448 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1449 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1450 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1451 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1452 FFtables = calcControls['FFtables'] 1453 BLtables = calcControls['BLtables'] 1454 Flack = 1.0 1455 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 1456 Flack = 1.-2.*parmDict[phfx+'Flack'] 1457 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1458 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 1459 ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast) 1460 modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 1461 FF = np.zeros(len(Tdata)) 1462 if 'NC' in calcControls[hfx+'histType']: 1463 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1464 elif 'X' in calcControls[hfx+'histType']: 1465 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1466 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1467 Uij = np.array(G2lat.U6toUij(Uijdata)).T 1468 bij = Mast*Uij 1469 blkSize = 32 #no. of reflections in a block 1470 nRef = refDict['RefList'].shape[0] 1471 if not len(refDict['FF']): 1472 if 'N' in calcControls[hfx+'histType']: 1473 dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's 1474 refDict['FF']['El'] = dat.keys() 1475 refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values() 1476 else: 1477 dat = G2el.getFFvalues(FFtables,0.) 1478 refDict['FF']['El'] = dat.keys() 1479 refDict['FF']['FF'] = np.ones((nRef,len(dat))) 1480 for iref,ref in enumerate(refDict['RefList']): 1481 SQ = 1./(2.*ref[5])**2 1482 dat = G2el.getFFvalues(FFtables,SQ) 1483 refDict['FF']['FF'][iref] *= dat.values() 1484 time0 = time.time() 1485 #reflection processing begins here - big arrays! 1486 iBeg = 0 1487 while iBeg < nRef: 1488 iFin = min(iBeg+blkSize,nRef) 1489 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1490 H = refl.T[:4] #array(blkSize,4) 1491 HP = H[:3]+modQ[:,nxs]*H[3:] #projected hklm to hkl 1492 SQ = 1./(2.*refl.T[5])**2 #array(blkSize) 1493 SQfactor = 4.0*SQ*twopisq #ditto prev. 1494 Uniq = np.inner(H.T,SSGMT) 1495 UniqP = np.inner(HP.T,SGMT) 1496 Phi = np.inner(H.T,SSGT) 1497 if SGInv: #if centro - expand HKL sets 1498 Uniq = np.hstack((Uniq,-Uniq)) 1499 Phi = np.hstack((Phi,-Phi)) 1500 UniqP = np.hstack((UniqP,-UniqP)) 1501 if 'T' in calcControls[hfx+'histType']: 1502 if 'P' in calcControls[hfx+'histType']: 1503 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14]) 1504 else: 1505 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 1506 FP = np.repeat(FP.T,Uniq.shape[1],axis=0) 1507 FPP = np.repeat(FPP.T,Uniq.shape[1],axis=0) 1508 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]) 1509 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1510 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1],axis=0) 1511 phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs]) 1512 sinp = np.sin(phase) 1513 cosp = np.cos(phase) 1514 biso = -SQfactor*Uisodata[:,nxs] 1515 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T 1516 HbH = -np.sum(UniqP[:,:,nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1) #use hklt proj to hkl 1517 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1518 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1] #refBlk x ops x atoms 1519 if 'T' in calcControls[hfx+'histType']: 1520 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 1521 fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) 1522 else: 1523 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 1524 fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) 1525 GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms 1526 fag = fa*GfpuA[0]-fb*GfpuA[1] #real; 2 x refBlk x sym x atoms 1527 fbg = fb*GfpuA[0]+fa*GfpuA[1] 1528 fas = np.sum(np.sum(fag,axis=-1),axis=-1) #2 x refBlk; sum sym & atoms 1529 fbs = np.sum(np.sum(fbg,axis=-1),axis=-1) 1530 # GSASIIpath.IPyBreak() 1531 if 'P' in calcControls[hfx+'histType']: 1532 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 #square of sums 1533 # refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) 1534 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1535 else: 1536 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 #square of sums 1537 refl.T[8] = np.copy(refl.T[10]) 1538 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1539 iBeg += blkSize 1540 print 'nRef %d time %.4f\r'%(nRef,time.time()-time0) 1541 1542 def SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1543 ''' 1544 Compute super structure factors for all h,k,l,m for phase - twins only 1430 1545 puts the result, F^2, in each ref[8+im] in refList 1431 works on blocks of 100reflections for speed1546 works on blocks of 32 reflections for speed 1432 1547 input: 1433 1548 … … 1496 1611 iFin = min(iBeg+blkSize,nRef) 1497 1612 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1498 H = refl.T[:4] #array(blkSize,4) 1499 HP = H[:3]+modQ[:,nxs]*H[3:] #projected hklm to hkl 1500 # HP = np.squeeze(np.inner(HP.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4) 1501 # TwMask = np.any(HP,axis=-1) 1502 # if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i] 1503 # for ir in range(blkSize): 1504 # iref = ir+iBeg 1505 # if iref in TwDict: 1506 # for i in TwDict[iref]: 1507 # for n in range(NTL): 1508 # HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1509 # TwMask = np.any(HP,axis=-1) 1613 H = refl[:,:4] #array(blkSize,4) 1614 H3 = refl[:,:3] 1615 HP = H[:,:3]+modQ[nxs,:]*H[:,3:] #projected hklm to hkl 1616 HP = np.inner(HP,TwinLaw) #array(blkSize,nTwins,4) 1617 H3 = np.inner(H3,TwinLaw) 1618 TwMask = np.any(HP,axis=-1) 1619 if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i] 1620 for ir in range(blkSize): 1621 iref = ir+iBeg 1622 if iref in TwDict: 1623 for i in TwDict[iref]: 1624 for n in range(NTL): 1625 HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1626 H3[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1627 TwMask = np.any(HP,axis=-1) 1510 1628 SQ = 1./(2.*refl.T[5])**2 #array(blkSize) 1511 1629 SQfactor = 4.0*SQ*twopisq #ditto prev. 1512 Uniq = np.inner(H.T,SSGMT) 1513 UniqP = np.inner(HP.T,SGMT) 1514 Phi = np.inner(H.T,SSGT) 1630 Uniq = np.inner(H,SSGMT) 1631 Uniq3 = np.inner(H3,SGMT) 1632 UniqP = np.inner(HP,SGMT) 1633 Phi = np.inner(H,SSGT) 1515 1634 if SGInv: #if centro - expand HKL sets 1516 1635 Uniq = np.hstack((Uniq,-Uniq)) 1636 Uniq3 = np.hstack((Uniq3,-Uniq3)) 1517 1637 Phi = np.hstack((Phi,-Phi)) 1518 1638 UniqP = np.hstack((UniqP,-UniqP)) … … 1527 1647 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1528 1648 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0) 1529 phase = twopi*(np.inner(Uniq [:,:,:3],(dXdata.T+Xdata.T))-Phi[:,:,nxs])1649 phase = twopi*(np.inner(Uniq3,(dXdata.T+Xdata.T))-Phi[:,nxs,:,nxs]) 1530 1650 sinp = np.sin(phase) 1531 1651 cosp = np.cos(phase) 1532 1652 biso = -SQfactor*Uisodata[:,nxs] 1533 1653 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T 1534 HbH = -np.sum(UniqP[:,:, nxs,:]*np.inner(UniqP[:,:,:],bij),axis=-1) #use hklt proj to hkl1654 HbH = -np.sum(UniqP[:,:,:,nxs]*np.inner(UniqP[:,:,:],bij),axis=-1) #use hklt proj to hkl 1535 1655 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1536 1656 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1] #refBlk x ops x atoms 1657 # GSASIIpath.IPyBreak() 1537 1658 if 'T' in calcControls[hfx+'histType']: 1538 1659 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) … … 1541 1662 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 1542 1663 fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) 1543 GfpuA = G2mth.Modulation (Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms1664 GfpuA = G2mth.ModulationTw(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms 1544 1665 fag = fa*GfpuA[0]-fb*GfpuA[1] #real; 2 x refBlk x sym x atoms 1545 1666 fbg = fb*GfpuA[0]+fa*GfpuA[1] 1546 1667 fas = np.sum(np.sum(fag,axis=-1),axis=-1) #2 x refBlk; sum sym & atoms 1547 1668 fbs = np.sum(np.sum(fbg,axis=-1),axis=-1) 1548 # GSASIIpath.IPyBreak() 1549 if 'P' in calcControls[hfx+'histType']: 1550 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 #square of sums 1551 # refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) 1552 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1553 else: 1554 if len(TwinLaw) > 1: 1555 refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2 #FcT from primary twin element 1556 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+ \ 1557 np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins 1558 refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f" 1559 else: 1560 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 #square of sums 1561 refl.T[8] = np.copy(refl.T[10]) 1562 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1669 refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2 #FcT from primary twin element 1670 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fas,axis=0)**2,axis=-1)+ \ 1671 np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins 1672 refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f" 1563 1673 iBeg += blkSize 1564 1674 print 'nRef %d time %.4f\r'%(nRef,time.time()-time0) … … 1698 1808 dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1699 1809 dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1700 # GSASIIpath.IPyBreak()1701 1810 if not SGData['SGInv']: #Flack derivative 1702 1811 dfadfl = np.sum(-FPP*Tcorr*sinp) … … 1705 1814 dfadfl = 1.0 1706 1815 dfbdfl = 1.0 1816 # GSASIIpath.IPyBreak() 1707 1817 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 1708 SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin])1709 SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin])1818 SA = fas[0]+fas[1] #float = A+A' 1819 SB = fbs[0]+fbs[1] #float = B+B' 1710 1820 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro? 1711 1821 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) … … 1788 1898 1789 1899 # GSASIIpath.IPyBreak() 1900 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 1790 1901 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1791 1902 dFdvDict[phfx+'BabU'] = dFdbab.T[1] … … 3120 3231 time0 = time.time() 3121 3232 if im: 3122 SStructureFactor(refDict, im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)3233 SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3123 3234 else: 3124 3235 StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict) … … 3828 3939 phfx = '%d:%d:'%(Phase['pId'],hId) 3829 3940 SGData = Phase['General']['SGData'] 3941 TwinLaw = calcControls[phfx+'TwinLaw'] 3830 3942 im = 0 3831 3943 if parmDict[phfx+'Scale'] < 0.: … … 3840 3952 time0 = time.time() 3841 3953 if im: 3842 SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3954 if len(TwinLaw) > 1: 3955 SStructureFactorTw(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3956 else: 3957 SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3843 3958 else: 3844 3959 StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset
for help on using the changeset viewer.