Changeset 3754 for trunk/GSASIImath.py


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/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)
Note: See TracChangeset for help on using the changeset viewer.