Test genmagstr

[rebeccafair]

There's not so many lines of code in this one, but lots of branching which has made it a bit tricky. I've found a fairly big list of things while doing this that we might want to fix, document or make an issue for:

• k and S are regarded as 'soft' parameters so they don't cause an error if they are the wrong shape (I assume this is to allow flexibility e.g. in helical allows single or multiple spins) . Although this can lead to silently setting k or S to [] if they are the wrong shape, this can happen easily if you get the S/k dimensions the wrong way round e.g. with swobj.genmagstr('mode', 'direct', 'S', ones(1, 3), 'k', ones(3, 1)). Should we warn/error if this happens?
• random mode actually reads k even though it says it doesn't read anything in the docstring. It doesn't actually do anything with it, just sets it in the output struct. Is this intended? Should I update the docstring?
• direct mode doesn't use n even though it says it does in the docstring. Should I update this?
• If tile mode is chosen and no S is provided, it just sets the mode to random. Isn't this just the same as choosing random in the first place? We should either warn that this is happening or require S to be provided with tile.
• I guess the if ~cmplxS && ~strcmpi(param.mode,'fourier') && ~strcmpi(param.mode,'direct') && any(k(:)) line is converting the input spins into the Fourier components? The cmplxS check is a bit strange, is that implicitly checking whether S is already a Fourier component? But a user could just input a complex S?
• tile mode sums over the third dimension of S (which is the multiple-k dimension), but without actually requiring k as input, is this the behaviour we want? (see test_tile_multik)
• rotate mode uses k from existing magnetic structure, and silently ignores the k input, is this desired?
• extend mode is undocumented and sets the mode to tile do we want to keep this or remove it?
• In rotate if the first spin is parallel to n and phi=0 we get nRot=[0 0 0] and S = NaN NaN NaN. We should probably error in this case?
• If nExt is a scalar, nExt is automatically determined from k. But what is the point in this? You still have to provide the correct number of spins so you could figure out nExt yourself? See test test_direct_multik_scalar_nExt.
• In rotate mode, if phi is imaginary it appears to rotate first spin to [1 0 0] (I don't really understand the code here though). This is undocumented. Should we document or remove this? See test_rotate_complex_phi.

Probably a lot of this could be fixed by some better validation, which could be made easier if each mode was in its own function with its own validation and genmagstr was just wrapper, but that's probably for another time.

[codecov-commenter]

Codecov Report

Merging #97 (4b57897) into master (66ac715) will decrease coverage by 0.41%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #97      +/-   ##
==========================================
- Coverage   26.40%   25.99%   -0.42%     
==========================================
  Files         235      235              
  Lines       15690    15738      +48     
==========================================
- Hits         4143     4091      -52     
- Misses      11547    11647     +100     

Impacted Files	Coverage Δ
swfiles/@spinw/genmagstr.m	100.00% <100.00%> (+41.44%)	⬆️
swfiles/gm_planar.m	0.00% <0.00%> (-84.62%)	⬇️
swfiles/@spinw/optmagstr.m	15.23% <0.00%> (-71.43%)	⬇️
external/sw_fminsearchbnd.m	0.00% <0.00%> (-61.71%)	⬇️
swfiles/sw_cartesian.m	48.27% <0.00%> (-17.25%)	⬇️
swfiles/@spinw/optmagsteep.m	69.94% <0.00%> (-13.12%)	⬇️
swfiles/sw_readparam.m	83.05% <0.00%> (+1.69%)	⬆️
swfiles/sw_nvect.m	73.91% <0.00%> (+21.73%)	⬆️
swfiles/gm_spherical3d.m	57.69% <0.00%> (+57.69%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[github-actions]

Unit Test Results

       4 files  ±    0       72 suites  +4   7m 48s ⏱️ + 2m 5s
   311 tests +  66     311 ✔️ +  66  0 💤 ±0  0 ❌ ±0 
1 232 runs  +262  1 232 ✔️ +262  0 💤 ±0  0 ❌ ±0 

Results for commit f3c8aa2. ± Comparison against base commit 66ac715.

♻️ This comment has been updated with latest results.

[mducle]

Thanks for doing this! It mostly looks good, but needs a few changes as you suggested in your main comments.

I didn't actually know that SpinW supported multi-k input in genmagstr and looking at the code here I'm not sure that it's doing things correctly, as noted in the comments... but I think it would be good to merge these sets of tests so that we know what the code is doing now - we can then investigate the multi-k and non-unity k and nExt combinations in another code (and probably then have to change some of these tests...).

---

As for your comments:

• k and S are regarded as 'soft' parameters so they don't cause an error if they are the wrong shape (I assume this is to allow flexibility e.g. in helical allows single or multiple spins) . Although this can lead to silently setting k or S to [] if they are the wrong shape, this can happen easily if you get the S/k dimensions the wrong way round e.g. with swobj.genmagstr('mode', 'direct', 'S', ones(1, 3), 'k', ones(3, 1)). Should we warn/error if this happens?

Yes!

• random mode actually reads k even though it says it doesn't read anything in the docstring. It doesn't actually do anything with it, just sets it in the output struct. Is this intended? Should I update the docstring?

Yes! What this means is that it generates a random structure in the first unit cell and then it will apply the rotation frame with the default n normal vector to this. Spin wave calculations will use the k vector given - so the users should be aware of this. (Note also that the n is also read in and used if k is given - it basically behaves like helical except that S is randomly generated instead being specified).

• direct mode doesn't use n even though it says it does in the docstring. Should I update this?

Yes, please!

• If tile mode is chosen and no S is provided, it just sets the mode to random. Isn't this just the same as choosing random in the first place? We should either warn that this is happening or require S to be provided with tile.

I think we should require that an S be provided with tile (replace line 308 with an error).

• I guess the if ~cmplxS && ~strcmpi(param.mode,'fourier') && ~strcmpi(param.mode,'direct') && any(k(:)) line is converting the input spins into the Fourier components? The cmplxS check is a bit strange, is that implicitly checking whether S is already a Fourier component? But a user could just input a complex S?

Yeah, it seems to be a way to implictly avoid using the n vector since this is only used if the check passes (the n vector is not otherwise used in helical mode). It's not documented or used in any examples though, so I'm minded to introduce a check that helical mode requires real S (or a warning that a complex S implies ignoring n if n is specified in helical mode).

• tile mode sums over the third dimension of S (which is the multiple-k dimension), but without actually requiring k as input, is this the behaviour we want? (see test_tile_multik)

I'm not sure... I think we keep the current test and behaviour but revisit it when we tackle multiple-k values (in the check of magstr) and basically in addressing this user request issue

• rotate mode uses k from existing magnetic structure, and silently ignores the k input, is this desired?

I think that rotate mode is meant to be used only after a magnetic structure has already been defined by another mode (certainly this is the case in all the tutorials it is used - 15, 21 (optmagsteep automatically generates a random structure if nothing is defined) and 30.

So I think we should raise an error if no structure has been defined.

Changing the k vector changes the magnetic structure drastically - a uniform rotation of the spin should mean that the spin wave dispersion is still the same (but the intensity at a given q-point might change). Changing k would change the spin wave dispersion. I think that this is not the purpose of the rotate mode (which is mostly used to make the magnetic structure figures look prettier by aligning the moments to an crystallographic axis or other such cosmetic changes). So, I think that if the use specifies a k with the rotate mode it should be an error (or at least a warning that the k vector will be ignored).

• extend mode is undocumented and sets the mode to tile do we want to keep this or remove it?

As discussed I think this should be maintained for now with a deprecation warning on the use of extend which is not used anywhere in the examples / tutorials. I think this was introduced later to be more inline with the syntax of other functions such as intmatrix, but was never propagated. Since tile is more often used I think we should deprecate extend instead of tile in this context.

• In rotate if the first spin is parallel to n and phi=0 we get nRot=[0 0 0] and S = NaN NaN NaN. We should probably error in this case?

Yes, I think this should raise an error.

• If nExt is a scalar, nExt is automatically determined from k. But what is the point in this? You still have to provide the correct number of spins so you could figure out nExt yourself? See test test_direct_multik_scalar_nExt.

This is used to convert a structure from a rotating frame to a supercell and is used in tutorial 18 (there is a formatting error in this tutorial we should fix). So the idea is that you first create a structure with helical or fourier and do an incommensurate spin wave calculation then want to check if it is correct/equivalent by converting to a super cell with a scalar nExt and doing a commensurate calculation. In some cases of complex structures, the incommensurate calculation might not be correct...

• In rotate mode, if phi is imaginary it appears to rotate first spin to [1 0 0] (I don't really understand the code here though). This is undocumented. Should we document or remove this? See test_rotate_complex_phi.

I don't understand this behaviour either... it seems to be an undocumented short hand to align the structure so that components the first spin perpendicular to n is parallel to [1 0 0] but it doesn't work if the first spin is parallel to n.

I think this might just be a weird short hand which Sandor introduced for his own use which we can just remove (it's not documented anywhere or used in the code anywhere else...). So maybe we can just remove this code and raise an error if an imaginary / complex input is given for phi or phid - and see if any users complain ;)

[RichardWaiteSTFC]

Apologies if this has already been covered, but in helical mode does it make sense for n to be parallel to S for example

testCase.swobj = spinw;
testCase.swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]);
testCase.swobj.addatom('r', [0 0 0],'S', 1, 'label', 'MNi2');
testCase.swobj.gencoupling('maxDistance', 7);
testCase.swobj.addmatrix('value', -eye(3), 'label', 'Ja');
testCase.swobj.addcoupling('mat', 'Ja', 'bond', 1);
testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', [1,0,0], 'k', [1/3 0 0])
testCase.swobj.plot('range', [3, 1, 1])

produces

I didn't really know what to expect, but it looks like maybe it has to zero the moment in cells that aren't an integer number of wavelengths along the modulation direction to stay consistent with k? If that makes sense...

[mducle]

@RichardWaiteSTFC

So, for helical structures, before it gets to the helical block in line 376, it adds a imaginary component perpendicular to the real component in line 342 in the basisvector F, which as explained here results in a helical structure.

Now, if n is parallel to S for some S then that perpendicular imaginary component is zero. In your example there is only one spin which is parallel to n so the basisvector is entirely real, giving an amplitude modulated spin structure. This does not necessarily result in zero moments on sites which are not integral number of unit cells of wavelengths of k. E.g. if in your example you then added:

mgstr = testCase.swobj.magstr('nExt', [3 1 1])
mgstr.S

You'd see that the intermediate spins are half of the main length rather than zero (they should in principle follow a sinusoidal dependence, S*cos(i/n*pi) where i is the index and n is the wavelength 1/k).

Edit: I'm not sure if we should warn users about this or not...

[RichardWaiteSTFC]

@RichardWaiteSTFC

So, for helical structures, before it gets to the helical block in line 376, it adds a imaginary component perpendicular to the real component in line 342 in the basisvector F, which as explained here results in a helical structure.

Now, if n is parallel to S for some S then that perpendicular imaginary component is zero. In your example there is only one spin which is parallel to n so the basisvector is entirely real, giving an amplitude modulated spin structure. This does not necessarily result in zero moments on sites which are not integral number of unit cells of wavelengths of k. E.g. if in your example you then added:

mgstr = testCase.swobj.magstr('nExt', [3 1 1])
mgstr.S

You'd see that the intermediate spins are half of the main length rather than zero (they should in principle follow a sinusoidal dependence, S*cos(i/n*pi) where i is the index and n is the wavelength 1/k).

Edit: I'm not sure if we should warn users about this or not...

Thanks @mducle I think I get it... but it can lead to some unexpected (at least to me) results! So if I understand correctly, if n has even a component parallel to the direction of the spin S (in this case there is only one per unit cell) then it results in a magnitude modulation. I think reading the docs one might expect
testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', [1,0,1], 'k', [1/3 0 0])
to result in a structure where the spins rotate around the 101 direction but in fact the spin rotates around the 001 direction (and there is a modulation in the moment size)

[mducle]

@RichardWaiteSTFC I think that by definition, in a helical magnet the moments must lie on a plane. In SpinW, n is taken as normal to this plane so if you specify an n which has a component parallel to any of your spins then there is an inconsistency which can either be resolved by ignoring the component parallel to S of n or the component of S parallel to n. SpinW choses the former, e.g. in your example, it silently makes n=[0,0,1] instead of the [1,0,1] you specified. The alternative is to keep n=[1,0,1] but to move the spin so that it is S=[0,1,0] (or similar) instead, which I think is a worst solution that changing n silently...

Maybe we should emit a warning if there is a component of n parallel to S to say that it would be ignored for that spin?

Now, if what you really want with S=[1,0,0] and n=[1,0,1] is a conical structure with a constant component parallel to [1,0,0] (e.g. a net magnetic moment) but with a fluctuating component along [0,1,0], then unfortunately you have to specify this by hand in SpinW. Also, this structure cannot be described in the rotating coordinate frame (I think!) so we cannot use the fast incommensurate algorithm to compute its spin wave spectrum.

[RichardWaiteSTFC]

@RichardWaiteSTFC I think that by definition, in a helical magnet the moments must lie on a plane.

Good point - the structure I was talking about wasn't helical in the sense that it is normally understood to be rotating in a plane orthogonal to S!

In SpinW, n is taken as normal to this plane so if you specify an n which has a component parallel to any of your spins then there is an inconsistency which can either be resolved by ignoring the component parallel to S of n or the component of S parallel to n. SpinW choses the former, e.g. in your example, it silently makes n=[0,0,1] instead of the [1,0,1] you specified. The alternative is to keep n=[1,0,1] but to move the spin so that it is S=[0,1,0] (or similar) instead, which I think is a worst solution that changing n silently...

Agreed changing the S would be weird - although it isn't simply ignoring the component of n // S (which would be nice and easy to explain) - as you point out above it modulates the moment.

Maybe we should emit a warning if there is a component of n parallel to S to say that it would be ignored for that spin?

I think a warning would be helpful, just to say the plane of the moments will be normal to the vector n minus the component parallel to S and the moment size will be modulated. I don't propose stopping support for defining a helical structure with a modulated moment size (although I'm not sure how intuitive the helical interface is for that, you'd probably use one of the other modes of genmagstr...).

[RichardWaiteSTFC]

Really minor but would it be possible to clarify the docstring the say n will be normalised to a unit vector (if that is what is meant here)

spinw/swfiles/@spinw/genmagstr.m

Lines 139 to 142 in 8c2a818

	% `'n'`
	% : Normal vector to the spin rotation plane for single-k magnetic
	% structures, stored in a 3-element row vector. Default value `[0 0 1]`. The
	% coordinate system of the vector is determined by `unit`.

[RichardWaiteSTFC]

The code coverage indicates this branch hasn;t been hit

spinw/swfiles/@spinw/genmagstr.m

Lines 329 to 332 in 8c2a818

	if size(param.n,1) ~= nK
	error('spinw:genmagstr:WrongInput',['The number of normal vectors has'...
	' to be equal to the number of k-vectors!'])
	end

But it looks like it is impossible to throw this error due to the size validation in the input to sw_readparam - perhaps it's worth deleting?

[RichardWaiteSTFC]

There seems to be validation around ensuring the supercell (if using) has to match the k (e.g. k=[1/3, 0, 0] is not compatible with nExt=[2,0,0]) - but it is possible to chose an artificially large supercell e.g.
testCase.swobj.genmagstr('mode', 'fourier','S', [1; 1i; 0], 'k', [1/3 0 0], 'nExt', [6, 1, 1])
testCase.swobj.genmagstr('mode', 'fourier','S', [1; 1i; 0], 'k', [1/3 0 0], 'nExt', [3, 1, 3])
Is it possible to throw a warning - at least for the latter case when the unit cell is enlarged along a direction for which the component of the modulation is 0?

[RichardWaiteSTFC]

Should probably throw a warning if k is provided with tile mode to say it is ignored - this sort of thing might apply to other modes with ignored parameters

[RichardWaiteSTFC]

There seems to be a problem if you specify a complex n for a helical structure - obviously this is silly...but perhaps we need to throw an error? If you specify an imag. vector it is OK (because post normalisation to unit vector it becomes real) - but a complex one causes issues

testCase.swobj = spinw;
testCase.swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]);
testCase.swobj.addatom('r', [0 0 0],'S', 1, 'label', 'MNi2');
testCase.swobj.gencoupling('maxDistance', 7);
testCase.swobj.addmatrix('value', -eye(3), 'label', 'Ja');
testCase.swobj.addcoupling('mat', 'Ja', 'bond', 1);
testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', [1,1i,0], 'k', [1/3 0 0])
mgstr = testCase.swobj.magstr('nExt', [3 1 1]);

The call to magstr just hangs and when I ctrl-c I get this

In spinw/magstr (line 130)
nUn =
sw_uniquetol(reshape(cross(real(obj.mag_str.F(:,:,kInc)),imag(obj.mag_str.F(:,:,kInc))),3,[])); 

Note that the n unit vector in genmagstr becomes non-finite
n = bsxfunsym(@rdivide,param.n,sqrt(sum(param.n.^2,2))) % [Inf + 0.0000i, 0.0000 + Infi, NaN + 0.0000i]

[RichardWaiteSTFC]

Great job, I think it's probably worth adding a few things:

1. A test with epsilon explicitly provided (see suggestion in attached comments)
2. A test with norm=false (i.e. taking the magnitude of the spin from S input) e.g.

testCase.swobj = spinw;
testCase.swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]);
testCase.swobj.addatom('r', [0 0 0],'S', 2, 'label', 'MNi2');
testCase.swobj.gencoupling('maxDistance', 7);
testCase.swobj.addmatrix('value', -eye(3), 'label', 'Ja');
testCase.swobj.addcoupling('mat', 'Ja', 'bond', 1);
testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', [0,1,0], 'k', [1/3 0 0], 'norm', false)
mgstr = testCase.swobj.magstr();
S = mgstr.S  % [1,0,0]
testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', [0,1,0], 'k', [1/3 0 0], 'norm', true)
mgstr = testCase.swobj.magstr();
S = mgstr.S  % [2,0,0]

