Changeset 2478


Ignore:
Timestamp:
Sep 23, 2016 1:43:24 PM (6 years ago)
Author:
vondreele
Message:

work up to magnetic structure factor calcs.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIElem.py

    r2473 r2478  
    6060    return FFtable
    6161   
    62 def GetMFtable(atomTypes):
     62def GetMFtable(atomTypes,Landeg):
    6363    ''' returns a dictionary of magnetic form factor data for atom types found in atomTypes
    6464
    6565    :param list atomTypes: list of atom types
     66    :param list Landeg: Lande g factors for atomTypes
    6667    :return: FFtable, dictionary of form factor data; key is atom type
    6768
    6869    '''
    69     FFtable = {}
    70     for El in atomTypes:
    71         FFs = GetMagFormFacCoeff(getElSym(El))
    72         for item in FFs:
     70    MFtable = {}
     71    for El,gfac in zip(atomTypes,Landeg):
     72        MFs = GetMagFormFacCoeff(getElSym(El))
     73        for item in MFs:
    7374            if item['Symbol'] == El.upper():
    74                 FFtable[El] = item
    75     return FFtable
     75                item['gfac'] = gfac
     76                MFtable[El] = item
     77    return MFtable
    7678   
    7779def GetBLtable(General):
    … …  
    122124    return BLvals
    123125       
     126def getMFvalues(MFtables,SQ,ifList=False):
     127    'Needs a doc string'
     128    if ifList:
     129        MFvals = []
     130        for El in MFtables:
     131            MFvals.append(MagScatFac(MFtables[El],SQ)[0])
     132    else:
     133        MFvals = {}
     134        for El in MFtables:
     135            MFvals[El] = MagScatFac(MFtables[El],SQ)[0]
     136    return MFvals
     137   
    124138def GetFFC5(ElSym):
    125139    '''Get 5 term form factor and Compton scattering data
    … …  
    327341    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
    328342       
    329 def MagScatFac(El, SQ,gfac):
     343def MagScatFac(El, SQ):
    330344    """compute value of form factor
    331345
    … …  
    343357    MMF = np.sum(mfa[:,np.newaxis]*np.exp(mt)[:],axis=0)+El['mfc']
    344358    NMF = np.sum(nfa[:,np.newaxis]*np.exp(nt)[:],axis=0)+El['nfc']
    345     return MMF+(2.0/gfac-1.0)*NMF
     359    return MMF+(2.0/El['gfac']-1.0)*NMF
    346360       
    347361def BlenResCW(Els,BLtables,wave):

  • trunk/GSASIIstrIO.py

    r2474 r2478  
    190190    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    191191    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
    192     Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
     192    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
     193        GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
    193194    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
    194195    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
    … …  
    10971098        FFtables.update(FFtable)
    10981099        BLtables.update(BLtable)
     1100        phaseDict[pfx+'isMag'] = False
    10991101        if General['Type'] == 'magnetic':
    1100             MFtable = G2el.GetMFtable(General['AtomTypes'])
     1102            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
    11011103            MFtables.update(MFtable)
     1104            phaseDict[pfx+'isMag'] = True
    11021105        Atoms = PhaseData[name]['Atoms']
    11031106        if Atoms and not General.get('doPawley'):
    … …  
    13551358            phaseVary += pawleyVary
    13561359               
    1357     return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave
     1360    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
    13581361   
    13591362def cellFill(pfx,SGData,parmDict,sigDict):
    … …  
    22762279                    Uniq = []
    22772280                    Phi = []
     2281                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
    22782282                    if Phases[phase]['General'].get('Modulated',False):
    22792283                        ifSuper = True
    … …  
    22832287                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    22842288                            mul *= 2      # for powder overlap of Friedel pairs
    2285                             if m or not ext:
     2289                            if m or not ext or useExt:
    22862290                                if 'C' in inst['Type'][0]:
    22872291                                    pos = G2lat.Dsp2pos(inst,d)
    … …  
    23062310                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    23072311                            mul *= 2      # for powder overlap of Friedel pairs
    2308                             if ext:
     2312                            if ext and not useExt:
    23092313                                continue
    23102314                            if 'C' in inst['Type'][0]:

  • trunk/GSASIIstrMain.py

    r2466 r2478  
    159159    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    160160    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
    161     Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
     161    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
     162        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
    162163    calcControls['atomIndx'] = atomIndx
    163164    calcControls['Natoms'] = Natoms
    164165    calcControls['FFtables'] = FFtables
    165166    calcControls['BLtables'] = BLtables
     167    calcControls['MFtables'] = MFtables
    166168    calcControls['maxSSwave'] = maxSSwave
    167169    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
    … …  
    264266    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
    265267    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
    266     Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
     268    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
     269        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
    267270    for item in phaseVary:
    268271        if '::A0' in item:
    … …  
    294297        calcControls['FFtables'] = FFtables
    295298        calcControls['BLtables'] = BLtables
     299        calcControls['MFtables'] = MFtables
    296300        calcControls['maxSSwave'] = maxSSwave
    297301        Histo = {histogram:Histograms[histogram],}

  • trunk/GSASIIstrMath.py

    r2405 r2478  
    593593    Uisodata = np.zeros(Natoms)
    594594    Uijdata = np.zeros((6,Natoms))
     595    Gdata = np.zeros((3,Natoms))
    595596    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
    596597        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
    597598        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
    598599        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
    599         'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
     600        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5],
     601        'AMx':Gdata[0],'AMy':Gdata[1],'AMz':Gdata[2],}
    600602    for iatm in range(Natoms):
    601603        for key in keys:
    … …  
    604606                keys[key][iatm] = parmDict[parm]
    605607    Fdata = np.where(Fdata,Fdata,1.e-8)         #avoid divide by zero in derivative calc.
    606     return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
     608    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata
    607609   
    608610def GetAtomSSFXU(pfx,calcControls,parmDict):
    … …  
    668670    FFtables = calcControls['FFtables']
    669671    BLtables = calcControls['BLtables']
     672    MFtables = calcControls['MFtables']
    670673    Flack = 1.0
    671674    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    … …  
    679682        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
    680683        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    681     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    682     FF = np.zeros(len(Tdata))
     684    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     685        GetAtomFXU(pfx,calcControls,parmDict)
    683686    if 'NC' in calcControls[hfx+'histType']:
    684687        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
    … …  
    691694    nRef = refDict['RefList'].shape[0]
    692695    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
     696        SQ = 1./(2.*refDict['RefList'].T[4])**2
    693697        if 'N' in calcControls[hfx+'histType']:
    694698            dat = G2el.getBLvalues(BLtables)
    695699            refDict['FF']['El'] = dat.keys()
    696             refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
     700            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
     701            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
     702            for iel,El in enumerate(refDict['FF']['El']):
     703                if El in MFtables:
     704                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    697705        else:       #'X'
    698706            dat = G2el.getFFvalues(FFtables,0.)
    699707            refDict['FF']['El'] = dat.keys()
    700             refDict['FF']['FF'] = np.ones((nRef,len(dat)))
    701             for iref,ref in enumerate(refDict['RefList']):
    702                 SQ = 1./(2.*ref[4])**2
    703                 dat = G2el.getFFvalues(FFtables,SQ)
    704                 refDict['FF']['FF'][iref] *= dat.values()
     708            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
     709            for iel,El in enumerate(refDict['FF']['El']):
     710                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
     711#    GSASIIpath.IPyBreak()
    705712#reflection processing begins here - big arrays!
    706713    iBeg = 0
    … …  
    731738        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    732739        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    733         FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
    734740        Uniq = np.inner(H,SGMT)
    735741        Phi = np.inner(H,SGT)
    … …  
    742748        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    743749        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
     750        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     751        if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
     752            MF = np.repeat(refDict['FF']['MF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     753            #calc fam & fbm here
    744754        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    745755            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    … …  
    754764            fas[1] *= 0.
    755765        if 'PWDR' in calcControls[hfx+'histType']:     #PWDR: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
    756             refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)    
     766            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
    757767            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    758768        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    … …  
    805815#    nTwin = len(TwinLaw)       
    806816#    nRef = len(refDict['RefList'])
    807 #    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     817#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     818#        GetAtomFXU(pfx,calcControls,parmDict)
    808819#    mSize = len(Mdata)
    809820#    FF = np.zeros(len(Tdata))
    … …  
    10011012    BLtables = calcControls['BLtables']
    10021013    nRef = len(refDict['RefList'])
    1003     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1014    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1015        GetAtomFXU(pfx,calcControls,parmDict)
    10041016    mSize = len(Mdata)
    10051017    FF = np.zeros(len(Tdata))
    … …  
    11711183#    nTwin = len(TwinLaw)       
    11721184#    nRef = len(refDict['RefList'])
    1173 #    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1185#    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1186#        GetAtomFXU(pfx,calcControls,parmDict)
    11741187#    mSize = len(Mdata)
    11751188#    FF = np.zeros(len(Tdata))
    … …  
    13171330    nTwin = len(TwinLaw)       
    13181331    nRef = len(refDict['RefList'])
    1319     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1332    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1333        GetAtomFXU(pfx,calcControls,parmDict)
    13201334    mSize = len(Mdata)
    13211335    FF = np.zeros(len(Tdata))
    … …  
    14761490    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    14771491        Flack = 1.-2.*parmDict[phfx+'Flack']
    1478     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1492    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1493        GetAtomFXU(pfx,calcControls,parmDict)
    14791494    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    14801495    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
    … …  
    14911506    nRef = refDict['RefList'].shape[0]
    14921507    if not len(refDict['FF']):
     1508        SQ = 1./(2.*refDict['RefList'].T[5])**2
    14931509        if 'N' in calcControls[hfx+'histType']:
    1494             dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1510            dat = G2el.getBLvalues(BLtables)
    14951511            refDict['FF']['El'] = dat.keys()
    1496             refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
     1512            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
     1513            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
     1514            for iel,El in enumerate(refDict['FF']['El']):
     1515                if El in MFtables:
     1516                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    14971517        else:
    1498             dat = G2el.getFFvalues(FFtables,0.)       
     1518            dat = G2el.getFFvalues(FFtables,0.)
    14991519            refDict['FF']['El'] = dat.keys()
    1500             refDict['FF']['FF'] = np.ones((nRef,len(dat)))
    1501             for iref,ref in enumerate(refDict['RefList']):
    1502                 SQ = 1./(2.*ref[5])**2
    1503                 dat = G2el.getFFvalues(FFtables,SQ)
    1504                 refDict['FF']['FF'][iref] *= dat.values()
     1520            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
     1521            for iel,El in enumerate(refDict['FF']['El']):
     1522                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
    15051523    time0 = time.time()
    15061524#reflection processing begins here - big arrays!
    … …  
    15991617        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
    16001618        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    1601     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1619    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1620        GetAtomFXU(pfx,calcControls,parmDict)
    16021621    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    16031622    ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast)
    … …  
    16131632    blkSize = 32       #no. of reflections in a block
    16141633    nRef = refDict['RefList'].shape[0]
    1615     if not len(refDict['FF']):
     1634    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
     1635        SQ = 1./(2.*refDict['RefList'].T[5])**2
    16161636        if 'N' in calcControls[hfx+'histType']:
    1617             dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1637            dat = G2el.getBLvalues(BLtables)
    16181638            refDict['FF']['El'] = dat.keys()
    1619             refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
     1639            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()
     1640            refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
     1641            for iel,El in enumerate(refDict['FF']['El']):
     1642                if El in MFtables:
     1643                    refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    16201644        else:
    1621             dat = G2el.getFFvalues(FFtables,0.)       
     1645            dat = G2el.getFFvalues(FFtables,0.)
    16221646            refDict['FF']['El'] = dat.keys()
    1623             refDict['FF']['FF'] = np.ones((nRef,len(dat)))
    1624             for iref,ref in enumerate(refDict['RefList']):
    1625                 SQ = 1./(2.*ref[5])**2
    1626                 dat = G2el.getFFvalues(FFtables,SQ)
    1627                 refDict['FF']['FF'][iref] *= dat.values()
     1647            refDict['FF']['FF'] = np.zeros((nRef,len(dat)))
     1648            for iel,El in enumerate(refDict['FF']['El']):
     1649                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
    16281650    time0 = time.time()
    16291651#reflection processing begins here - big arrays!
    … …  
    17251747    BLtables = calcControls['BLtables']
    17261748    nRef = len(refDict['RefList'])
    1727     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1749    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1750        GetAtomFXU(pfx,calcControls,parmDict)
    17281751    mSize = len(Mdata)  #no. atoms
    17291752    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    … …  
    19371960    BLtables = calcControls['BLtables']
    19381961    nRef = len(refDict['RefList'])
    1939     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1962    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     1963        GetAtomFXU(pfx,calcControls,parmDict)
    19401964    mSize = len(Mdata)  #no. atoms
    19411965    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    … …  
    21672191    nTwin = len(TwinLaw)       
    21682192    nRef = len(refDict['RefList'])
    2169     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     2193    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     2194        GetAtomFXU(pfx,calcControls,parmDict)
    21702195    mSize = len(Mdata)  #no. atoms
    21712196    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.