Changeset 1921 for trunk/GSASIIstrMath.py
 Timestamp:
 Jul 6, 2015 12:48:05 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIstrMath.py
r1919 r1921 631 631 return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze() 632 632 633 def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):634 ''' Not Used: here only for comparison the StructureFactor2  faster version635 Compute structure factors for all h,k,l for phase636 puts the result, F^2, in each ref[8] in refList637 input:638 639 :param dict refDict: where640 'RefList' list where each ref = h,k,l,m,d,...641 'FF' dict of form factors  filed in below642 :param np.array G: reciprocal metric tensor643 :param str pfx: phase id string644 :param dict SGData: space group info. dictionary output from SpcGroup645 :param dict calcControls:646 :param dict ParmDict:647 648 '''649 phfx = pfx.split(':')[0]+hfx650 ast = np.sqrt(np.diag(G))651 Mast = twopisq*np.multiply.outer(ast,ast)652 SGMT = np.array([ops[0].T for ops in SGData['SGOps']])653 SGT = np.array([ops[1] for ops in SGData['SGOps']])654 FFtables = calcControls['FFtables']655 BLtables = calcControls['BLtables']656 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)657 FF = np.zeros(len(Tdata))658 if 'NC' in calcControls[hfx+'histType']:659 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])660 else:661 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])662 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])663 Uij = np.array(G2lat.U6toUij(Uijdata))664 bij = Mast*Uij.T665 if not len(refDict['FF']):666 if 'N' in calcControls[hfx+'histType']:667 dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's668 else:669 dat = G2el.getFFvalues(FFtables,0.)670 refDict['FF']['El'] = dat.keys()671 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))672 for iref,refl in enumerate(refDict['RefList']):673 if 'NT' in calcControls[hfx+'histType']:674 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14])675 fbs = np.array([0,0])676 H = refl[:3]677 SQ = 1./(2.*refl[4])**2678 SQfactor = 4.0*SQ*twopisq679 Bab = parmDict[phfx+'BabA']*np.exp(parmDict[phfx+'BabU']*SQfactor)680 if not np.any(refDict['FF']['FF'][iref]): #no form factors  1st time thru StructureFactor681 if 'N' in calcControls[hfx+'histType']:682 dat = G2el.getBLvalues(BLtables)683 refDict['FF']['FF'][iref] = dat.values()684 else: #'X'685 dat = G2el.getFFvalues(FFtables,SQ)686 refDict['FF']['FF'][iref] = dat.values()687 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])688 FF = refDict['FF']['FF'][iref][Tindx]689 Uniq = np.inner(H,SGMT)690 Phi = np.inner(H,SGT)691 phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])692 sinp = np.sin(phase)693 cosp = np.cos(phase)694 biso = SQfactor*Uisodata695 Tiso = np.where(biso<1.,np.exp(biso),1.0)696 HbH = np.array([np.inner(h,np.inner(bij,h)) for h in Uniq])697 Tuij = np.where(HbH<1.,np.exp(HbH),1.0)698 Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)699 fa = np.array([(FF+FPBab)*cosp*Tcorr,FPP*sinp*Tcorr])700 fas = np.sum(np.sum(fa,axis=1),axis=1) #real701 if not SGData['SGInv']:702 fb = np.array([(FF+FPBab)*sinp*Tcorr,FPP*cosp*Tcorr])703 fbs = np.sum(np.sum(fb,axis=1),axis=1)704 fasq = fas**2705 fbsq = fbs**2 #imaginary706 refl[9] = np.sum(fasq)+np.sum(fbsq)707 refl[10] = atan2d(fbs[0],fas[0])708 709 633 def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 710 634 ''' … … 1038 962 dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)] 1039 963 dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)] 1040 dFdtw[iref] = 2.*TwMask*np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2964 dFdtw[iref] = TwMask*np.sum(fas,axis=0)**2+TwMask*np.sum(fbs,axis=0)**2 1041 965 else: #these are good for no twin single crystals 1042 966 dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \ … … 1080 1004 if nTwin > 1: 1081 1005 for i in range(nTwin): 1082 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] /4.1006 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1083 1007 return dFdvDict 1084 1008
Note: See TracChangeset
for help on using the changeset viewer.