Changeset 4054


Ignore:
Timestamp:
Jul 6, 2019 7:00:41 PM (2 years ago)
Author:
vondreele
Message:

changes to incommensurate magnetism math

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r4025 r4054  
    14021402    return Mmod,MmodA,MmodB    #Ntau,Nops,Natm,,Mxyz; sum,sin & cos parts
    14031403       
     1404def MagMod2(XYZ,modQ,MSSdata,SGData,SSGData):
     1405    '''
     1406    this needs to make magnetic moment modulations & magnitudes
     1407    '''
     1408    Am = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
     1409    Bm = np.array(MSSdata[3:]).T  #...cos pos mods
     1410    nWaves = Am.shape[1]
     1411    if not nWaves:
     1412        return 0.0,0.0
     1413    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
     1414    Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
     1415    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1416    if SGData['SGInv']:
     1417        SGMT = np.vstack((SGMT,-SGMT))
     1418        Sinv =np.vstack((Sinv,-Sinv))
     1419        SGT = np.vstack((SGT,-SGT))
     1420    SGMT = np.vstack([SGMT for cen in SGData['SGCen']])     #Nops,3,3
     1421    Sinv = np.vstack([Sinv for cen in SGData['SGCen']])     #Nops,4,4
     1422    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
     1423    detSM = nl.det(SGMT)            #Nops
     1424    epsinv = Sinv[:,3,3]            #Nops
     1425    kdr = np.inner(XYZ,modQ).T      #Nops,Natm
     1426    phase = kdr+(epsinv*(np.inner(modQ,SGT[:,:3])-SGT[:,3]))[:,nxs]     #Nops,Natm
     1427   
     1428    psin = np.sin(twopi*phase)      #Nops,Natm
     1429    pcos = np.cos(twopi*phase)
     1430    MmodB = np.sum(Bm[nxs,:,:,:]*pcos[:,:,nxs,nxs],axis=2)      #Nops,Natm,3
     1431    MmodA = np.sum(Am[nxs,:,:,:]*psin[:,:,nxs,nxs],axis=2)
     1432#    if SGData['SGGray']:
     1433#        MmodA = -np.sum(SGMT[:,nxs,:,:]*MmodA[:,:,nxs,:],axis=-1)*detSM[:,nxs,nxs]*SGData['SpnFlp'][:,nxs,nxs]
     1434#        MmodB = -np.sum(SGMT[:,nxs,:,:]*MmodB[:,:,nxs,:],axis=-1)*detSM[:,nxs,nxs]*SGData['SpnFlp'][:,nxs,nxs]
     1435#    else:
     1436    MmodA = np.sum(SGMT[:,nxs,:,:]*MmodA[:,:,nxs,:],axis=-1)*SGData['MagMom'][:,nxs,nxs]
     1437    MmodB = np.sum(SGMT[:,nxs,:,:]*MmodB[:,:,nxs,:],axis=-1)*SGData['MagMom'][:,nxs,nxs]
     1438    Mmod = MmodA+MmodB
     1439    return Mmod,MmodA,MmodB    #Nops,Natm,,Mxyz; sum,sin & cos parts
     1440       
    14041441def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    14051442    '''
     
    14311468    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    14321469   
    1433 def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
    1434     '''
    1435     H: array nRefBlk x ops X hklt
    1436     HP: array nRefBlk x ops X hklt proj to hkl
    1437     nWaves: list number of waves for frac, pos, uij & mag
    1438     Fmod: array 2 x atoms x waves    (sin,cos terms)
    1439     Xmod: array atoms X 3 X ngl
    1440     Umod: array atoms x 3x3 x ngl
    1441     glTau,glWt: arrays Gauss-Lorentzian pos & wts
    1442     '''
    1443    
    1444     if nWaves[2]:       #uij (adp) waves
    1445         if len(HP.shape) > 2:
    1446             HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
    1447         else:
    1448             HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
    1449     else:
    1450         HbH = 1.0
    1451     HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
    1452     if len(H.shape) > 2:
    1453         D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
    1454         HdotXD = twopi*(HdotX+D[:,:,nxs,:])
    1455     else:
    1456         D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
    1457         HdotXD = twopi*(HdotX+D[:,nxs,:])
    1458     M = np.swapaxes(Mmod,1,2)
    1459     cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
    1460     sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
    1461     return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    1462 
     1470#def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
     1471#    '''
     1472#    H: array nRefBlk x ops X hklt
     1473#    HP: array nRefBlk x ops X hklt proj to hkl
     1474#    nWaves: list number of waves for frac, pos, uij & mag
     1475#    Fmod: array 2 x atoms x waves    (sin,cos terms)
     1476#    Xmod: array atoms X 3 X ngl
     1477#    Umod: array atoms x 3x3 x ngl
     1478#    glTau,glWt: arrays Gauss-Lorentzian pos & wts
     1479#    '''
     1480#   
     1481#    if nWaves[2]:       #uij (adp) waves
     1482#        if len(HP.shape) > 2:
     1483#            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1484#        else:
     1485#            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1486#    else:
     1487#        HbH = 1.0
     1488#    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
     1489#    if len(H.shape) > 2:
     1490#        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
     1491#        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
     1492#    else:
     1493#        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
     1494#        HdotXD = twopi*(HdotX+D[:,nxs,:])
     1495#    M = np.swapaxes(Mmod,1,2)
     1496#    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
     1497#    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
     1498#    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
     1499#
    14631500def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    14641501    '''
  • trunk/GSASIIstrMath.py

    r4051 r4054  
    15101510       
    15111511        mXYZ = np.array([[xyz[0] for xyz in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T])%1. #Natn,Nop,xyz
    1512         Tmag,TmagA,TmagB = G2mth.MagMod(ngl,mXYZ,modQ,MSSdata,SGData,SSGData)   #Ntau,Nops,Natm,Mxyz-Tmag matches drawing moments
     1512        Tmag,TmagA,TmagB = G2mth.MagMod2(mXYZ,modQ,MSSdata,SGData,SSGData)   #Nops,Natm,Mxyz-Tmag matches drawing moments
    15131513       
    15141514        if not SGData['SGGray']:    #for fixed Mx,My,Mz
     
    16121612            M = np.array(np.abs(H[3]),dtype=np.int)-1
    16131613       
    1614             fam = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*cosm[i,nxs,:,:,nxs]-    \
    1615                 np.sign(H[3,i])*TmagA*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])
    1616             fbm = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*sinm[i,nxs,:,:,nxs]+    \
    1617                 np.sign(H[3,i])*TmagA*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])
     1614            fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*cosm[i,:,:,nxs]-    \
     1615                np.sign(H[3,i])*TmagA*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
     1616            fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(TmagB*sinm[i,:,:,nxs]+    \
     1617                np.sign(H[3,i])*TmagA*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
    16181618                       
    16191619            famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
    16201620            fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
    16211621           
    1622             fa = np.sum(famq,axis=-1)**2-np.sum(eM.T[:,nxs,:]*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
    1623             fb = np.sum(fbmq,axis=-1)**2-np.sum(eM.T[:,nxs,:]*fbmq,axis=-1)**2
    1624            
    1625             fas = np.sum(fa*glWt,axis=1)
    1626             fbs = np.sum(fb*glWt,axis=1)
     1622            fas = np.sum(famq,axis=-1)**2-np.sum(eM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
     1623            fbs = np.sum(fbmq,axis=-1)**2-np.sum(eM.T*fbmq,axis=-1)**2
    16271624           
    16281625            refl.T[10] = np.where(H[3],fas+fbs,fas0+fbs0)
Note: See TracChangeset for help on using the changeset viewer.