Changeset 3867 for trunk/GSASIImath.py
 Timestamp:
 Apr 2, 2019 3:16:24 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r3865 r3867 1362 1362 ''' 1363 1363 Bm = np.array(MSSdata[:3]).T #atoms x waves x sin pos mods 1364 Am = np.array(MSSdata[3:]).T 1364 Am = np.array(MSSdata[3:]).T #...cos pos mods 1365 1365 nWaves = Am.shape[1] 1366 ngl = 201367 1366 tau = np.arange(ngl)/ngl 1368 1367 if not nWaves: 1369 1368 return 0.0,0.0 1370 1369 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']]) 1371 1371 SGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1372 RI = np.array([ops[0][3,3] for ops in SSGData['SSGOps']])1373 1372 if SGData['SGInv']: 1374 1373 SGMT = np.vstack((SGMT,SGMT)) 1374 Sinv =np.vstack((Sinv,Sinv)) 1375 1375 SGT = np.vstack((SGT,SGT)) 1376 RI = np.concatenate((RI,RI))1377 1376 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']]) 1379 1378 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,:,:]*(XYZSGT[:,:3][nxs,:,:]),axis=1).T 1384 tauT = TA[nxs,:,:] + epsinv[nxs,:,nxs]*(tau[:,nxs,nxs]SGT[:,3][nxs,:,nxs]+phi[nxs,:,:]) 1381 1385 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 1387 1387 psin = np.sin(twopi*phase) 1388 1388 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] 1390 1394 return Mmod #Ntau,Nops,Natm,,Mxyz 1391 1392 1393 1394 ##works for DyMn6Ge6 but not MnWO41395 # nCen = len(SSGData['SSGCen'])1396 # nEqv = XYZ.shape[1] #full no. of equivalent pos incl centering1397 # mEqv = nEqv//nCen #no. operators not with centering; includes inversion1398 # MmodAA = 0; MmodBA = 01399 # if nWaves:1400 # phaseA = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ))).T #= 2pimk.r Nops,Natm,Nwave1401 # if nCen > 0:1402 # phshp = phaseA.shape1403 # 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 ##fails1413 # modQp = np.zeros(4); modQp[:3] = modQ; modQp[3] = 1.1414 # toff = RI*np.inner(modQp,SGT)1415 # MmodA = 0; MmodB = 01416 # if nWaves:1417 # modind = np.arange(nWaves)+1.1418 # phase = twopi*(modind[:,nxs,nxs]*(np.inner(XYZ,modQ)toff)).T1419 # 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 = RACRBS1426 # MmodB = RAS+RBC1427 # MmodA *= thdetR[:,nxs,nxs,nxs]1428 # MmodB *= thdetR[:,nxs,nxs,nxs]1429 #1430 #from 38121431 # MmodA = 0; MmodB = 01432 # 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,Nwave1436 # 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,Mxyz1439 1395 1440 1396 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
Note: See TracChangeset
for help on using the changeset viewer.