Changeset 3865


Ignore:
Timestamp:
Mar 31, 2019 6:59:15 PM (3 years ago)
Author:
vondreele
Message:

a new attempt at incommensurate mag. str. fctrs.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3864 r3865  
    13561356    return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt
    13571357
    1358 def MagMod(XYZ,modQ,MSSdata,SGData,SSGData):
    1359    
     1358def MagMod(ngl,XYZ,modQ,MSSdata,SGData,SSGData):
     1359    '''
     1360    this needs to make magnetic moment modulations & magnitudes as
     1361    fxn of ngl tau points
     1362    '''
     1363    Bm = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
     1364    Am = np.array(MSSdata[3:]).T   #...cos pos mods
     1365    nWaves = Am.shape[1]
     1366    ngl = 20
     1367    tau = np.arange(ngl)/ngl
     1368    if not nWaves:
     1369        return 0.0,0.0
    13601370    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
    13611371    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     
    13691379    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
    13701380    RI = np.hstack([RI for cen in SSGData['SSGCen']])
    1371     Bm = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
    1372     Am = np.array(MSSdata[3:]).T   #...cos pos mods
    1373     nWaves = Am.shape[1]
    1374 
    1375 #works for DyMn6Ge6 but not MnWO4
    1376     nCen = len(SSGData['SSGCen'])
    1377     nEqv = XYZ.shape[1]         #full no. of equivalent pos incl centering
    1378     mEqv = nEqv//nCen           #no. operators not with centering; includes inversion
    1379     MmodAA = 0; MmodBA = 0
    1380     if nWaves:
    1381         modind = np.arange(nWaves)+1.
    1382         phaseA = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ))).T         #= 2pimk.r Nops,Natm,Nwave
    1383         if nCen > 0:
    1384             phshp = phaseA.shape
    1385             phaseA = np.reshape(phaseA,(nCen,mEqv,phshp[1],-1))
    1386             for ic,cen in enumerate(SSGData['SSGCen']):
    1387                 phaseA[ic] += twopi*cen[3]/2.
    1388             phaseA = np.reshape(phaseA,phshp)               
    1389         psinA = np.sin(phaseA)
    1390         pcosA = np.cos(phaseA)
    1391         MmodAA = Am[nxs,:,:,:]*pcosA[:,:,:,nxs]-Bm[nxs,:,:,:]*psinA[:,:,:,nxs]
    1392         MmodBA = Am[nxs,:,:,:]*psinA[:,:,:,nxs]+Bm[nxs,:,:,:]*pcosA[:,:,:,nxs]   
    1393    
    1394 #fails   
     1381    modind = np.arange(nWaves)+1.
    13951382    modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = -1.
    13961383    toff = RI*np.inner(modQp,SGT)
    1397     MmodA = 0; MmodB = 0
    1398     if nWaves:
    1399         modind = np.arange(nWaves)+1.
    1400         phase = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ)-toff)).T
    1401         psin = np.sin(phase)
    1402         pcos = np.cos(phase)
    1403         RAM = np.rollaxis(np.inner(Am,SGMT),2)
    1404         RBM = np.rollaxis(np.inner(Bm,SGMT),2)
    1405         RAC = RAM*pcos[:,:,:,nxs]
    1406         RBS = RBM*psin[:,:,:,nxs]
    1407         RAS = RAM*psin[:,:,:,nxs]
    1408         RBC = RBM*pcos[:,:,:,nxs]
    1409         MmodA = RAC-RBS
    1410         MmodB = RAS+RBC
    1411         MmodA *= thdetR[:,nxs,nxs,nxs]
    1412         MmodB *= thdetR[:,nxs,nxs,nxs]
    1413 
    1414     return MmodAA,MmodBA    #cos & sin Nops,Natm,Nwaves,Mxyz
     1384    phase = (modind[:,nxs,nxs]*(np.inner(XYZ,modQ))).T     #Nops,Natm,Nwave
     1385#    phase = (thdetR[nxs,nxs,:]*(modind[:,nxs,nxs]*np.inner(XYZ,modQ))).T     #Nops,Natm,Nwave
     1386    phase = phase[nxs,:,:,:]+tau[:,nxs,nxs,nxs]                 #Ntau,Nops,Natm,Nwave
     1387    psin = np.sin(twopi*phase)
     1388    pcos = np.cos(twopi*phase)
     1389    Mmod = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,:,nxs]+Bm[nxs,nxs,:,:,:]*psin[:,:,:,:,nxs],axis=3)
     1390    return Mmod    #Ntau,Nops,Natm,,Mxyz
     1391
     1392
     1393
     1394##works for DyMn6Ge6 but not MnWO4
     1395#    nCen = len(SSGData['SSGCen'])
     1396#    nEqv = XYZ.shape[1]         #full no. of equivalent pos incl centering
     1397#    mEqv = nEqv//nCen           #no. operators not with centering; includes inversion
     1398#    MmodAA = 0; MmodBA = 0
     1399#    if nWaves:
     1400#        phaseA = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ))).T         #= 2pimk.r Nops,Natm,Nwave
     1401#        if nCen > 0:
     1402#            phshp = phaseA.shape
     1403#            phaseA = np.reshape(phaseA,(nCen,mEqv,phshp[1],-1))
     1404#            for ic,cen in enumerate(SSGData['SSGCen']):
     1405#                phaseA[ic] += twopi*cen[3]/2.
     1406#            phaseA = np.reshape(phaseA,phshp)               
     1407#        psinA = np.sin(phaseA)
     1408#        pcosA = np.cos(phaseA)
     1409#        MmodAA = Am[nxs,:,:,:]*pcosA[:,:,:,nxs]-Bm[nxs,:,:,:]*psinA[:,:,:,nxs]
     1410#        MmodBA = Am[nxs,:,:,:]*psinA[:,:,:,nxs]+Bm[nxs,:,:,:]*pcosA[:,:,:,nxs]   
     1411#   
     1412##fails   
     1413#    modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = -1.
     1414#    toff = RI*np.inner(modQp,SGT)
     1415#    MmodA = 0; MmodB = 0
     1416#    if nWaves:
     1417#        modind = np.arange(nWaves)+1.
     1418#        phase = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ)-toff)).T
     1419#        RAM = np.rollaxis(np.inner(Am,SGMT),2)
     1420#        RBM = np.rollaxis(np.inner(Bm,SGMT),2)
     1421#        RAC = RAM*pcos[:,:,:,nxs]
     1422#        RBS = RBM*psin[:,:,:,nxs]
     1423#        RAS = RAM*psin[:,:,:,nxs]
     1424#        RBC = RBM*pcos[:,:,:,nxs]
     1425#        MmodA = RAC-RBS
     1426#        MmodB = RAS+RBC
     1427#        MmodA *= thdetR[:,nxs,nxs,nxs]
     1428#        MmodB *= thdetR[:,nxs,nxs,nxs]
     1429#
     1430#from 3812
     1431#    MmodA = 0; MmodB = 0
     1432#    if nWaves:
     1433#        modind = np.arange(nWaves)+1.
     1434#        phase = np.sum(twopi*XYZ[:,:,nxs,:]*modind[nxs,nxs,:,nxs]*modQ[nxs,nxs,nxs,:],axis=-1)
     1435#        phase = np.swapaxes(phase,0,1)  #Nops,Natm,Nwave
     1436#        MmodA = Am[nxs,:,:,:]*np.cos(phase[:,:,:,nxs])-Bm[nxs,:,:,:]*np.sin(phase[:,:,:,nxs])
     1437#        MmodB = Am[nxs,:,:,:]*np.sin(phase[:,:,:,nxs])+Bm[nxs,:,:,:]*np.cos(phase[:,:,:,nxs])
     1438#    return MmodA,MmodB    #cos & sin Nops,Natm,Nwaves,Mxyz
    14151439       
    14161440def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
  • trunk/GSASIIstrMath.py

    r3864 r3865  
    15091509    if parmDict[pfx+'isMag']:       #This part correct for making modulated mag moments on equiv atoms
    15101510       
    1511         GSdata = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
     1511        GSdata0 = np.inner(Gdata.T,np.swapaxes(SGMT,1,2))  #apply sym. ops.--> Natm,Nops,Nxyz
    15121512        if SGData['SGInv'] and not SGData['SGFixed']:   #inversion if any
    1513             GSdata = np.hstack((GSdata,-GSdata))     
    1514         GSdata = np.hstack([GSdata for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
    1515         GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
    1516         GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
    1517         GSdata = np.inner(GSdata,uAmat.T)   #--> cartesian
     1513            GSdata0 = np.hstack((GSdata0,-GSdata0))     
     1514        GSdata0 = np.hstack([GSdata0 for cen in SSCen])        #dup over cell centering - Natm,Nops,Mxyz
     1515        GSdata0 = SGData['MagMom'][nxs,:,nxs]*GSdata0   #flip vectors according to spin flip * det(opM)
     1516        GSdata0 = np.swapaxes(GSdata0,0,1)    #Nop,Natm,Mxyz
     1517        GSdata0 = np.inner(GSdata0,uAmat.T)   #--> cartesian
    15181518       
    15191519        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
    1520         MmodA,MmodB = G2mth.MagMod(mXYZ,modQ,MSSdata,SGData,SSGData)   #Re cos/Im sin,Nops,Natm,Nwaves,Mxyz
    1521         MmodA = np.inner(MmodA,uAmat.T)   #make cartesian
    1522         MmodB = np.inner(MmodB,uAmat.T)
     1520        Mmod,SMag = G2mth.MagMod(ngl,mXYZ,modQ,MSSdata,SGData,SSGData)   #Re cos/Im sin,Nops,Natm,Nwaves,Mxyz
     1521#        MmodA = np.inner(MmodA,uAmat.T)   #make cartesian
     1522#        MmodB = np.inner(MmodB,uAmat.T)
     1523
     1524#from #3812
     1525#1st try at this
     1526#        GSdata = Gdata.T[:,nxs,:]+Mmod    #Natm,ntau,Mxyz
     1527#        GSdata = np.inner(GSdata,SGMT).T  #apply sym. ops.--> Mxyz,Nops,Ntau,Natm
     1528#        if SGData['SGInv'] and not SGData['SGFixed']:
     1529#            GSdata = np.hstack((GSdata,-GSdata))       #inversion if any
     1530#        GSdata = np.hstack([GSdata for icen in range(Ncen)])        #dup over cell centering
     1531#        GSdata = SGData['MagMom'][nxs,:,nxs,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
     1532#        Kdata = np.inner(GSdata.T,uAmat).T     #Cartesian unit vectors
     1533#        SMag = np.sqrt(np.sum(Kdata**2,axis=0))         #Nops,Ntau,Natm
     1534#        Kmean = np.mean(SMag,axis=0)    #normalization --> unit vectors
     1535#        Kdata /= Kmean[nxs,nxs,:,:]      #mxyz,nops,ntau,natm
    15231536       
    15241537    FF = np.zeros(len(Tdata))
     
    15981611            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #& normalize
    15991612#for fixed moments --> m=0 reflections                       
    1600             fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
    1601             fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]   
     1613            fam0 = TMcorr[:,nxs,:,nxs]*GSdata0[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
     1614            fbm0 = TMcorr[:,nxs,:,nxs]*GSdata0[nxs,:,:,:]*sinm[:,:,:,nxs]   
    16021615                       
    16031616            famq0 = np.sum(np.sum(fam0,axis=-2),axis=-2)        #Nref,Mxyz; sum ops & atoms
     
    16071620            fbs0 = np.sum(fbmq0,axis=-1)**2-np.sum(HM.T*fbmq0,axis=-1)**2
    16081621#for modulated moments --> m != 0 reflections
    1609             M = np.array(np.abs(H[3]),dtype=np.int)-1
    1610             fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*cosm[i,:,:,nxs]-np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
    1611             fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*sinm[i,:,:,nxs]+np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
    1612                        
    1613             famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
    1614             fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
     1622#            M = np.array(np.abs(H[3]),dtype=np.int)-1
     1623#            fam = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*cosm[i,:,:,nxs]-np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*sinm[i,:,:,nxs]),0.) for i in range(mRef)])
     1624#            fbm = TMcorr[:,nxs,:,nxs]*np.array([np.where(M[i]>=0,(MmodA[:,:,M[i],:]*sinm[i,:,:,nxs]+np.sign(H[3])[i,nxs,nxs,nxs]*MmodB[:,:,M[i],:]*cosm[i,:,:,nxs]),0.) for i in range(mRef)])
     1625#                       
     1626#            famq = np.sum(np.sum(fam/2.,axis=-2),axis=-2)      #Nref,Mxyz; sum ops & atoms
     1627#            fbmq = np.sum(np.sum(fbm/2.,axis=-2),axis=-2)
     1628#           
     1629#            fas = np.sum(famq,axis=-1)**2-np.sum(HM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
     1630#            fbs = np.sum(fbmq,axis=-1)**2-np.sum(HM.T*fbmq,axis=-1)**2
     1631#from #3812
     1632            D = twopi*H.T[:,3:]*glTau[nxs,:]
     1633            mphase = phase[:,:,nxs,:]+D[:,nxs,:,nxs]
     1634            mphase = np.array([mphase+twopi*np.inner(cen,HP.T)[:,nxs,nxs,nxs] for cen in SGData['SGCen']])
     1635            mphase = np.concatenate(mphase,axis=1)    #remove extra axis; Nref,Nop,Natm
     1636            sinm = np.swapaxes(np.sin(mphase),2,3)    #--> Nref,Nop,Natm,Ntau
     1637            cosm = np.swapaxes(np.cos(mphase),2,3)                               #ditto
    16151638           
    1616             fas = np.sum(famq,axis=-1)**2-np.sum(HM.T*famq,axis=-1)**2      #mag intensity calc F^2-(e.F)^2
    1617             fbs = np.sum(fbmq,axis=-1)**2-np.sum(HM.T*fbmq,axis=-1)**2
     1639            HM = np.inner(Bmat,HP.T)                             #put into cartesian space
     1640            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
     1641            eDotK = np.sum(HM[:,:,nxs,nxs,nxs]*Kdata[:,nxs,:,:,:],axis=0)
     1642            Q = HM[:,:,nxs,nxs,nxs,nxs]*eDotK[nxs,:,:,:,:]-Kdata[:,nxs,:,:,:] #Mxyz,Nref,Nop,Ntau,Natm
     1643
     1644            fam = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*cosm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])   #Mxyz,Nref,Nop,Natm,Ntau,ReIm
     1645            fbm = (Q*TMcorr[nxs,:,nxs,:,nxs,nxs]*sinm[nxs,:,:,:,:,nxs]*SMag[nxs,nxs,:,:,:,nxs])
     1646           
     1647            fams = np.sum(np.sum(np.sum(fam,axis=-1),axis=2),axis=2)      #xyz,Nref,ntau; sum ops & atoms
     1648            fbms = np.sum(np.sum(np.sum(fbm,axis=-1),axis=2),axis=2)      #ditto
     1649           
     1650            fas = np.sum(fams*glWt[nxs,nxs,:],axis=-1)
     1651            fbs = np.sum(fbms*glWt[nxs,nxs,:],axis=-1)
    16181652                       
    16191653            refl.T[10] = np.where(H[3],fas+fbs,fas0+fbs0)
Note: See TracChangeset for help on using the changeset viewer.