Changeset 3865 for trunk/GSASIImath.py
 Timestamp:
 Mar 31, 2019 6:59:15 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r3864 r3865 1356 1356 return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt 1357 1357 1358 def MagMod(XYZ,modQ,MSSdata,SGData,SSGData): 1359 1358 def 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 1360 1370 SGMT = np.array([ops[0] for ops in SGData['SGOps']]) #not .T!! 1361 1371 SGT = np.array([ops[1] for ops in SSGData['SSGOps']]) … … 1369 1379 SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1. 1370 1380 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. 1395 1382 modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = 1. 1396 1383 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 = RACRBS 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 = RACRBS 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 1415 1439 1416 1440 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
Note: See TracChangeset
for help on using the changeset viewer.