Changeset 3754


Ignore:
Timestamp:
Dec 6, 2018 11:52:54 AM (4 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

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r3738 r3754  
    13041304    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods as betaij
    13051305    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods as betaij
    1306     nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]]
     1306    Am = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
     1307    Bm = np.array(MSSdata[3:]).T   #...cos pos mods
     1308    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]]
    13071309    if nWaves[0]:
    13081310        tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl
    … …  
    13431345    else:
    13441346        Umod = 1.0
    1345     Mmod = 1.0
     1347    if nWaves[3]:
     1348        tauM = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Mwaves x ngl
     1349        MmodA = Am[:,:,:,nxs]*np.sin(twopi*tauM[nxs,:,nxs,:]) #atoms X waves X 3 X ngl
     1350        MmodB = Bm[:,:,:,nxs]*np.cos(twopi*tauM[nxs,:,nxs,:]) #ditto
     1351        Mmod = np.sum(MmodA+MmodB,axis=1)                #atoms X 3 X ngl; sum waves
     1352        Mmod = np.swapaxes(Mmod,1,2)
     1353    else:
     1354        Mmod = 1.0
    13461355    return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt
    13471356       
    1348 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
     1357def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
    13491358    '''
    13501359    H: array nRefBlk x ops X hklt
    … …  
    13531362    Xmod: array atoms X 3 X ngl
    13541363    Umod: array atoms x 3x3 x ngl
     1364    Mmod: array atoms X 3 X ngl         #TODO: add mag moment modulation math
    13551365    glTau,glWt: arrays Gauss-Lorentzian pos & wts
    13561366    '''
    … …  
    14171427    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    14181428    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1419     nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]]
     1429    Am = np.array(MSSdata[:3]).T   #...cos pos mods x waves x atoms
     1430    Bm = np.array(MSSdata[3:]).T   #...cos pos mods
     1431    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]]
    14201432    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
    14211433    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
    … …  
    14631475        UmodA = 0.
    14641476        UmodB = 0.
    1465     return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
     1477    if nWaves[3]:
     1478        tauM = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
     1479        StauM = np.ones_like(Am)[:,:,:,nxs]*np.sin(twopi*tauM)[nxs,:,nxs,:]   #also dMmodA/dAm
     1480        CtauM = np.ones_like(Bm)[:,:,:,nxs]*np.cos(twopi*tauM)[nxs,:,nxs,:]   #also dMmodB/dBm
     1481        MmodA = Am[:,:,:,nxs]*StauM #atoms x waves x 3x3 x ngl
     1482        MmodB = Bm[:,:,:,nxs]*CtauM #ditto
     1483    else:
     1484        StauM = 1.0
     1485        CtauM = 1.0
     1486        MmodA = 0.
     1487        MmodB = 0.
     1488    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB,[StauM,CtauM],MmodA+MmodB
    14661489   
    14671490def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
    … …  
    16611684            tauT = G2spc.getTauT(tau,sop,ssop,drxyz,modul)[-1]
    16621685            tauT *= icent       #invert wave on -1
     1686#            print(tau,tauT,opr,G2spc.MT2text(sop).replace(' ',''),G2spc.SSMT2text(ssop).replace(' ',''))
    16631687            wave = np.zeros(3)
    16641688            uwave = np.zeros(6)

  • trunk/GSASIIspc.py

    r3753 r3754  
    35603560    else:
    35613561        newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC]
    3562 #        print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom)
     3562#        print(len(SGOps),Ax[0],iAx,nC,nA,NA,MT2text([M,T]).replace(' ',''),SGData['SpnFlp'][NA],Mom,newMom)
     3563#    print(Mom,newMom,MT2text([M,T]),)
    35633564    return newMom
    35643565       

  • trunk/GSASIIstrIO.py

    r3752 r3754  
    23952395                                ext = G2spc.checkMagextc([h,k,l],SGData)
    23962396                            mul *= 2      # for powder overlap of Friedel pairs
    2397                             if ext:# and not useExt:
     2397                            if ext and not useExt:
    23982398                                continue
    23992399                            if 'C' in inst['Type'][0]:

  • 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.