Changeset 3865 for trunk/GSASIImath.py


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

a new attempt at incommensurate mag. str. fctrs.

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