3. Add a warning (and test) for when n has a component parallel to S in helical mode
4. A warning if n is complex (I don't think it's worth doing such a check for all other parameters - it's just this does cause a crash when magstr is called)

[rebeccafair]

Should probably throw a warning if k is provided with tile mode to say it is ignored - this sort of thing might apply to other modes with ignored parameters

Done in d8d54f3, have extended it to other parameters, probably overkill but wasn't too difficult to implement.

There seems to be a problem if you specify a complex n for a helical structure - obviously this is silly...but perhaps we need to throw an error? If you specify an imag. vector it is OK (because post normalisation to unit vector it becomes real) - but a complex one causes issues

Done in ca2c7af, again probably overkill to check other parameters too but it seemed to not be much more effort...

1. A test with epsilon explicitly provided (see suggestion in attached comments)

Ah good spot I missed that! Done.

2. A test with norm=false (i.e. taking the magnitude of the spin from S input) e.g.

Oops, completely forgot about that parameter, done.

3. Add a warning (and test) for when n has a component parallel to S in helical mode

Done.

4. A warning if n is complex (I don't think it's worth doing such a check for all other parameters - it's just this does cause a crash when magstr is called)

It errors rather than warning to avoid the crash, but also done.

[rebeccafair]

I didn't realise varargin can also be a struct (e.g. as called from optmagstr), which caused the tests to fail. I've fixed it to allow a struct, and for simplicity slightly reverted ca2c7af to only check for complex n specifically, doing any more was overkill anyway.

In addition to varargin potentially being a struct, if it is called from an optimisation function like optmagstr it also contains a lot of extra parameters to do with the optimisation. The optmagstr etc. functions get around this by calling warning('off','sw_readparam:UnreadInput') to avoid sw_readparam causing lots of ugly errors. It also has the full list of parameters S, k, nExt etc. which, depending on the mode, can also cause lots of ugly warnings to be triggered. To get around this, I have renamed the new warning spinw:genmagstr:UnusedInput to spinw:genmagstr:UnreadInput, which seems to mysteriously be covered by warning('off','sw_readparam:UnreadInput') and stops the warnings from being triggered.

[RichardWaiteSTFC]

Looks great, thanks for taking on this tough one!