Ignore:
Timestamp:
Dec 6, 2018 11:52:54 AM (3 years ago)
Author:
vondreele
Message:

begin work on incommensurate mag. structure factors
fix problem of magnetic reflection extinctions
replace analytical derivatives of magnetic moments with numerical ones using deltM = 0.0001

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r3595 r3754  
    10591059        iBeg += blkSize
    10601060#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
     1061    return copy.deepcopy(refDict['RefList'])
    10611062
     1063def MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     1064    '''Compute magnetic structure factor derivatives numerically on blocks of reflections - for powders/nontwins only
     1065    input:
     1066   
     1067    :param dict refDict: where
     1068        'RefList' list where each ref = h,k,l,it,d,...
     1069        'FF' dict of form factors - filled in below
     1070    :param np.array G:      reciprocal metric tensor
     1071    :param str hfx:    histogram id string
     1072    :param str pfx:    phase id string
     1073    :param dict SGData: space group info. dictionary output from SpcGroup
     1074    :param dict calcControls:
     1075    :param dict parmDict:
     1076   
     1077    :returns: dict dFdvDict: dictionary of magnetic derivatives
     1078    '''
     1079   
     1080    trefDict = copy.deepcopy(refDict)
     1081    dM = 0.0001
     1082    dFdvDict = {}
     1083    for parm in parmDict:
     1084        if 'AM' in parm:
     1085            parmDict[parm] += dM
     1086            prefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
     1087            parmDict[parm] -= 2*dM
     1088            mrefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict)
     1089            parmDict[parm] += dM
     1090            dFdvDict[parm] = (prefList[:,9]-mrefList[:,9])/(2.*dM)
     1091    return dFdvDict
     1092           
    10621093def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    1063     '''Compute magnetic structure factor derivatives on blocks of reflections - for powders/nontwins only
     1094    '''Compute nonmagnetic structure factor derivatives on blocks of reflections in magnetic structures - for powders/nontwins only
    10641095    input:
    10651096   
     
    10981129    mSize = len(Mdata)
    10991130    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
    1100     dMdm = np.inner(Gdata.T,Ginv).T/Mag
    11011131    Gones = np.ones_like(Gdata)
    11021132    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     
    11121142    VGi = np.sqrt(nl.det(Ginv))
    11131143    Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag       #make unit vectors in Cartesian space
    1114     dkdG = (np.inner(Gones.T,uAmat).T*VGi)/Mag
    1115     dkdm = dkdG-Kdata*dMdm[:,nxs,:]/Mag[nxs,:,:]
    1116     dFdMx = np.zeros((nRef,mSize,3))
    11171144    Uij = np.array(G2lat.U6toUij(Uijdata))
    11181145    bij = Mast*Uij.T
     
    11201147    dFdfr = np.zeros((nRef,mSize))
    11211148    dFdx = np.zeros((nRef,mSize,3))
    1122     dFdMx = np.zeros((3,nRef,mSize))
    11231149    dFdui = np.zeros((nRef,mSize))
    11241150    dFdua = np.zeros((nRef,mSize,6))
     
    11651191        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
    11661192        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    1167         dqdk = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T     #Mxyz**2,Nref
    1168         dqdm = dqdk[:,:,:,nxs,nxs]*dkdm[:,nxs,nxs,:,:]   #Mxyz**2,Nref,Nops,Natms
    1169         dmx = Q*dMdm[:,nxs,nxs,:]
    1170         dmx = dmx[nxs,:,:,:,:]+dqdm*Mag[nxs,nxs,nxs,:,:]
    1171         dmx /= 2.
    11721193       
    11731194        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
     
    11801201        dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK
    11811202        dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK
    1182         dfadmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*cosm[nxs,nxs,:,:,:],axis=3)
    11831203        dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK
    11841204        dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK
     
    11861206        dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem
    11871207        dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2)
    1188         dfbdmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*sinm[nxs,nxs,:,:,:],axis=3)
    11891208        dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms)
    11901209        dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2)
     
    11921211        dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok
    11931212        dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok
    1194         dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)          #problems
    11951213        dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok
    11961214        dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok
     
    12031221        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
    12041222        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
    1205         dFdvDict[pfx+'AMx:'+str(i)] = dFdMx[0,:,i]
    1206         dFdvDict[pfx+'AMy:'+str(i)] = dFdMx[1,:,i]
    1207         dFdvDict[pfx+'AMz:'+str(i)] = dFdMx[2,:,i]
    12081223        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
    12091224        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     
    15251540                fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    15261541                fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1527             GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms
     1542            GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x refBlk x sym X atoms
    15281543            fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms
    15291544            fbg = fb*GfpuA[0]+fa*GfpuA[1]
     
    17161731    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    17171732    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
    1718     waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
     1733    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5]
    17191734    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    17201735    FF = np.zeros(len(Tdata))
     
    17801795        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
    17811796        fotp = FPP*Tcorr            #ops x atoms
    1782         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
     1797        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms
    17831798        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    17841799        # GfpuA is 2 x ops x atoms
     
    19301945    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    19311946    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
    1932     waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
     1947    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5]
    19331948    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    19341949    FF = np.zeros(len(Tdata))
     
    20052020        fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms
    20062021        fotp = FPP*Tcorr            #ops x atoms
    2007         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
     2022        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms
    20082023        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    20092024        # GfpuA is 2 x ops x atoms
     
    21602175    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    21612176    ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
    2162     waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)
     2177    waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5]
    21632178    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    21642179    FF = np.zeros(len(Tdata))
     
    22322247        fot = (FF+FP-Bab)*Tcorr     #ops x atoms
    22332248        fotp = FPP*Tcorr            #ops x atoms
    2234         GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms
     2249        GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms
    22352250        dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)
    22362251        # GfpuA is 2 x ops x atoms
     
    34943509                dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    34953510            else:
    3496                 if Phase['General']['Type'] == 'magnetic':
     3511                if Phase['General']['Type'] == 'magnetic': 
    34973512                    dFdvDict = MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3513                    dFdvDict.update(MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict))
    34983514                else:
    34993515                    dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.