This notebook contains an mPDF analysis of the DMS data (all three compositions) taken on NOMAD in May 2015.
In [1]:
### Import a bunch of stuff.
import scipy
from scipy import interpolate
from scipy.optimize.minpack import curve_fit
import numpy as np
import matplotlib.pyplot as plt
import sys
import random
import lmfit ## important note: I had to install ipywidgets to get lmfit to work inside the ipython environment
from lmfit import Parameters, minimize, fit_report
from mcalculator import *
%matplotlib notebook
In [2]:
### Define lots of modules
def fitfunc(x,lowerb,upperb,scale1,scale2,theta,phi,width,damp,rSr,Sr,para):
'''
x: meaningless array, can be simply np.array([0]), just needed to make the curve_fit module work
scale1: scale factor of the correlated part of d(r) (coming from the ideal mPDF)
scale2: scale factor of the "paramagnetic" part of d(r)
width: smoothing factor when calculating the mPDF
theta, phi: angles giving direction of "up"-spin in the cubic coordinate system
damp: full-width half max of overall gaussian envelope applied to mPDF.
'''
svec = np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
svec=S*svec
spins[:]=svec
[r,fr]=calculateMPDF(atoms,spins,uclist,rstep,rcalcmax,width) ### ideal f(r)
rDr,Dr = cv(r,fr,rSr,Sr) ### correlated term in d(r)
Dr = scale1*Dr
Dr[:len(para)] += scale2*para ### adding paramagnetic term
th=(np.sin(qmax*r)-np.sin(qmin*r))/np.pi/r ### convolution to simulate effects of finite qmin and qmax
rDrcv, Drcv = cv(rDr,Dr,r,th)
dampf=np.exp(-(rDrcv/2.0/damp)**2)
Drcv=dampf*Drcv
#rDrcv, Drcv = rDr,Dr
return Drcv[np.logical_and(rDrcv>lowerb+0.5*rstep,rDrcv<=upperb+0.5*rstep)]
def residual(pars, x, data=None):
vals = pars.valuesdict()
scalePara=vals['scalePara']
scaleCorr=vals['scaleCorr']
width=vals['width']
damp=vals['damp']
theta=vals['theta']
phi=vals['phi']
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
model = fitfunc(x,rmin,rmax,scaleCorr,scalePara,theta,phi,width,damp,rSr,Sr,para)
if data is None:
return model
return (model - data)
def PDFdiffGrabber(inName,outName,scale=1.0):
'''
inName = file name of the .fgr exported fit file from PDFgui
outName = desired file name for the experimental mPDF (just the difference curve from the calculated structural PDF)
scale = float giving refined scale factor from PDFgui for normalization purposes
'''
allcols = np.loadtxt(inName,unpack=True,comments='#',skiprows=14)
r,diff=allcols[0],allcols[4]
np.savetxt(outName,np.transpose((r,diff/scale)))
def getRlat(a,b,c):
return 2*np.pi*np.cross(b,c)/(np.dot(a,np.cross(b,c))),2*np.pi*np.cross(c,a)/(np.dot(b,np.cross(c,a))),2*np.pi*np.cross(a,b)/(np.dot(c,np.cross(a,b)))
First, I will do the mPDF refinements for the x=0.2,y=0.15 data set.
In [3]:
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_'
inSuffix='K.fgr'
outPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','050','100','130','150','200','250','300']
temps,scaleFactors=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
for i in range(len(tempStrings)):
PDFdiffGrabber(inPrefix+tempStrings[i]+inSuffix,outPrefix+tempStrings[i]+outSuffix,scale=scaleFactors[i])
In [4]:
### Set some important values for calculation and structure
rstep=0.01
rmin=0.5 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
In [5]:
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
for i in [0,1,2,3,4,5,6,7]:
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0,0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
spincell=S*np.array([svec,svec,svec,svec])
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
spins=np.empty([len(latos)*len(spincell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
spins[index] = spincell[j]
index+=1
### Perform the calculations.
uclist=np.array([0]) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
print 'Refining...'
fit_params = Parameters()
fit_params.add('scalePara', value=np.random.uniform(0,1.0),min=0,max=20.0)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('scaleCorr', value=np.random.uniform(0,1.0),min=0.05,max=2.5)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('width', value=0.2,vary=False)
fit_params.add('damp', value=20,min=0,max=100)
fit_params.add('theta',value=np.arccos(np.random.uniform(-1,1)),min=0.000001,max=np.pi)#value=theta0,vary=False)# leave spins fixed for now
fit_params.add('phi',value=np.random.uniform(-1.0*np.pi, np.pi),min=-1.0*np.pi,max=np.pi)#value=phi0,vary=False)
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
out = minimize(residual, fit_params, args=(x,), kws={'data':data})
fit = residual(fit_params, x)
print 'Results for '+tempStrings[i]+'K:'
print fit_report(fit_params)
refinedVals=fit_params.valuesdict()
scaleCorr=refinedVals['scaleCorr']
scalePara=refinedVals['scalePara']
width=refinedVals['width']
damp=refinedVals['damp']
theta=refinedVals['theta']
phi=refinedVals['phi']
scaleCorrs[i]=scaleCorr
scaleParas[i]=scalePara
damps[i]=damp
widths[i]=width
thetas[i]=theta
phis[i]=phi
finalSpin=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
print 'Final refined spin: '+str(finalSpin)+'\n'
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
diff = data-fit
chisq=np.sum((diff)**2/len(diff))
print chisq
chisqs[i]=chisq
spins[:]=S*finalSpin
[r,fr]=calculateMPDF(atoms,spins,uclist,rstep,rcalcmax,width) ### ideal f(r)
rfull=expr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
datafull=expDr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
compexpDr=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
offset = 1.25*np.abs(np.min(data))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(rfull,datafull,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,fit,'r-',lw=2)
ax.plot(rcomp,np.zeros_like(rcomp)-offset,'k-',rcomp,diff-offset,'g-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
plt.show()
###### provide options to save data
fitstring=' Experimental data: '+dataFileName+\
'\n Chi-squared: '+str(chisq)+\
'\n Correlated scale: '+str(scaleCorr)+\
'\n Paramagnetic scale: '+str(scalePara)+\
'\n Broadening factor: '+str(width)+\
'\n Theta: '+str(theta)+\
'\n Phi: '+str(phi)+\
'\n Damp: '+str(damp)+\
'\n '+\
'\n Column format: r, Obs., Calc., Diff.'
savefile='x02y015/mPDF/mPDFfit_I4mmm_0-20_ferroFree_'+tempStrings[i]+'K_1.txt'
np.savetxt(savefile,np.column_stack((rcomp,data,fit,diff)),header=fitstring)
In [9]:
plt.close('all')
#np.savetxt('scale_correlated.txt',np.transpose((temps,scaleCorrs)))
In [10]:
### Examine the temperature-dependent results a little more closely
goodidxs=np.array([0,1,2,4,5,7])
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps[goodidxs],scaleParas[goodidxs],'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Paramagnetic Scale')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps[goodidxs],scaleCorrs[goodidxs],'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Scale Factor')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleCorrs/scaleParas,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Scale Ratio')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,damps,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Correlation length')
Sxs=np.sin(thetas)*np.cos(phis)
Sys=np.sin(thetas)*np.sin(phis)
Szs=np.cos(thetas)
for i in range(len(Szs)):
if Szs[i] < 0:
Sxs[i]=-1.0*Sxs[i]
Sys[i]=-1.0*Sys[i]
Szs[i]=-1.0*Szs[i]
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Szs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_z$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sxs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_x$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sys,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_y$')
Out[10]:
In [ ]:
plt.close('all')
Now I'll do the same analysis for the x=0.0, y=0.15 sample.
In [ ]:
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x00y015/mPDF/exportedPDFfits/I4mmm_0-20_'
inSuffix='K.fgr'
outPrefix='x00y015/mPDF/exportedPDFfits/I4mmm_0-20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','030','150','300']
temps,scaleFactors=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
for i in range(len(tempStrings)):
PDFdiffGrabber(inPrefix+tempStrings[i]+inSuffix,outPrefix+tempStrings[i]+outSuffix,scale=scaleFactors[i])
In [ ]:
### Set some important values for calculation and structure
rstep=0.01
rmin=0.2 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
In [ ]:
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
for i in range(len(temps)):
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0,0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
spincell=S*np.array([svec,svec,svec,svec])
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
spins=np.empty([len(latos)*len(spincell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
spins[index] = spincell[j]
index+=1
### Perform the calculations.
uclist=np.array([0]) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
print 'Refining...'
fit_params = Parameters()
fit_params.add('scalePara', value=np.random.uniform(0,1.0),min=0,max=20.0)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('scaleCorr', value=np.random.uniform(0,1.0),min=0.05,max=2.5)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('width', value=0.2,vary=False)
fit_params.add('damp', value=5,min=0,max=100)
fit_params.add('theta',value=np.arccos(np.random.uniform(-1,1)),min=0.000001,max=np.pi)#value=theta0,vary=False)# leave spins fixed for now
fit_params.add('phi',value=np.random.uniform(-1.0*np.pi, np.pi),min=-1.0*np.pi,max=np.pi)#value=phi0,vary=False)
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
out = minimize(residual, fit_params, args=(x,), kws={'data':data})
fit = residual(fit_params, x)
print 'Results for '+tempStrings[i]+'K:'
print fit_report(fit_params)
refinedVals=fit_params.valuesdict()
scaleCorr=refinedVals['scaleCorr']
scalePara=refinedVals['scalePara']
width=refinedVals['width']
damp=refinedVals['damp']
theta=refinedVals['theta']
phi=refinedVals['phi']
scaleCorrs[i]=scaleCorr
scaleParas[i]=scalePara
damps[i]=damp
widths[i]=width
thetas[i]=theta
phis[i]=phi
finalSpin=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
print 'Final refined spin: '+str(finalSpin)+'\n'
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
diff = data-fit
chisq=np.sum((diff)**2/len(diff))
print chisq
chisqs[i]=chisq
spins[:]=S*finalSpin
[r,fr]=calculateMPDF(atoms,spins,uclist,rstep,rcalcmax,width) ### ideal f(r)
rfull=expr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
datafull=expDr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
compexpDr=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
offset = 1.25*np.abs(np.min(data))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(rfull,datafull,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,fit,'r-',lw=2)
ax.plot(rcomp,np.zeros_like(rcomp)-offset,'k-',rcomp,diff-offset,'g-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
#plt.show()
###### provide options to save data
fitstring=' Experimental data: '+dataFileName+\
'\n Chi-squared: '+str(chisq)+\
'\n Correlated scale: '+str(scaleCorr)+\
'\n Paramagnetic scale: '+str(scalePara)+\
'\n Broadening factor: '+str(width)+\
'\n Theta: '+str(theta)+\
'\n Phi: '+str(phi)+\
'\n Damp: '+str(damp)+\
'\n '+\
'\n Column format: r, Obs., Calc., Diff.'
savefile='x00y015/mPDF/mPDFfit_I4mmm_0-20_ferroFree_'+tempStrings[i]+'K_0.txt'
np.savetxt(savefile,np.column_stack((rcomp,data,fit,diff)),header=fitstring)
In [ ]:
plt.close('all')
In [ ]:
### Examine the temperature-dependent results a little more closely
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleParas,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Paramagnetic Scale')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleCorrs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Correlated Scale')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleCorrs/scaleParas,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Scale Ratio')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,damps,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Correlation length')
Sxs=np.sin(thetas)*np.cos(phis)
Sys=np.sin(thetas)*np.sin(phis)
Szs=np.cos(thetas)
for i in range(len(Szs)):
if Szs[i] < 0:
Sxs[i]=-1.0*Sxs[i]
Sys[i]=-1.0*Sys[i]
Szs[i]=-1.0*Szs[i]
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Szs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_z$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sxs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_x$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sys,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_y$')
In [ ]:
plt.close('all')
Finally, we do the analysis for the x=0.4, y=0.15 sample.
In [ ]:
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x04y015/mPDF/exportedPDFfits/I4mmm_0-20_'
inSuffix='K.fgr'
outPrefix='x04y015/mPDF/exportedPDFfits/I4mmm_0-20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','100','230']
temps,scaleFactors=np.loadtxt('x04y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
for i in range(len(tempStrings)):
PDFdiffGrabber(inPrefix+tempStrings[i]+inSuffix,outPrefix+tempStrings[i]+outSuffix,scale=scaleFactors[i])
In [ ]:
### Set some important values for calculation and structure
rstep=0.01
rmin=0.7 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
In [ ]:
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x04y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x04y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
for i in range(len(temps)):
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0,0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
spincell=S*np.array([svec,svec,svec,svec])
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
spins=np.empty([len(latos)*len(spincell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
spins[index] = spincell[j]
index+=1
### Perform the calculations.
uclist=np.array([0]) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
print 'Refining...'
fit_params = Parameters()
fit_params.add('scalePara', value=np.random.uniform(0,1.0),min=0,max=20.0)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('scaleCorr', value=np.random.uniform(0,1.0),min=0.05,max=2.5)#value=np.random.uniform(5,25),min=0,max=30)
fit_params.add('width', value=0.2,vary=False)
fit_params.add('damp', value=5,min=0,max=100)
fit_params.add('theta',value=np.arccos(np.random.uniform(-1,1)),min=0.000001,max=np.pi)#value=theta0,vary=False)# leave spins fixed for now
fit_params.add('phi',value=np.random.uniform(-1.0*np.pi, np.pi),min=-1.0*np.pi,max=np.pi)#value=phi0,vary=False)
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
out = minimize(residual, fit_params, args=(x,), kws={'data':data})
fit = residual(fit_params, x)
print 'Results for '+tempStrings[i]+'K:'
print fit_report(fit_params)
refinedVals=fit_params.valuesdict()
scaleCorr=refinedVals['scaleCorr']
scalePara=refinedVals['scalePara']
width=refinedVals['width']
damp=refinedVals['damp']
theta=refinedVals['theta']
phi=refinedVals['phi']
scaleCorrs[i]=scaleCorr
scaleParas[i]=scalePara
damps[i]=damp
widths[i]=width
thetas[i]=theta
phis[i]=phi
finalSpin=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
print 'Final refined spin: '+str(finalSpin)+'\n'
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
diff = data-fit
chisq=np.sum((diff)**2/len(diff))
print chisq
chisqs[i]=chisq
spins[:]=S*finalSpin
[r,fr]=calculateMPDF(atoms,spins,uclist,rstep,rcalcmax,width) ### ideal f(r)
rfull=expr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
datafull=expDr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
compexpDr=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
offset = 1.25*np.abs(np.min(data))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(rfull,datafull,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,fit,'r-',lw=2)
ax.plot(rcomp,np.zeros_like(rcomp)-offset,'k-',rcomp,diff-offset,'g-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
#plt.show()
###### provide options to save data
fitstring=' Experimental data: '+dataFileName+\
'\n Chi-squared: '+str(chisq)+\
'\n Correlated scale: '+str(scaleCorr)+\
'\n Paramagnetic scale: '+str(scalePara)+\
'\n Broadening factor: '+str(width)+\
'\n Theta: '+str(theta)+\
'\n Phi: '+str(phi)+\
'\n Damp: '+str(damp)+\
'\n '+\
'\n Column format: r, Obs., Calc., Diff.'
savefile='x04y015/mPDF/mPDFfit_I4mmm_0-20_ferroFree_'+tempStrings[i]+'K_0.txt'
np.savetxt(savefile,np.column_stack((rcomp,data,fit,diff)),header=fitstring)
In [11]:
plt.close('all')
In [ ]:
### Examine the temperature-dependent results a little more closely
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleParas,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Paramagnetic Scale')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleCorrs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Correlated Scale')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,scaleCorrs/scaleParas,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Scale Ratio')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,damps,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Correlation length')
Sxs=np.sin(thetas)*np.cos(phis)
Sys=np.sin(thetas)*np.sin(phis)
Szs=np.cos(thetas)
for i in range(len(Szs)):
if Szs[i] < 0:
Sxs[i]=-1.0*Sxs[i]
Sys[i]=-1.0*Sys[i]
Szs[i]=-1.0*Szs[i]
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Szs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_z$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sxs,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_x$')
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(temps,Sys,'bo')
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('S$_y$')
In [ ]:
plt.close('all')
In [3]:
### Define lots of modules
def shellcalc(x,lowerb,upperb,scaleNN1,scaleNN2,scaleNN3,scaleNN4,scaleNN5,scaleNN6,scalePara,width,theta,phi,damp,dampNN1,dampNN2,dampNN3,dampNN4,dampNN5,dampNN6):
'''
x: meaningless array, can be simply np.array([0]), just needed to make the curve_fit module work
scale1: scale factor of the correlated part of d(r) (coming from the ideal mPDF)
scale2: scale factor of the "paramagnetic" part of d(r)
width: smoothing factor when calculating the mPDF
theta, phi: angles giving direction of "up"-spin in the cubic coordinate system
damp: full-width half max of overall gaussian envelope applied to mPDF.
structureScale: option factor by which to uniformly scale the positions of all the atoms.
'''
svec = np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
svec=S*svec
spinsNN1[:]=svec
spinsNN2[:]=svec
spinsNN3[:]=svec
spinsNN4[:]=svec
spinsNN5[:]=svec
spinsNN6[:]=svec
[r,frNN1]=calculateMPDF(atomsNN1,spinsNN1,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN1 = cv(r,frNN1,rSr,Sr) ### correlated term in d(r)
DrNN1 = scaleNN1*DrNN1*np.exp(-1.0*rDr/dampNN1)
frNN1=scaleNN1*frNN1*np.exp(-1.0*r/dampNN1)
[r,frNN2]=calculateMPDF(atomsNN2,spinsNN2,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN2 = cv(r,frNN2,rSr,Sr) ### correlated term in d(r)
DrNN2 = scaleNN2*DrNN2*np.exp(-1.0*rDr/dampNN2)
frNN2 = scaleNN2*frNN2*np.exp(-1.0*r/dampNN2)
[r,frNN3]=calculateMPDF(atomsNN3,spinsNN3,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN3 = cv(r,frNN3,rSr,Sr) ### correlated term in d(r)
DrNN3 = scaleNN3*DrNN3*np.exp(-1.0*rDr/dampNN3)
frNN3 = scaleNN3*frNN3*np.exp(-1.0*r/dampNN3)
[r,frNN4]=calculateMPDF(atomsNN4,spinsNN4,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN4 = cv(r,frNN4,rSr,Sr) ### correlated term in d(r)
DrNN4 = scaleNN4*DrNN4*np.exp(-1.0*rDr/dampNN4)
frNN4 = scaleNN4*frNN4*np.exp(-1.0*r/dampNN4)
[r,frNN5]=calculateMPDF(atomsNN5,spinsNN5,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN5 = cv(r,frNN5,rSr,Sr) ### correlated term in d(r)
DrNN5 = scaleNN5*DrNN5*np.exp(-1.0*rDr/dampNN5)
frNN5 = scaleNN5*frNN5*np.exp(-1.0*r/dampNN5)
[r,frNN6]=calculateMPDF(atomsNN6,spinsNN6,np.array([0]),rstep,rcalcmax,width) ### ideal f(r)
rDr,DrNN6 = cv(r,frNN6,rSr,Sr) ### correlated term in d(r)
DrNN6 = scaleNN6*DrNN6*np.exp(-1.0*rDr/dampNN6)
frNN6 = scaleNN6*frNN6*np.exp(-1.0*r/dampNN6)
Dr = DrNN1+DrNN2+DrNN3+DrNN4+DrNN5+DrNN6
fr = frNN1+frNN2+frNN3+frNN4+frNN5+frNN6
Dr[:len(para)] += scalePara*para ### adding paramagnetic term
th=(np.sin(qmax*r)-np.sin(qmin*r))/np.pi/r ### convolution to simulate effects of finite qmin and qmax
rDrcv, Drcv = cv(rDr,Dr,r,th)
rfrcv, frcv = cv(r,fr,r,th)
dampD=np.exp(-1.0*rDrcv/damp)
dampf=np.exp(-1.0*rfrcv/damp)
Drcv=dampD*Drcv
frcv=dampf*frcv
#rDrcv, Drcv = rDr,Dr
return Drcv[np.logical_and(rDrcv>lowerb+0.5*rstep,rDrcv<=upperb+0.5*rstep)],rfrcv[np.logical_and(rfrcv>lowerb+0.5*rstep,rfrcv<=upperb+0.5*rstep)],frcv[np.logical_and(rfrcv>lowerb+0.5*rstep,rfrcv<=upperb+0.5*rstep)]
def shellResidual(pars, x, data=None):
vals = pars.valuesdict()
scaleNN1=vals['scaleNN1']
scaleNN2=vals['scaleNN2']
scaleNN3=vals['scaleNN3']
scaleNN4=vals['scaleNN4']
scaleNN5=vals['scaleNN5']
scaleNN6=vals['scaleNN6']
dampNN1=vals['dampNN1']
dampNN2=vals['dampNN2']
dampNN3=vals['dampNN3']
dampNN4=vals['dampNN4']
dampNN5=vals['dampNN5']
dampNN6=vals['dampNN6']
scalePara=vals['scalePara']
damp=vals['damp']
theta=vals['theta']
phi=vals['phi']
width=vals['width']
model = shellcalc(x,rmin,rmax,scaleNN1,scaleNN2,scaleNN3,scaleNN4,scaleNN5,scaleNN6,scalePara,width,theta,phi,damp,dampNN1,dampNN2,dampNN3,dampNN4,dampNN5,dampNN6)[0]
if data is None:
return model
return (model - data)
def neighborFinder(shell,atoms,spins):
idxs=np.abs(datomsmag-shellvals[shell])<0.0025*shellvals[shell]
goodAs=atoms[idxs]
goodAs=np.concatenate((np.array([atoms[0]]),goodAs))
goodSs=spins[idxs]
goodSs=np.concatenate((np.array([spins[0]]),goodSs))
return goodAs,goodSs
In [4]:
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_'
inSuffix='K.fgr'
outPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','050','100','130','150','200','250','300']
temps,scaleFactors=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
### Set some important values for calculation and structure
rstep=0.01
rmin=0.5 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
for i in range(len(temps)):
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0,0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
spincell=S*np.array([svec,svec,svec,svec])
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
spins=np.empty([len(latos)*len(spincell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
spins[index] = spincell[j]
index+=1
datoms = atoms-atoms[0]
datomsmag = np.sqrt(np.sum((datoms)**2,axis=1).reshape(datoms.shape[0],1)).ravel()
datomsmag = np.around(datomsmag,decimals=5)
shellvals = np.unique(datomsmag)
### Extract neighbor shells
atomsNN1,spinsNN1=neighborFinder(1,atoms,spins)
atomsNN2,spinsNN2=neighborFinder(2,atoms,spins)
atomsNN3,spinsNN3=neighborFinder(3,atoms,spins)
atomsNN4,spinsNN4=neighborFinder(4,atoms,spins)
atomsNN5,spinsNN5=neighborFinder(5,atoms,spins)
atomsNN6,spinsNN6=neighborFinder(6,atoms,spins)
uclist=np.array([0]) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
print 'Refining...'
fit_params = Parameters()
fit_params.add('scalePara', value=0.45,min=0,max=20.0)
fit_params.add('scaleNN1', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN2', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN3', value=np.random.uniform(0,1.0),min=0.0,max=3)
fit_params.add('scaleNN4', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN5', value=0,vary=False)#np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN6', value=0,vary=False)#np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('dampNN1',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN2',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN3',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN4',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN5',value=1000,vary=False)#np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN6',value=1000,vary=False)#np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('width', value=0.2,vary=False)
fit_params.add('damp', value=100,vary=False)#min=0.1,max=100)#vary=False)
fit_params.add('theta',value=theta0,vary=False)#np.arccos(np.random.uniform(-1,1)),min=0.000001,max=np.pi)
fit_params.add('phi',value=theta0,vary=False)#np.random.uniform(-1.0*np.pi, np.pi),min=-1.0*np.pi,max=np.pi)
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
out = minimize(shellResidual, fit_params, args=(x,), kws={'data':data})
fit = shellResidual(fit_params, x)
print 'Results for '+tempStrings[i]+'K:'
print fit_report(fit_params)
refinedVals=fit_params.valuesdict()
scaleNN1=refinedVals['scaleNN1']
scaleNN2=refinedVals['scaleNN2']
scaleNN3=refinedVals['scaleNN3']
scaleNN4=refinedVals['scaleNN4']
scaleNN5=refinedVals['scaleNN5']
scaleNN6=refinedVals['scaleNN6']
scalePara=refinedVals['scalePara']
width=refinedVals['width']
damp=refinedVals['damp']
dampNN1=refinedVals['dampNN1']
dampNN2=refinedVals['dampNN2']
dampNN3=refinedVals['dampNN3']
dampNN4=refinedVals['dampNN4']
dampNN5=refinedVals['dampNN5']
dampNN6=refinedVals['dampNN6']
theta=refinedVals['theta']
phi=refinedVals['phi']
# scaleCorrs[i]=scaleCorr
# scaleParas[i]=scalePara
# damps[i]=damp
# widths[i]=width
# thetas[i]=theta
# phis[i]=phi
finalSpin=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
print 'Final refined spin: '+str(finalSpin)+'\n'
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
diff = data-fit
chisq=np.sum((diff)**2/len(diff))
print chisq
chisqs[i]=chisq
spins[:]=S*finalSpin
[r,fr]=shellcalc(x,rmin,rmax,scaleNN1,scaleNN2,scaleNN3,scaleNN4,scaleNN5,scaleNN6,scalePara,width,theta,phi,damp,dampNN1,dampNN2,dampNN3,dampNN4,dampNN5,dampNN6)[1:]
rfull=expr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
datafull=expDr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
compexpDr=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
offset = 1.25*np.abs(np.min(data))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(rfull,datafull,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,fit,'r-',lw=2)
ax.plot(rcomp,np.zeros_like(rcomp)-offset,'k-',rcomp,diff-offset,'g-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
plt.show()
###### provide options to save data
fitstring=' Experimental data: '+dataFileName+\
'\n Chi-squared: '+str(chisq)+\
'\n NN1 scale, damp: '+str(scaleNN1)+'\t'+str(dampNN1)+\
'\n NN2 scale, damp: '+str(scaleNN2)+'\t'+str(dampNN2)+\
'\n NN3 scale, damp: '+str(scaleNN3)+'\t'+str(dampNN3)+\
'\n NN4 scale, damp: '+str(scaleNN4)+'\t'+str(dampNN4)+\
'\n NN5 scale, damp: '+str(scaleNN5)+'\t'+str(dampNN5)+\
'\n NN6 scale, damp: '+str(scaleNN6)+'\t'+str(dampNN6)+\
'\n Paramagnetic scale: '+str(scalePara)+\
'\n Broadening factor: '+str(width)+\
'\n Theta: '+str(theta)+\
'\n Phi: '+str(phi)+\
'\n Damp: '+str(damp)+\
'\n '+\
'\n Column format: r, Obs., Calc., Diff.'
savefile='x02y015/mPDF/mPDFfit_I4mmm_0-20_ferroFree_shell_'+tempStrings[i]+'K_1.txt'
np.savetxt(savefile,np.column_stack((rcomp,data,fit,diff)),header=fitstring)
In [ ]:
plt.close('all')
So, the shell-by-shell SRO model does a slightly better job than the LRO model (mostly by getting the 4th NN), but only if I allow independent damping parameters for each shell, i.e. just having independent scale factors alone is not enough. This may be physically justifiable, since the contribution from every coordination shell contains a low-r part (although it disappears as 1/r_{ij}^3).
Do a double-check with the lower qmax data set.
In [5]:
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_qmax20_'
inSuffix='K.fgr'
outPrefix='x02y015/mPDF/exportedPDFfits/I4mmm_0-20_qmax20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','050','100','130','150','200','250','300']
temps,scaleFactors=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax20full_I4mmm_Tscan_scale.dat',unpack=True)
for i in range(len(tempStrings)):
PDFdiffGrabber(inPrefix+tempStrings[i]+inSuffix,outPrefix+tempStrings[i]+outSuffix,scale=scaleFactors[i])
### Set some important values for calculation and structure
rstep=0.01
rmin=0.5 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x02y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
for i in range(len(temps)):
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0,0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
spincell=S*np.array([svec,svec,svec,svec])
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
spins=np.empty([len(latos)*len(spincell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
spins[index] = spincell[j]
index+=1
datoms = atoms-atoms[0]
datomsmag = np.sqrt(np.sum((datoms)**2,axis=1).reshape(datoms.shape[0],1)).ravel()
datomsmag = np.around(datomsmag,decimals=5)
shellvals = np.unique(datomsmag)
### Extract neighbor shells
atomsNN1,spinsNN1=neighborFinder(1,atoms,spins)
atomsNN2,spinsNN2=neighborFinder(2,atoms,spins)
atomsNN3,spinsNN3=neighborFinder(3,atoms,spins)
atomsNN4,spinsNN4=neighborFinder(4,atoms,spins)
atomsNN5,spinsNN5=neighborFinder(5,atoms,spins)
atomsNN6,spinsNN6=neighborFinder(6,atoms,spins)
uclist=np.array([0]) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
print 'Refining...'
fit_params = Parameters()
fit_params.add('scalePara', value=0.45,min=0,max=20.0)
fit_params.add('scaleNN1', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN2', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN3', value=np.random.uniform(0,1.0),min=0.0,max=3)
fit_params.add('scaleNN4', value=np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN5', value=0,vary=False)#np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('scaleNN6', value=0,vary=False)#np.random.uniform(0,1.0),min=0.0,max=20)
fit_params.add('dampNN1',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN2',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN3',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN4',value=np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN5',value=1000,vary=False)#np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('dampNN6',value=1000,vary=False)#np.random.uniform(1,5),min=0.0,max=20)
fit_params.add('width', value=0.2,vary=False)
fit_params.add('damp', value=100,vary=False)#min=0.1,max=100)#vary=False)
fit_params.add('theta',value=theta0,vary=False)#np.arccos(np.random.uniform(-1,1)),min=0.000001,max=np.pi)
fit_params.add('phi',value=theta0,vary=False)#np.random.uniform(-1.0*np.pi, np.pi),min=-1.0*np.pi,max=np.pi)
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
out = minimize(shellResidual, fit_params, args=(x,), kws={'data':data})
fit = shellResidual(fit_params, x)
print 'Results for '+tempStrings[i]+'K:'
print fit_report(fit_params)
refinedVals=fit_params.valuesdict()
scaleNN1=refinedVals['scaleNN1']
scaleNN2=refinedVals['scaleNN2']
scaleNN3=refinedVals['scaleNN3']
scaleNN4=refinedVals['scaleNN4']
scaleNN5=refinedVals['scaleNN5']
scaleNN6=refinedVals['scaleNN6']
scalePara=refinedVals['scalePara']
width=refinedVals['width']
damp=refinedVals['damp']
dampNN1=refinedVals['dampNN1']
dampNN2=refinedVals['dampNN2']
dampNN3=refinedVals['dampNN3']
dampNN4=refinedVals['dampNN4']
dampNN5=refinedVals['dampNN5']
dampNN6=refinedVals['dampNN6']
theta=refinedVals['theta']
phi=refinedVals['phi']
# scaleCorrs[i]=scaleCorr
# scaleParas[i]=scalePara
# damps[i]=damp
# widths[i]=width
# thetas[i]=theta
# phis[i]=phi
finalSpin=np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
print 'Final refined spin: '+str(finalSpin)+'\n'
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
diff = data-fit
chisq=np.sum((diff)**2/len(diff))
print chisq
chisqs[i]=chisq
spins[:]=S*finalSpin
[r,fr]=shellcalc(x,rmin,rmax,scaleNN1,scaleNN2,scaleNN3,scaleNN4,scaleNN5,scaleNN6,scalePara,width,theta,phi,damp,dampNN1,dampNN2,dampNN3,dampNN4,dampNN5,dampNN6)[1:]
rfull=expr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
datafull=expDr[np.logical_and(expr>rcalcmin+0.5*rstep,expr<=rcalcmax+0.5*rstep)]
compexpDr=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
offset = 1.25*np.abs(np.min(data))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(rfull,datafull,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,fit,'r-',lw=2)
ax.plot(rcomp,np.zeros_like(rcomp)-offset,'k-',rcomp,diff-offset,'g-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
plt.show()
###### provide options to save data
fitstring=' Experimental data: '+dataFileName+\
'\n Chi-squared: '+str(chisq)+\
'\n NN1 scale, damp: '+str(scaleNN1)+'\t'+str(dampNN1)+\
'\n NN2 scale, damp: '+str(scaleNN2)+'\t'+str(dampNN2)+\
'\n NN3 scale, damp: '+str(scaleNN3)+'\t'+str(dampNN3)+\
'\n NN4 scale, damp: '+str(scaleNN4)+'\t'+str(dampNN4)+\
'\n NN5 scale, damp: '+str(scaleNN5)+'\t'+str(dampNN5)+\
'\n NN6 scale, damp: '+str(scaleNN6)+'\t'+str(dampNN6)+\
'\n Paramagnetic scale: '+str(scalePara)+\
'\n Broadening factor: '+str(width)+\
'\n Theta: '+str(theta)+\
'\n Phi: '+str(phi)+\
'\n Damp: '+str(damp)+\
'\n '+\
'\n Column format: r, Obs., Calc., Diff.'
savefile='x02y015/mPDF/mPDFfit_I4mmm_0-20_qmax20_ferroFree_shell_'+tempStrings[i]+'K_1.txt'
np.savetxt(savefile,np.column_stack((rcomp,data,fit,diff)),header=fitstring)
In [4]:
### Define lots of modules
def fitfuncAF(x,lowerb,upperb,scale1,scale2,theta,phi,width,damp,rSr,Sr,para,Q,atoms):
'''
x: meaningless array, can be simply np.array([0]), just needed to make the curve_fit module work
scale1: scale factor of the correlated part of d(r) (coming from the ideal mPDF)
scale2: scale factor of the "paramagnetic" part of d(r)
width: smoothing factor when calculating the mPDF
theta, phi: angles giving direction of "up"-spin in the cubic coordinate system
damp: full-width half max of overall gaussian envelope applied to mPDF.
'''
svec = np.array([np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)])
svec=S*svec
signs=np.cos(np.dot(atoms-atoms[0],Q))
spins=svec*signs.reshape(-1,1)
[r,fr]=calculateMPDF(atoms,spins,uclist,rstep,rcalcmax,width) ### ideal f(r)
rDr,Dr = cv(r,fr,rSr,Sr) ### correlated term in d(r)
Dr = scale1*Dr
Dr[:len(para)] += scale2*para ### adding paramagnetic term
th=(np.sin(qmax*r)-np.sin(qmin*r))/np.pi/r ### convolution to simulate effects of finite qmin and qmax
rDrcv, Drcv = cv(rDr,Dr,r,th)
dampf=np.exp(-(rDrcv/2.0/damp)**2)
Drcv=dampf*Drcv
#rDrcv, Drcv = rDr,Dr
return Drcv[np.logical_and(rDrcv>lowerb+0.5*rstep,rDrcv<=upperb+0.5*rstep)],r[np.logical_and(r>lowerb+0.5*rstep,r<=upperb+0.5*rstep)],fr[np.logical_and(r>lowerb+0.5*rstep,r<=upperb+0.5*rstep)]
def residualAF(pars, x, Q, atoms, data=None):
vals = pars.valuesdict()
scalePara=vals['scalePara']
scaleCorr=vals['scaleCorr']
width=vals['width']
damp=vals['damp']
theta=vals['theta']
phi=vals['phi']
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
model = fitfuncAF(x,rmin,rmax,scaleCorr,scalePara,theta,phi,width,damp,rSr,Sr,para,Q,atoms)[0]
if data is None:
return model
return (model - data)
In [5]:
### Do some pure calculations
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='x00y015/mPDF/exportedPDFfits/I4mmm_0-20_'
inSuffix='K.fgr'
outPrefix='x00y015/mPDF/exportedPDFfits/I4mmm_0-20_'
outSuffix='_scaleNorm.diff'
tempStrings=['002','030','150','300']
temps,scaleFactors=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
### Set some important values for calculation and structure
rstep=0.01
rmin=0.2 # rmin and rmax set the data range for the fit
rmax=18.5
rcalcmin = 0 # The calculation is extended slightly beyond the fit range to avoid artifacts at boundaries.
rcalcmax=20
qmin,qmax=0.0,31.41
[dq,q1,q2] = [0.01,0.00,10.00]
q=np.arange(q1,q2,dq)
ffMag=j0calc(q,[0.422,17.684,0.5948,6.005,0.0043,-0.609,-0.0219]) ### magnetic form factor for Ni2+
r1,r2,dr=-5.0,5.0,0.01
rsr, sr=costransform(q,ffMag,rmin=r1,rmax=r2,rstep=dr)
sr = np.sqrt(np.pi/2.0)*sr
rSr,Sr = cv(rsr,sr,rsr,sr)
para = -1.0*np.sqrt(2.0*np.pi)*np.gradient(Sr,rSr[1]-rSr[0]) ### paramagnetic term in d(r)
### Prepare the tetragonal structure for the calculations
### lattice parameters of magnetic unit cell, and some basis transformation stuff
e1=np.array([1,0,0])
e1n=e1/np.linalg.norm(e1)
e2=np.array([0,1,0])
e2n=e2/np.linalg.norm(e2)
e3=np.array([0,0,1])
e3n=e3/np.linalg.norm(e3)
temps,aVals=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_a-lat.dat',unpack=True)
temps,cVals=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_c-lat.dat',unpack=True)
Basis=np.array([e1,e2,e3])
IBasis=np.linalg.inv(Basis)
BasisN=np.array([e1n,e2n,e3n])
### positions of atoms in magnetic unit cell
basis=np.array([[0.5,0,0.75],[0,0.5,0.75],[0.5,0,0.25],[0,0.5,0.25]])
scaleParas=np.zeros_like(temps)
scaleCorrs=np.zeros_like(temps)
thetas=np.zeros_like(temps)
phis=np.zeros_like(temps)
damps=np.zeros_like(temps)
widths=np.zeros_like(temps)
chisqs=np.zeros_like(temps)
### loop through the temperatures to do the refinements
i=0
a=aVals[i]
b=aVals[i]
c=cVals[i]
latpars=np.array([[a],[b],[c]])
full=BasisN*latpars
cell=np.dot(IBasis,full)
atomcell=np.dot(basis,cell)
### spin orientations in same order as atomic positions.
svec=np.array([0.0,0.0,1])
svec=svec/np.linalg.norm(svec)
svecmag=np.sqrt((svec*svec).sum())
theta0=np.arccos(svec[2]/svecmag)
phi0=np.arctan2(svec[1],svec[0])
S=2.5
### how big to make the box
radius=1.5*rcalcmax
dim1=np.round(radius/np.linalg.norm(cell[0]))
dim2=np.round(radius/np.linalg.norm(cell[1]))
dim3=np.round(radius/np.linalg.norm(cell[2]))
### generate the coordinates of each unit cell
latos=np.dot(np.mgrid[-dim1:dim1+1,-dim2:dim2+1,-dim3:dim3+1].transpose().ravel().reshape((2*dim1+1)*(2*dim2+1)*(2*dim3+1),3),cell)
### select points within a desired radius from origin
latos=latos[np.where(np.apply_along_axis(np.linalg.norm,1,latos)<=(rcalcmax+10.0))]
## rearrange latos array so that [0,0,0] is the first one (for convenience)
latos[np.where(np.all(latos==[0,0,0],axis=1))]=latos[0]
latos[0]=np.array([0,0,0])
### create list of all atomic positions and spin directions
atoms=np.empty([len(latos)*len(atomcell),3])
index=0
for k in range(len(latos)):
for j in range(len(atomcell)):
atoms[index]=latos[k]+atomcell[j]
index+=1
### create spins using propagation vector
astar,bstar,cstar=getRlat(cell[0],cell[1],cell[2])
Q=0.0*astar+1.0*bstar+1.0*cstar
signs=np.cos(np.dot(atoms-atoms[0],Q))
spins=S*svec*signs.reshape(-1,1)
### Perform the calculations.
uclist=np.arange(1) # indices of atoms which will be used as centers in the mPDF calculation
x=np.array([0]) # empty array to be used as placeholder in the calculations
dataFileName=outPrefix+tempStrings[i]+outSuffix
expr,expDr = np.loadtxt(dataFileName,unpack=True) # grab the experimental data.
data=expDr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
rcomp = expr[np.logical_and(expr>rmin+0.5*rstep,expr<=rmax+0.5*rstep)]
print 'Calculating...'
scaleCorr,scalePara=0.075,0.35
theta,phi=theta0,phi0
width,damp=0.2,10
calcDr,rfr,fr=fitfuncAF(x,rmin,rmax,scaleCorr,scalePara,theta,phi,width,damp,rSr,Sr,para,Q,atoms)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(expr,expDr,marker='o',mfc='none',mec='b',linestyle='none')
ax.plot(rcomp,calcDr,'r-',lw=2)
#ax.plot(rfr,fr,'k-')
#ax.plot(r,fr,'b-')
ax.set_xlim(xmin=rcalcmin,xmax=rcalcmax)
ax.set_xlabel('r ($\AA$)')
ax.set_ylabel('d(r) ($\AA^{-2}$)')
Out[5]:
In [6]:
plt.close('all')
In [7]:
astar
Out[7]:
So we see that the standard AF pattern does not work at all; that first nearest neighbor peak really does seem to be ferromagnetic.
In [22]:
### Do some pure calculations
### Extract and save the experimental mPDF from the exported PDF fit file
inPrefix='SQfiles/NOM_9999_Ba122K02Mn015_'
inSuffix='K_SQ.dat'
tempStrings=['300','250','200','150','130','100','50','2']
#temps,scaleFactors=np.loadtxt('x00y015/mPDF/exportedPDFfits/qmax31lorch_I4mmm_Tscan_scale.dat',unpack=True)
qmin,qmax=0.0,10
fig=plt.figure()
ax=fig.add_subplot(111)
#q300,sq300,sqerr300=np.loadtxt(inPrefix+tempStrings[0]+inSuffix,unpack=True)
#for i in [2,7]:
# dataFileName=inPrefix+tempStrings[i]+inSuffix
# q,sq,sqerr=np.loadtxt(dataFileName,unpack=True)
# ax.plot(q,sq,label=tempStrings[i],linestyle='-')
q200,sq200,sqerr200=np.loadtxt(inPrefix+tempStrings[2]+inSuffix,unpack=True)
q2,sq2,sqerr2=np.loadtxt(inPrefix+tempStrings[7]+inSuffix,unpack=True)
ax.plot(q200,sq200,label='200 K',color='r',linestyle='-')
ax.plot(q,sq2,label='2 K',color='b',linestyle='-')
ax.plot(q,sq2-sq200-3.5,label='Low - High',color='g',linestyle='-')
ax.plot(q,np.zeros_like(q)-3.5,'k-')
plt.legend(loc=0)
ax.set_xlim(xmin=qmin,xmax=qmax)
ax.set_xlabel('Q ($\AA^{-1}$)')
ax.set_ylabel('S (arb. units)')
plt.savefig('sqdata.png',format='PNG',dpi=150,bbox_inches='tight')
In [20]:
plt.close('all')
In [ ]: