McStasMcXtrace / iFit

a simple library to analyze data (with McCode and Phonons/DFT hooks). :warning: this project has been moved to https://gitlab.com/soleil-data-treatment/soleil-software-projects/remote-desktop
http://ifit.mccode.org
Other
5 stars 5 forks source link

ResLibCal conv with e.g. SpinW #160

Closed farhi closed 6 years ago

farhi commented 6 years ago

Bug found by H. Jacobsen Feb 27th 2018:

Hi Emmanuel

I am trying to convolve a SpinW with the resolution of IN20, and am running into problems. I am using the newest versions of both programs. My script is attached along with the .mat file that is generated. I am also getting problems when I run the sqw_sine3d example in http://ifit.mccode.org/Applications/ResLibCal/doc/ResLibCal.html#mozTocId244259. I get problems both with and without the resolution convolution.

For this I used Matlab R2017a on Windows.

I could also use some help settings the convolution up properly, if you don't mind, as it is a bit confusing to me how the instrument parameters are to be included in the calculation.

Thanks!

Henrik

which produces the output:

Creating the bond list (maxDistance = 8 Å, nCell = 3x3x2)... ...220 bonds are retained out of 2240 generated! Warning: Unregognised input parameter "unitS"!

In sw_readparam (line 143) In spinw/genmagstr (line 221) In SpinW_resolution_convolution (line 16) sqw_spinw: Model Cu1 built using SpinW. Ground state energy: -38.8075 [mev/spin]

  • S. Toth and B. Lake, J. Phys.: Condens. Matter 27, 166002 (2015). 26-Feb-2018 23:04:05: Loaded ResLibCal from C:\Users\Henrik\AppData\Roaming\MathWorks\MATLAB\R2017a\ResLibCal.ini ResLibCal_tas_conv4d: Using reference frame: rlu and method: Popovici (ResCal5) [R] rlu lattice frame, xyz=[a b c*] conv(Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw], ResLibCal('struct('method','Pop...')) [ResLibCal_tas_conv4d] Calculating COMMENSURATE spin wave spectra (nMagExt = 16, nHkl = 200, nTwin = 1)... No magnetic form factor is included in the calculated structure factor. No g-tensor is included in the calculated structure factor.

Error: Could not evaluate Expression in model Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw] iF541785 this = iFunc (methods,doc,plot,code) 4D model: "Sqw_spinw spin-wave dispersion(HKL) [s..." Expression: ax = {x y z t};check_vector = cellfun(@(c)(~isempty(c) && length(c)... Description: A HKL spin-wave dispersion from SpinW package with Gaussian energy width: Cu1 Jperp Jb Jch Jprime D Tag: 'iF541785' Date: '26-Feb-2018 23:04:05' Name: 'Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw]' Parameters: {1×8 cell} Guess: [0.3000 0 1 3.6500 1.1400 143.6000 2.9300 0.9800] Constraint: [1×1 struct] Dimension: 4 ParameterValues: [8×1 double] UserData: [1×1 struct] Duration: 0

Parameters (8): p( 1)= Gamma=0.1 % entered as: energy broadening around spin-wave modes [meV] p( 2)= Temperature=200 % entered as: [K] p( 3)= Amplitude=1 p( 4)= Jperp=-3.65 p( 5)= Jb=-1.14 p( 6)= Jch=143.6 p( 7)= Jprime=-2.93 p( 8)= D=-0.98 '' 'ax = {x y z t};' 'check_vector = cellfun(@(c)(~isempty(c) && length(c) == numel(c)), ax);' 'check_numel = cellfun(@numel, ax);' 'check_orient = zeros(1,4);' 'index = find(check_vector);' 'check_orient(index) = cellfun(@(c)find(size(c)==numel(c)), ax(index));' 'HKL = []; is_event = false; resize_me = [];' 'if all(check_vector) && all(check_numel == check_numel(1)) && all(check_orient == check_orient(1));' 'is_event = true;'

... 'spec = spinwave(this.UserData.spinw, HKL');' 'spec = sw_neutron(spec);' 'this.UserData.maxFreq=max(spec.omega(:));' 'spec = sw_egrid(spec,'component',this.UserData.component,'Evect',t(:)', 'T', p(2));' 'spec = sw_instrument(spec,'dE',p(1),'ki',this.UserData.ki);' 'signal = reshape(spec.swConv,resize_me([4 1:3])); signal=p(3)*permute(signal,[2:4 1]);'

Error using spinw/spinwave (line 819) Hamiltonian matrix is not positive definite, probably the magnetic structure is wrong! For approximate diagonalization try the param.hermit=false option

ans =

'J:\Dropbox\Postdoc\CuO\resolution\iFunc_feval_error'

Warning: Escaped character '\D' is not valid. See 'doc sprintf' for supported special characters.

In iFunc/feval>iFunc_feval_expr (line 429) In iFunc/feval (line 343) In iFunc/feval>iFunc_feval_expr (line 413) In iFunc/feval (line 343) In iData>iData_iFunc2iData (line 247) In iData (line 173) In SpinW_resolution_convolution (line 51)

I attached the config file that I (think that) I used.

Also, it looks like the sign of the coupling parameters J gets positive disregarding how they are input (-J). This could be the issue.

The log indicates that the actual error is:

    Error using spinw/spinwave (line 819)
    Hamiltonian matrix is not positive definite, probably the magnetic structure is wrong! For approximate diagonalization try the param.hermit=false option

It might be that the parameters are read wrong - it seems that the negative signs for some of the interactions disappear. Here's some of the output from iFit: Parameters (8): p( 1)= Gamma=0.3 % entered as: energy broadening around spin-wave modes [meV] p( 2)= Temperature=0 % entered as: [K] p( 3)= Amplitude=1 p( 4)= Jperp=3.65 p( 5)= Jb=1.14 p( 6)= Jch=143.6 p( 7)= Jprime=2.93 p( 8)= D=0.98
And here are the parameters as I set them: Jch=143.6; Jperp=-3.65; Jb=-1.14; Jprime=-2.93; D=-0.98; t.ParameterValues=[0.1 200 1 Jperp Jb Jch Jprime D];

farhi commented 6 years ago

Message from David Voneshen - UKRI STFC david.voneshen@stfc.ac.uk on Aug 24th 2018 2:48PM:

Hi Emmanuel,

I am trying to use ResLibCal (latest) and spinw (R1321) within MATLAB 2015b to model some data we recently collected on IN8. My issues seem to be almost identical to the iFit bug #160 reported by Henrik Jacobsen earlier this year.

Unlike Henrik I am able to use ResLibCal and the simple example for sqw_sine3d with some minor modifications. My sine3d script looks like

s=sqw_sine3d(1.); % create a 4D S(q,w) for a cubic pure material, using default parameters
t=ResLibCal(s);         % convolute it with a TAS resolution, and open ResLibCal.
w=linspace(0.5,2.0,50); qh=1.1*ones(size(w)); qk=0*qh; ql=0*qh; % an E-scan around the [1 0 0] Bragg
signal0=iData(s, [], qh,qk,ql,w);   % the ideal model evaluated along the HKLE line
signal1=iData(t, [], qh,qk,ql,w);   % the convoluted model evaluated along the HKLE line
plot(w,signal0.signal,w,signal1.signal*120)

When trying to use spinW however I seem to get the exact same issue as Henrik. The minimum code required to reproduce my problem is

s=sqw_spinw()
t=ResLibCal(s);
w=linspace(0,10,101);
qh=1.15*ones(size(w)); qk=0*qh; ql=0*qh;
cut=iData(t, t.p, qh,qk,ql,w)
plot(w,cut.signal)

Which produces the following output

1.10 iFit/iData (Aug. 22, 2017) by E.Farhi, P. Willendrup and Y.Debab, (c) ILL/DS/CS <farhi@ill.eu> EUPL. Aug. 22, 2017
Opening iData documentation from
  web C:\Users\tyo59427\Documents\projects\spin-cal\Magnetite\iFit\ifit-1.10\Objects\@iData\..\..\\Docs\Models.html#mozTocId908192
Calculation is finished in 00:00:00 (hh:mm:ss).
sqw_spinw: Model  atom_1 built using SpinW.
Ground state energy: -4 [mev/spin]
* S. Toth and B. Lake, J. Phys.: Condens. Matter 27, 166002 (2015).
s =  iFunc (methods,doc,plot,more...) 4D model:

    [Tag] [Dim]                                [Model] [Parameters 'p']
iF488838   4  Sqw_spinw spin-wav. A HKL spi. ax = ... Gamma Temperature Amplitude J1
ResLibCal_tas_conv4d: Using reference frame: rlu and method: unknown
[R] rlu lattice frame, xyz=[a* b* c*]
conv(Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw], ResLibCal('struct('Title','ResL...')) [ResLibCal_tas_conv4d]
Calculation is finished in 00:00:00 (hh:mm:ss).
Error: Could not evaluate Expression in model Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw] iF488838
this = iFunc (methods,doc,plot,code) 4D model: "Sqw_spinw spin-wave dispersion(HKL) [s..."
         Expression: ax = {x y z t};check_vector = cellfun(@(c)(~isempty(c) && length(c)...
        Description: A HKL spin-wave dispersion from SpinW package with Gaussian energy width:  atom_1 J1
                Tag: 'iF488838'
               Date: '24-Aug-2018 11:38:50'
               Name: 'Sqw_spinw spin-wave dispersion(HKL) [sqw_spinw]'
         Parameters: {'Gamma energy broadening around spin-wave modes [meV]'  'Temperature [K]'  'Amplitude'  'J1'}
              Guess: [0.3000 0 1 2]
         Constraint: [1x1 struct]
          Dimension: 4
    ParameterValues: [4x1 double]
           UserData: [1x1 struct]
           Duration: 0

Parameters (4):
  p(  1)=               Gamma=0.3   % entered as: energy broadening around spin-wave modes [meV]
  p(  2)=         Temperature=0   % entered as: [K]
  p(  3)=           Amplitude=1
  p(  4)=                  J1=2
    ''
    'ax = {x y z t};'
    'check_vector = cellfun(@(c)(~isempty(c) && length(c) == numel(c)), ax);'
    'check_numel  = cellfun(@numel, ax);'
    'check_orient = zeros(1,4);'
    'index        = find(check_vector);'
    'check_orient(index) = cellfun(@(c)find(size(c)==numel(c)), ax(index));'
    'HKL = []; is_event = false; resize_me = [];'
    'if all(check_vector) && all(check_numel == check_numel(1)) && all(check_orient == check_orient(1));'
    'is_event = true;'

...
    'spec = spinwave(this.UserData.spinw, HKL');'
    'spec = sw_neutron(spec);'
    'this.UserData.maxFreq=max(spec.omega(:));'
    'spec = sw_egrid(spec,'component',this.UserData.component,'Evect',t(:)', 'T', p(2));'
    'spec = sw_instrument(spec,'dE',p(1),'ki',this.UserData.ki);'
    'signal = reshape(spec.swConv,resize_me([4 1:3])); signal=p(3)*permute(signal,[2:4 1]);'

Index exceeds matrix dimensions.

Warning: Control Character '\U' is not valid. See 'doc sprintf' for control characters valid in the format string.
> In iFunc/feval>iFunc_feval_expr (line 428)
  In iFunc/feval (line 343)
  In iFunc/feval>iFunc_feval_expr (line 413)
  In iFunc/feval (line 343)
  In iData>iData_iFunc2iData (line 247)
  In iData (line 173)
Warning: iData/setaxis: the Alias Axis_1 used to define axis rank 1 does not exist in object out iD592717 "".
      Defining it.
cut =  iData (methods,doc,plot,plot(log),values,more...) 4D object:

    [Tag] [Dimension]                                     [Title] [Last command]          [Label/DisplayName]
iD592717     [1 101] 'conv(Sqw_spinw spin-wave di... "Data ..."' iD592717=setalias(iD... conv(Sqw_s/conv(Sqw_s...

I can if you wish also send my actual spinw model and code but the error is the same. I’ve also tried using the internal version of spinw but it produces the same error. Is there a fix for this? Or will I need to move to use a different package?

Best wishes, David Voneshen P.S. the requested smile J

farhi commented 6 years ago

Nailed down error:

comes from the last 'resize' call in Script/Models/Factory/sqw_spinw, last in-line:

signal = reshape(spec.swConv,resize_me([4 1:3])); signal=p(3)*permute(signal,[2:4 1]);

as in some cases, the ignal is not 4D. In the examples above, it is 2D. We thus change the line to take into account the actual dimensionality of the signel:

signal = reshape(spec.swConv,resize_me([end 1:(end-1)])); signal=p(3)*permute(signal,[2:ndims(signal) 1]);