How to fit powder inelastic neutron scattering data
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hello, here is my code, the upper part simulates the calculation model, and the lower part fits the experimental data. I referred to Tutorial 39. The issue I encountered is whether my data format is correct so that SpinW can read it and successfully perform the fitting. Attached is my data (fitting was done only for the default Q = 2.0).
cnso = spinw;
cnso.genlattice('lat_const',[9.18620 5.30370 11.90660],'angled',[90 104.9180 90],'spgr','C 2/c')
cnso.addatom('label','Ni','r',[0.41673; 0.25003; 0.00016],'S',[1],'color','Blue');
cnso.addatom('r',[0.5 0.67001; 0.90888 0.41913; 0.25 0.25001],'S',[0 0],'label',{'Cu' 'Cu'},'color',{'lightGray' 'lightGray'});
cnso.addatom('r',[0.75; 0.75; 0.5],'S',[0],'label',{'Sb'},'color',{'Orange'});
cnso.addatom('r',[0.56040 0.72020 0.61910; 0.91050 0.42920 0.41020; 0.41100 0.41070 0.08910],'S',[0 0 0],'label',{'O' 'O' 'O'},'color',{'Red' 'Red' 'Red'});
plot(cnso,'baseShift',[0;0;0])
swplot.zoom(1.3)
% We generate all bonds up to 8 Angstrom length.
cnso.gencoupling('maxDistance',8, 'forceNoSym', true);
% Kitaev term
cnso.addmatrix('label','Jxx','value',0,'color','r');
cnso.addmatrix('label','Jyy','value',0,'color','g');
cnso.addmatrix('label','Jzz','value',0,'color','b');
% Heisenberg terms
cnso.addmatrix('label','J1','value',-1.7,'color','gray');
cnso.addmatrix('label','J2','value',0,'color','orange');
cnso.addmatrix('label','J3','value',1.8,'color','cyan');
% add J1, J2 and J3 and JK couplings
cnso.addcoupling('mat','J1','bond',[1]);
cnso.addcoupling('mat','J2','bond',[2]);
cnso.addcoupling('mat','J3','bond',[4]);
% Plot all couplings.
plot(cnso,'range',[2 2 2],'atomMode','mag','cellMode','inside','atomLegend',false,'cellcolor','gray','bondMode','line','bondLinewidth0',2)
swplot.zoom(1.4)
snapnow
% add JJxx, Jyy and Jzz couplings
cnso.addcoupling('mat','Jxx','bond',1,'subidx',[5 6 7 8]);
cnso.addcoupling('mat','Jyy','bond',1,'subidx',[2 3 10 12]);
cnso.addcoupling('mat','Jzz','bond',1,'subidx',[1 4 9 11]);
% Plot Kitaev couplings only.
plot(cnso,'range',[2 2 2],'atomMode','mag','cellMode','inside','atomLegend',false,'cellcolor','gray','bondMode','line','bondLinewidth0',2)
swplot.zoom(1.4)
%Plot Kitaev term
cnso.addmatrix('label','Jxx','value',diag([1 0 0]),'color','r')
cnso.addmatrix('label','Jyy','value',diag([0 1 0]),'color','g')
cnso.addmatrix('label','Jzz','value',diag([0 0 1]),'color','b')
cnso.addmatrix('label','J1','value',0);
cnso.addmatrix('label','J2','value', 0);
cnso.addmatrix('label','J3','value', 0);
plot(cnso,'range',[0 2;0 2;-0.1 1.1],'atomMode','mag','bondRadius1',0.15,'bondMode','line','bondLineWidth','lin','bondLinewidth0',4,'atomLegend',false)
cnso.table('atom')
%Q scans
Qp{1} = [ -1; 0; 0];
Qp{2} = [ 0; 0; 0];
Qp{3} = [ 0; 1; 0];
Qp{4} = [ 1; 1; 0];
Qp{5} = [1/2; 1/2; 0];
Qp{6} = [ 0; 0; 0];
Qp{7} = 500;
%Spin wave spetrum
Jfun = @(x)cat(3,diag([-x(4) 0 0]),diag([0 -x(4) 0]),diag([0 0 -x(4)]),eye(3)*x(1),eye(3)*x(2),eye(3)*x(3));
J1 = -1.7; J2 = 0; J3 = 1.8; JK = 0;
cnso.matrix.mat = Jfun([J1 J2 J3 JK]);
cnso.genmagstr('mode','func','func',@gm_planar,'x0',[[3/2 1/2 1/2 3/2 3/2 1/2 1/2 3/2]*pi 0 0 0 0 0]);
specIJ = cnso.spinwave(Qp,'hermit',false);
specIJ = sw_neutron(specIJ);
cnsoPow = cnso.powspec(linspace(0.5,4.2,74),'Evect',linspace(0,31,32),'nRand',1e3,'hermit',false);
%figure
%sw_plotspec(cnsoPow,'axLim',[0 2.2],'colorbar',true)
%colormap(jet)
%title('Calculated Powder INS Spectrum: T = 0 K', 'FontSize', 16, 'FontName', 'Times New Roman')
%xlabel('Q (Å^{-1})', 'FontSize', 14, 'FontName', 'Times New Roman')
%ylabel('E (meV)', 'FontSize', 14, 'FontName', 'Times New Roman')
%set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
%-------------------------------------------------------------------------%
data = load('C:\Users\L26111170\Desktop\Ni2p0.txt');
Q = unique(data(:,1));
E = unique(data(:,2));
I = unique(data(:,3));
e = unique(data(:,4));
qm = unique(data(:,5));
qM = unique(data(:,6));
data = struct('x', E, 'y', I, 'e', e, 'qmin', qm, 'qmax', qM);
fit_func = @(obj, p) matparser(obj, 'param', p, 'mat', {'J1', 'J2', 'J3'}, 'init', true);
J1_guess = -1.0;
J2_guess = 1.0;
J3_guess = 1.0;
fitpow = sw_fitpowder(cnso, data, fit_func, [J1_guess, J2_guess, J3_guess]);
% 設定儀器參數
fitpow.powspec_args = struct('nRand', 1e3, 'hermit', true, 'fastmode', false);
fitpow.sw_instrument_args = struct('dQ', 0.5, 'dE', @(E) 0.1*sqrt(E), 'Ei', 25);
% energy range
fitpow.crop_energy_range(0.5, inf); % fit E > 0.5 meV
% parameter bounds
fitpow.set_model_parameter_bounds(1, -10, 0); % J1 < 0
fitpow.set_model_parameter_bounds(2, 0, 10);
fitpow.set_model_parameter_bounds(3, 0, 10);
[p_fit, chi2, stat] = fitpow.fit(); % fitting
disp(['J1 = ', num2str(p_fit(1)), ' meV']);
disp(['J2 = ', num2str(p_fit(2)), ' meV']);
disp(['J3 = ', num2str(p_fit(3)), ' meV']);
disp(['χ² = ', num2str(chi2)]);
% plot 2D
%fitpow.plot_result(p_fit, 'EdgeAlpha', 0.5);
% plot 1D cut
%q_cuts = [0.5, 1.0, 1.5];
%dq = 0.1;
%fitpow.plot_1d_cuts_of_2d_data(q_cuts-dq, q_cuts+dq, p_fit);
% residual calculation
%residual = (data.y - fitpow.model_intensity) ./ data.e;
%imagesc(Q, E, residual'); colorbar;
%xlabel('Q (Å⁻¹)'); ylabel('E (meV)'); title('殘差圖');
0 Commenti
Risposte (1)
Aryan
il 22 Ago 2025
Hi Wen-Chun,
I understand that you want to know the exact data format that 'sw_fitpowder' in SpinW expects.
As in your case for a 1-D constant-Q cut, the data struct should be formatted like below:
data = struct( ...
'x', [E1 E2 ... EN], % vector of energy values
'y', [I1 I2 ... IN], % vector of intensities (same length as x)
'e', [e1 e2 ... eN], % vector of errors (same length as x)
'qmin', q_lower, % scalar: min Q of the cut
'qmax', q_upper ); % scalar: max Q of the cut
In your current script, the mistake is using unique(...), which destroys the actual mapping between each energy bin and its measured intensity/error. You should instead keep the full vectors and assign qmin/qmax as scalars.
You can also refer to the below link to understand how the formatting should look like:
Hope it helps !
0 Commenti
Vedere anche
Categorie
Scopri di più su Quantum Mechanics in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!