# Changeset 3754

Ignore:
Timestamp:
Dec 6, 2018 11:52:54 AM (4 years ago)
Message:

begin work on incommensurate mag. structure factors
fix problem of magnetic reflection extinctions
replace analytical derivatives of magnetic moments with numerical ones using deltM = 0.0001

Location:
trunk
Files:
4 edited

Unmodified
Removed
• ## trunk/GSASIImath.py

 r3738 Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods as betaij Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods as betaij nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] Am = np.array(MSSdata[:3]).T   #atoms x waves x sin pos mods Bm = np.array(MSSdata[3:]).T   #...cos pos mods nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]] if nWaves[0]: tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl else: Umod = 1.0 Mmod = 1.0 if nWaves[3]: tauM = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Mwaves x ngl MmodA = Am[:,:,:,nxs]*np.sin(twopi*tauM[nxs,:,nxs,:]) #atoms X waves X 3 X ngl MmodB = Bm[:,:,:,nxs]*np.cos(twopi*tauM[nxs,:,nxs,:]) #ditto Mmod = np.sum(MmodA+MmodB,axis=1)                #atoms X 3 X ngl; sum waves Mmod = np.swapaxes(Mmod,1,2) else: Mmod = 1.0 return ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt): def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt): ''' H: array nRefBlk x ops X hklt Xmod: array atoms X 3 X ngl Umod: array atoms x 3x3 x ngl Mmod: array atoms X 3 X ngl         #TODO: add mag moment modulation math glTau,glWt: arrays Gauss-Lorentzian pos & wts ''' Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] Am = np.array(MSSdata[:3]).T   #...cos pos mods x waves x atoms Bm = np.array(MSSdata[3:]).T   #...cos pos mods nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1],Am.shape[1]] StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl)) UmodA = 0. UmodB = 0. return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB if nWaves[3]: tauM = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl StauM = np.ones_like(Am)[:,:,:,nxs]*np.sin(twopi*tauM)[nxs,:,nxs,:]   #also dMmodA/dAm CtauM = np.ones_like(Bm)[:,:,:,nxs]*np.cos(twopi*tauM)[nxs,:,nxs,:]   #also dMmodB/dBm MmodA = Am[:,:,:,nxs]*StauM #atoms x waves x 3x3 x ngl MmodB = Bm[:,:,:,nxs]*CtauM #ditto else: StauM = 1.0 CtauM = 1.0 MmodA = 0. MmodB = 0. return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB,[StauM,CtauM],MmodA+MmodB def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt): tauT = G2spc.getTauT(tau,sop,ssop,drxyz,modul)[-1] tauT *= icent       #invert wave on -1 #            print(tau,tauT,opr,G2spc.MT2text(sop).replace(' ',''),G2spc.SSMT2text(ssop).replace(' ','')) wave = np.zeros(3) uwave = np.zeros(6)
• ## trunk/GSASIIspc.py

 r3753 else: newMom = np.inner(Mom,M).T*nl.det(M)*SGData['SpnFlp'][NA+nC] #        print(len(SGOps),Ax[0],iAx,nC,nA,NA,SGData['SpnFlp'][NA],Mom,newMom) #        print(len(SGOps),Ax[0],iAx,nC,nA,NA,MT2text([M,T]).replace(' ',''),SGData['SpnFlp'][NA],Mom,newMom) #    print(Mom,newMom,MT2text([M,T]),) return newMom
• ## trunk/GSASIIstrIO.py

 r3752 ext = G2spc.checkMagextc([h,k,l],SGData) mul *= 2      # for powder overlap of Friedel pairs if ext:# and not useExt: if ext and not useExt: continue if 'C' in inst['Type'][0]:
• ## trunk/GSASIIstrMath.py

 r3595 iBeg += blkSize #    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize) return copy.deepcopy(refDict['RefList']) def MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): '''Compute magnetic structure factor derivatives numerically on blocks of reflections - for powders/nontwins only input: :param dict refDict: where 'RefList' list where each ref = h,k,l,it,d,... 'FF' dict of form factors - filled in below :param np.array G:      reciprocal metric tensor :param str hfx:    histogram id string :param str pfx:    phase id string :param dict SGData: space group info. dictionary output from SpcGroup :param dict calcControls: :param dict parmDict: :returns: dict dFdvDict: dictionary of magnetic derivatives ''' trefDict = copy.deepcopy(refDict) dM = 0.0001 dFdvDict = {} for parm in parmDict: if 'AM' in parm: parmDict[parm] += dM prefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict) parmDict[parm] -= 2*dM mrefList = MagStructureFactor2(trefDict,G,hfx,pfx,SGData,calcControls,parmDict) parmDict[parm] += dM dFdvDict[parm] = (prefList[:,9]-mrefList[:,9])/(2.*dM) return dFdvDict def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict): '''Compute magnetic structure factor derivatives on blocks of reflections - for powders/nontwins only '''Compute nonmagnetic structure factor derivatives on blocks of reflections in magnetic structures - for powders/nontwins only input: mSize = len(Mdata) Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T]) dMdm = np.inner(Gdata.T,Ginv).T/Mag Gones = np.ones_like(Gdata) Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops. VGi = np.sqrt(nl.det(Ginv)) Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag       #make unit vectors in Cartesian space dkdG = (np.inner(Gones.T,uAmat).T*VGi)/Mag dkdm = dkdG-Kdata*dMdm[:,nxs,:]/Mag[nxs,:,:] dFdMx = np.zeros((nRef,mSize,3)) Uij = np.array(G2lat.U6toUij(Uijdata)) bij = Mast*Uij.T dFdfr = np.zeros((nRef,mSize)) dFdx = np.zeros((nRef,mSize,3)) dFdMx = np.zeros((3,nRef,mSize)) dFdui = np.zeros((nRef,mSize)) dFdua = np.zeros((nRef,mSize,6)) eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0) Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK dqdk = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T     #Mxyz**2,Nref dqdm = dqdk[:,:,:,nxs,nxs]*dkdm[:,nxs,nxs,:,:]   #Mxyz**2,Nref,Nops,Natms dmx = Q*dMdm[:,nxs,nxs,:] dmx = dmx[nxs,:,:,:,:]+dqdm*Mag[nxs,nxs,nxs,:,:] dmx /= 2. fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm dfadfr = np.sum(fam/occ,axis=2)        #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2)          #deriv OK dfadmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*cosm[nxs,nxs,:,:,:],axis=3) dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms)  deriv OK dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2)            #deriv OK dfbdfr = np.sum(fbm/occ,axis=2)        #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2) dfbdmx = np.sum(dmx*TMcorr[nxs,nxs,:,nxs,:]*sinm[nxs,nxs,:,:,:],axis=3) dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms) dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2) dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/Nops,axis=0) #ok dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0)         #ok dFdMx[:,iBeg:iFin,:] = 2.*np.sum(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,axis=0)          #problems dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0)              #ok dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0)      #ok dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] dFdvDict[pfx+'AMx:'+str(i)] = dFdMx[0,:,i] dFdvDict[pfx+'AMy:'+str(i)] = dFdMx[1,:,i] dFdvDict[pfx+'AMz:'+str(i)] = dFdMx[2,:,i] dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x refBlk x sym X atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x refBlk x sym X atoms fag = fa*GfpuA[0]-fb*GfpuA[1]   #real; 2 x refBlk x sym x atoms fbg = fb*GfpuA[0]+fa*GfpuA[1] waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5] modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) FF = np.zeros(len(Tdata)) fot = (FF+FP-Bab)*Tcorr     #ops x atoms fotp = FPP*Tcorr            #ops x atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) # GfpuA is 2 x ops x atoms waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5] modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) FF = np.zeros(len(Tdata)) fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr     #ops x atoms fotp = FPP*Tcorr            #ops x atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) # GfpuA is 2 x ops x atoms waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) ngl,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast) waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,MSSdata,Mast)[:5] modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) FF = np.zeros(len(Tdata)) fot = (FF+FP-Bab)*Tcorr     #ops x atoms fotp = FPP*Tcorr            #ops x atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt) #2 x sym X atoms dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) # GfpuA is 2 x ops x atoms dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) else: if Phase['General']['Type'] == 'magnetic': if Phase['General']['Type'] == 'magnetic': dFdvDict = MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict) dFdvDict.update(MagStructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)) else: dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.