Ignore:
Timestamp:
Jul 6, 2015 12:48:05 PM (7 years ago)
Author:
vondreele
Message:

twin derivatives now OK.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1919 r1921  
    631631    return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()   
    632632   
    633 def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    634     ''' Not Used: here only for comparison the StructureFactor2 - faster version
    635     Compute structure factors for all h,k,l for phase
    636     puts the result, F^2, in each ref[8] in refList
    637     input:
    638    
    639     :param dict refDict: where
    640         'RefList' list where each ref = h,k,l,m,d,...
    641         'FF' dict of form factors - filed in below
    642     :param np.array G:      reciprocal metric tensor
    643     :param str pfx:    phase id string
    644     :param dict SGData: space group info. dictionary output from SpcGroup
    645     :param dict calcControls:
    646     :param dict ParmDict:
    647 
    648     '''       
    649     phfx = pfx.split(':')[0]+hfx
    650     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.T
    665     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's
    668         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])**2
    678         SQfactor = 4.0*SQ*twopisq
    679         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 StructureFactor
    681             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*Uisodata
    695         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+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
    700         fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    701         if not SGData['SGInv']:
    702             fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    703             fbs = np.sum(np.sum(fb,axis=1),axis=1)
    704         fasq = fas**2
    705         fbsq = fbs**2        #imaginary
    706         refl[9] = np.sum(fasq)+np.sum(fbsq)
    707         refl[10] = atan2d(fbs[0],fas[0])
    708    
    709633def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    710634    '''
     
    1038962                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)]
    1039963                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)**2
     964                dFdtw[iref] = TwMask*np.sum(fas,axis=0)**2+TwMask*np.sum(fbs,axis=0)**2
    1041965            else:   #these are good for no twin single crystals
    1042966                dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \
     
    10801004    if nTwin > 1:
    10811005        for i in range(nTwin):
    1082             dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]/4.
     1006            dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]
    10831007    return dFdvDict
    10841008   
Note: See TracChangeset for help on using the changeset viewer.