Changeset 3867 for trunk/GSASIImath.py


Ignore:
Timestamp:
Apr 2, 2019 3:16:24 PM (3 years ago)
Author:
vondreele
Message:

modulated mag moments now ok; not sf calc

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3865 r3867  
    13621362    '''
    13631363    Bm = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods
    1364     Am = np.array(MSSdata[3:]).T   #...cos pos mods
     1364    Am = np.array(MSSdata[3:]).T  #...cos pos mods
    13651365    nWaves = Am.shape[1]
    1366     ngl = 20
    13671366    tau = np.arange(ngl)/ngl
    13681367    if not nWaves:
    13691368        return 0.0,0.0
    13701369    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
     1370    Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
    13711371    SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    1372     RI = np.array([ops[0][3,3] for ops in SSGData['SSGOps']])
    13731372    if SGData['SGInv']:
    13741373        SGMT = np.vstack((SGMT,-SGMT))
     1374        Sinv =np.vstack((Sinv,-Sinv))
    13751375        SGT = np.vstack((SGT,-SGT))
    1376         RI = np.concatenate((RI,-RI))
    13771376    SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
    1378     thdetR = np.array([nl.det(op) for op in SGMT])*SGData['SpnFlp']
     1377    Sinv = np.vstack([Sinv for cen in SGData['SGCen']])
    13791378    SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
    1380     RI = np.hstack([RI for cen in SSGData['SSGCen']])
     1379    detSM = nl.det(SGMT)
     1380    mst = Sinv[:,3,:3]
     1381    epsinv = Sinv[:,3,3]
     1382    phi = np.inner(XYZ,modQ).T
     1383    TA = np.sum(mst[nxs,:,:]*(XYZ-SGT[:,:3][nxs,:,:]),axis=-1).T
     1384    tauT =  TA[nxs,:,:] + epsinv[nxs,:,nxs]*(tau[:,nxs,nxs]-SGT[:,3][nxs,:,nxs]+phi[nxs,:,:])
    13811385    modind = np.arange(nWaves)+1.
    1382     modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = -1.
    1383     toff = RI*np.inner(modQp,SGT)
    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
     1386    phase = (modind[:,nxs,nxs]*tauT)     #Nops,Natm,Nwave
    13871387    psin = np.sin(twopi*phase)
    13881388    pcos = np.cos(twopi*phase)
    1389     Mmod = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,:,nxs]+Bm[nxs,nxs,:,:,:]*psin[:,:,:,:,nxs],axis=3)
     1389    Mmod = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs]+Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)
     1390    if SGData['SGGray']:
     1391        Mmod = -np.sum(SGMT[nxs,:,nxs,:,:]*Mmod[:,:,:,nxs,:],axis=-1)*detSM[nxs,:,nxs,nxs]
     1392    else:
     1393        Mmod = np.sum(SGMT[nxs,:,nxs,:,:]*Mmod[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
    13901394    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
    14391395       
    14401396def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
Note: See TracChangeset for help on using the changeset viewer.