Bug in optmagsteep when assigning moments to sublattices

[henrikjacobsenfys]

If some spins are not couple to any other spins, optmagsteep can give errors. I have traced the error to the call to sw_fSub. In the example below, there are 64 spins in the magnetic cell, but only 32 of them are coupled to other spins. The highest index of these spins happens to be 60, and so the code that assigns spins to sublattices only assigns the first 60 spins:

SSc = param.fSub(SS.all(4:5,:),param.nExt);

This then leads to errors.

Here is a (somewhat) minimal example:

clear
close all

%% Bi2Fe4O9
s = spinw;
s.genlattice('lat_const',[7.99 8.43 6.01], 'sym', 'P b a m');% According to Duc et al.
 
s.addatom('r',[0.3540+0.5 0.33610 0.5000],'S', 5/2,'label','Fe1','color','black');
s.addatom('r',[0.0000+0.5 0.50000 0.2582],'S', 5/2,'label','Fe2','color','white');
 
%% Insert exchange values
 
s.gencoupling('maxDistance', 15)

s.addmatrix('label','J4','value', 1,'color','b');
s.addcoupling('mat', 'J4', 'bond', 5);
 
%% Magntic structure

% Rotation between each Fe1-Fe1 or Fe2-Fe2 spin is 90. But Fe1-Fe2 is 155
% in pentagon
M1=rotz(-90);
M2=rotz(155);
 
V1=[-0.2 1 0]';
V2=M2*V1;

% The white spins rotate in all unit cells (Fe2), the black (Fe1) stays the same 
S0 = [V1 -V1 M1*V1 -M1*V1 -V2 -M1*V2 -V2 -M1*V2];
S1 = [V1 -V1 M1*V1 -M1*V1 V2 M1*V2 V2 M1*V2];
 
s.genmagstr('mode', 'direct', 'nExt', [2 2 2], 'S', [S0 S1 S1 S0 -[S0 S1 S1 S0]]); 

% Optimize the spin orientation to lowest energy
s.optmagsteep('nRun', 50000)