SpinW / spinw

SpinW Matlab library for spin wave calculation
http://www.spinw.org
GNU General Public License v3.0
36 stars 15 forks source link

Error when calculating intensities in sw_plotspec #58

Closed henrikjacobsenfys closed 2 years ago

henrikjacobsenfys commented 3 years ago

The normalized intensity changes with the length of Evect, which I don't believe it should. It can be fixed by adding the following to line 636:

fG=fG/(spectra.Evect(2)-spectra.Evect(1));

(Assuming a uniform spacing of the E vector)

An example of how the calculation goes wrong can be seen with this code:

MnF2 = spinw; a = 4.873; c=3.130; MnF2.genlattice('lat_const',[a a c],'angled',[90 90 90],'spgr','P 42/m n m') MnF2.addatom('label','Mn2+','r',[1/2 1/2 1/2],'S',5/2) % MnF2.addatom('label','F-','r',[0 0 0],'S',0)

MnF2.gencoupling;

MnF2.addmatrix('label','J1','value',-0.026812,'color','b') %FM MnF2.addmatrix('label','J2','value',0.155852,'color','r') %AFM

MnF2.addmatrix('value',diag([0 0 -0.091/(5/2)^2]),'label','D')

MnF2.addcoupling('mat','J1','bond',1) MnF2.addcoupling('mat','J2','bond',2) MnF2.addaniso('D') MnF2.genmagstr('mode','direct','S',[0 0 ; 0 0 ; -1 1 ])

Q{1}=[-1 0 -1]; Q{2}=[-1 0 1]; Q{end+1}=81;

MnF2.unit.nformula=uint8(2); %Get result per formula unit

MnF2Spec = MnF2.spinwave(Q,'formfact',1,'gtensor',1); MnF2Spec = sw_neutron(MnF2Spec); MnF2Spec.norm=1;

figure MnF2Spec = sw_egrid(MnF2Spec,'Evect',linspace(0,8,137),'component','Sperp'); sw_plotspec(MnF2Spec,'axLim',[0 0.5],'mode',3,'dE',0.5); caxis([0 30]) ylim([-0.4 8])

figure MnF2Spec = sw_egrid(MnF2Spec,'Evect',linspace(0,8,1137),'component','Sperp'); sw_plotspec(MnF2Spec,'axLim',[0 0.5],'mode',3,'dE',0.5); caxis([0 30]) ylim([-0.4 8])

These figures should look identical, but they don't. In addition, the integrated intensity (along E) does not match Sperp. With the fix it does.

mducle commented 3 years ago

Yeah, I think the intensity is calculated in barns rather than barns/meV which is what it should be (although I think we don't explicitly say anywhere in the docs what the intensity unit is...). I'll have a dig...

henrikjacobsenfys commented 3 years ago

Yes, it's that, and that the convolution doesn't take the step size of Evect into account. Doubling the step size also doubles the integrated intensity with the current method. With my suggested fix the integrated intensity (calculated e.g. using trapz) of the signal matches the value stored in Sperp.

Also, I think that the unit for S should just be meV^-1, although I'm not certain.

Just to clarify, I suggest the following on line 634 in sw_plotspec

if param.norm fG = fG/sum(fG); fG=fG/(spectra.Evect(2)-spectra.Evect(1)); else fG = fG/max(fG(:)); end

mducle commented 3 years ago

So the intensity in SpinW is calculated in natural units (not barns). We should have an option to output it as a cross-section, specifically as mb/sr/meV/f.u., which involves multiplying by the total neutron-magnetic cross-section, dividing by the energy bin size and dividing by the number of formula units per unit cell Z.

The total neutron-magnetic cross-section is (r0/2)^2=291mb where r0=gn*re/2 and gn=-3.826 is the neutron g-factor and re=0.282e-14m is the classical electron radius.

henrikjacobsenfys commented 3 years ago

I see! Well, this is embarrassing, I have been telling people the wrong thing for a while now, I was sure it was in barns. Adding an option to output as a cross-section would be very useful.

Thanks for looking into this!

mducle commented 3 years ago

@henrikjacobsenfys yeah... sorry I only just realised just now whilst chasing up a question from a user.

GoldenwhaleD commented 3 years ago

@mducle Hi, sw_instrument(MnF2Spec ,'dE',0.5,'norm',1) seems to be able to convert the intensity to absolute unit (mb/meV/f.u.) to me. Could you check this? Thanks.

wardsimon commented 2 years ago

It should be:

 sw_instrument(MnF2Spec ,'dE',0.5,'norm',true)

to convert to mbarn. Note that a 'g' factor of 2 is applied which might not be the case for your system.

mducle commented 2 years ago

@GoldenwhaleD ah thanks for pointing this out.