How to fit powder inelastic neutron scattering data

4 views (last 30 days)
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('殘差圖');

Answers (1)

Aryan
Aryan on 22 Aug 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 !

Categories

Find more on Quantum Mechanics in Help Center and File Exchange

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!