How to fit powder inelastic neutron scattering data

[l26111170-lab]

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).

Ni2p0.txt

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('殘差圖');

[RichardWaiteSTFC]

Hi @l26111170-lab - thanks for getting in touch, sorry for the delay.

So one issue is that the min/max |Q| in the file are the same - were the data not integrated between some |Q| range or was it a constant |Q| scan?

If you want to do the fitting just with the calculation performed at |Q| = 2 Ang^-1 then I can suggest setting nQ = 1 in the sw_fitpowder constructor (it's the last argument so unfortunately you need to explicitly pass the background strategy as well - it should be "independent" as you have a single 1D cut). You might also need to set |Q| min < max - just pick some symmetrical region centred on 2 Ang^-1.

Another consideration is that the fit might not be well constrained with a single |Q| cut - is this a subset of your data?