# Changeset 3754 for trunk/GSASIIstrMath.py

Ignore:
Timestamp:
Dec 6, 2018 11:52:54 AM (3 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

File:
1 edited

Unmodified
Added
Removed
• ## 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.