From 4d02dbdbc44e380f13e892a80ed0a995a6c349fd Mon Sep 17 00:00:00 2001 From: Rebecca Fair Date: Wed, 17 Aug 2022 14:50:33 +0100 Subject: [PATCH] Test genmagstr (#97) * Add initial genmagstr test file * Add initial simple comm/incomm tests * Add zero nExt test * Add afm chain test * Add random structure tests * Check mag_str not magstr * Add random test with k * Add test with existing k * Add multiatom direct test * Add helical tests, fix multi-k bug In the if ~cmplxS... conditional any(k) returns a vector for multi-k structures, causing a MATLAB:nonLogicalConditional error, fix this by flattening k first * Add to helical tests - Add test with multiple k and n - Add test with nSpin = nMagAtom - Add test with nSpin = nMagExt * Add fourier tests * Add tile tests * Add rotate tests * Add custom func test * Add default func test * Add extend and unit tests * Remove unreachable code Test for number of magnetic atoms is already performed at beginning of function * Add scalar nExt test * Add more error tests * Add complex rotate and symbolic tests * Raise error if invalid S/k are provided, improve random tests * Raise error with tile and wrong number of input spins * Improve input documentation * Error if complex S is provided in helical mode * Raise error if rotate mode is used without initialisation * Emit deprecation warning if extend mode is used * Error if n and S parallel in rotate mode with no phi * Raise error if imaginary phi is used in rotate mode * Refactor tests - use TestParameter * Clarify random mode docstring * Warn if n || S for helical, and add tests * Small upates from review comments - Specify n is converted to a unit vector in docstring - Remove unreachable n/k size error check * Warn if nExt is too large in helical, fourier * Warn on unused inputs - Also move 'extend' -> 'tile' to the top so tile can be used in unused input checks - Also move valid mode check to in unused warning checks so an invalid mode isn't searched for in the struct * Raise error if nonsensical complex values are provided - Also move errors before warnings * Add test with epsilon * Put phi and phid into parametrised test * Add test for norm explicitly true/false * Only test for complex n Inspecting varargin is more complex than initially thought, sometimes it can be a struct (e.g. when called from optmagk) * Account for varargin being a struct When called from places like optmagstr, varargin is actually a struct. Also, the struct contains many other additional parameters to do with fitting, so before calling genmagstr in optmagstr, warning('off','sw_readparam:UnreadInput') is called. This stops sw_readparam from issuing warnings. This also stops the `spinw:genmagstr:UnreadInput` warning from being triggered, which is why its name has been changed to UnreadInput, to avoid ugly unecessary warnings. --- .../+unit_tests/unittest_spinw_genmagstr.m | 655 ++++++++++++++++++ CHANGELOG.rst | 29 +- swfiles/@spinw/genmagstr.m | 148 +++- 3 files changed, 793 insertions(+), 39 deletions(-) create mode 100644 +sw_tests/+unit_tests/unittest_spinw_genmagstr.m diff --git a/+sw_tests/+unit_tests/unittest_spinw_genmagstr.m b/+sw_tests/+unit_tests/unittest_spinw_genmagstr.m new file mode 100644 index 00000000..7ee83075 --- /dev/null +++ b/+sw_tests/+unit_tests/unittest_spinw_genmagstr.m @@ -0,0 +1,655 @@ +classdef unittest_spinw_genmagstr < sw_tests.unit_tests.unittest_super + + properties + swobj = []; + swobj_tri = []; + default_mag_str = struct('nExt', int32([1 1 1]), ... + 'k', [0; 0; 0], ... + 'F', [0; 1; 0]); + end + properties (TestParameter) + fm_chain_input_errors = { ... + % varargin, identifier + {{'mode', 'something'}, 'spinw:genmagstr:WrongMode'}; ... + {{'unit', 'something'}, 'spinw:genmagstr:WrongInput'}; ... + % nExt can't be zero + {{'nExt', [0 1 1]}, 'spinw:genmagstr:WrongInput'}; ... + % n must have dimensions nK, 3 + {{'n', ones(2, 3), 'k', [0 0 1/2]}, 'sw_readparam:ParameterSizeMismatch'}; ... + % S with direct must have same number of spins as atoms in nExt supercell + {{'mode', 'direct', 'S', [0; 1; 0], 'nExt', [2 1 1]}, 'spinw:genmagstr:WrongSpinSize'}; ... + % S with helical must have 1 spin, or same number of spins as + % atoms in unit or supercell + {{'mode', 'helical', 'S', [0 1; 1 0; 0 0]}, 'spinw:genmagstr:WrongNumberSpin'}; ... + {{'mode', 'direct', 'S', [1 0 0]}, 'spinw:genmagstr:WrongInput'}; ... + {{'mode', 'direct', 'k', [1/2; 0; 0]}, 'spinw:genmagstr:WrongInput'}; ... + {{'mode', 'direct', 'S', [1 0 0], 'k', [1/2; 0; 0]}, 'spinw:genmagstr:WrongInput'}; ... + {{'mode', 'tile', 'S', [0 0; 1 1; 1 1]}, 'spinw:genmagstr:WrongInput'}; ... + % S with helical must be real + {{'mode', 'helical', 'S', [1.5; 0; 1.5i]}, 'spinw:genmagstr:WrongInput'}; ... + % Rotate mode must first initialise a magnetic structure + {{'mode', 'rotate'}, 'spinw:genmagstr:WrongInput'} + }; + rotate_input_errors = { ... + % varargin, identifier + % rotation angle must be real (previously phi=i had special + % behaviour) + {{'mode', 'rotate', 'phi', i}, 'spinw:genmagstr:ComplexPhi'}; ... + % If no angle is supplied to rotate, the rotation axis is set + % orthogonal to both the first spin and n. If the first spin + % and n are parallel, this should therefore cause an error + {{'mode', 'rotate'}, 'spinw:genmagstr:InvalidRotation'}; ... + }; + ignored_inputs = { ... + % arguments + {'mode', 'random', 'S', [1; 0; 0], 'epsilon', 1}, ... + {'n', [0 1 1], 'mode', 'direct', 'S', [1; 0; 0]}, ... + {'mode', 'tile', 'k', [0 0 1/2], 'n', [0 0 1], 'S', [1; 0; 0]}, ... + {'mode', 'helical', 'k', [0 0 1/3], 'S', [1; 0; 0;], 'x0', []}, ... + {'mode', 'func', 'x0', [pi/2 -pi/4 0 0 [0 0 1/3] pi pi/2], 'unit', 'lu', 'next', [1 2 1]}, ... + {'mode', 'fourier', 'S', [1; i; 0], 'n', [1 0 0]} + }; + complex_n_input = {[1 1+i 0], [i 0 0]}; + rotate_phi_inputs = {{'phi', pi/4}, {'phid', 45}, {'phi', pi/4, 'phid', 90}}; + input_norm_output_F = { + % norm_bool, mag_str.F + {true, [2; 0; -2i]}, ... + {false, [1; 0; -1i]} + }; + end + methods (TestMethodSetup) + function setup_chain_model(testCase) + % Create a simple FM 1D chain model + testCase.swobj = spinw(); + testCase.swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]); + testCase.swobj.addatom('r', [0 0 0],'S', 1); + end + function setup_tri_model(testCase) + % Create a simple triangular lattice model + testCase.swobj_tri = spinw; + testCase.swobj_tri.genlattice('lat_const', [4 4 6], 'angled', [90 90 120]); + testCase.swobj_tri.addatom('r', [0 0 0], 'S', 3/2); + end + end + + methods (Test) + function test_invalid_fm_chain_input_raises_error(testCase, fm_chain_input_errors) + varargin = fm_chain_input_errors{1}; + identifier = fm_chain_input_errors{2}; + testCase.verifyError(... + @() testCase.swobj.genmagstr(varargin{:}), ... + identifier) + end + function test_invalid_rotate_input_raises_error(testCase, rotate_input_errors) + varargin = rotate_input_errors{1}; + identifier = rotate_input_errors{2}; + swobj = copy(testCase.swobj); + swobj.genmagstr('mode', 'direct', 'S', [0; 0; 1]); + testCase.verifyError(... + @() swobj.genmagstr(varargin{:}), ... + identifier) + end + function test_no_magnetic_atoms_raises_error(testCase) + swobj = spinw; + testCase.verifyError(... + @() swobj.genmagstr(), ... + 'spinw:genmagstr:NoMagAtom') + end + function test_complex_n_input_raises_error(testCase, complex_n_input) + testCase.verifyError(... + @() testCase.swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'n', complex_n_input), ... + 'spinw:genmagstr:WrongInput') + end + function test_tile_too_few_S_raises_error(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + testCase.verifyError(... + @() swobj.genmagstr('mode', 'tile', 'S', [0; 1; 1]), ... + 'spinw:genmagstr:WrongInput') + end + function test_ignored_input_warns(testCase, ignored_inputs) + testCase.verifyWarning(... + @() testCase.swobj.genmagstr(ignored_inputs{:}), ... + 'spinw:genmagstr:UnreadInput') + end + function test_rotate_ignored_input_warns(testCase) + swobj = copy(testCase.swobj); + % Need to initialise structure before rotating it + swobj.genmagstr('mode', 'direct', 'S', [1; 0; 0]); + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'rotate', 'k', [0 0 1/3]), ... + 'spinw:genmagstr:UnreadInput') + end + function test_helical_spin_size_incomm_with_nExt_warns(testCase) + swobj = copy(testCase.swobj); + nExt = [2 1 1]; + k = [1/3 0 0]; + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'helical', ... + 'S', [1 0; 0 1; 0 0], ... + 'k', k, ... + 'nExt', nExt), ... + 'spinw:genmagstr:UCExtNonSuff') + expected_mag_str = struct('nExt', int32(nExt), ... + 'k', k', ... + 'F', [1 -1i; 1i 1; 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_spin_size_incomm_with_epsilon_warns(testCase) + swobj = copy(testCase.swobj); + delta = 1e-6; + k = [(delta+1)/3 0 0]; + nExt = [3 1 1]; + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'helical', ... + 'S', [1; 0; 0], ... + 'k', k, ... + 'nExt', nExt, ... + 'epsilon', 0.99*delta), ... + 'spinw:genmagstr:UCExtNonSuff') + expected_mag_str = struct('nExt', int32(nExt), ... + 'k', k', ... + 'F', [1 -0.5-0.866i -0.5+0.866i; ... + 1i 0.866-0.5i -0.866-0.5i; ... + 0 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str, 'rel_tol', 1e-4); + end + function test_fourier_too_large_nExt_warns(testCase) + swobj = copy(testCase.swobj); + nExt = [6 1 1]; + k = [1/3 0 0]; + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'fourier', ... + 'S', [1; 1i; 0], ... + 'k', k, ... + 'nExt', nExt), ... + 'spinw:genmagstr:UCExtOver') + F_rep = [ 1 -0.5-1i*sqrt(3)/2 -0.5+1i*sqrt(3)/2; ... + 1i sqrt(3)/2-0.5i -sqrt(3)/2-0.5i; ... + 0 0 0]; + expected_mag_str = struct('nExt', int32(nExt), ... + 'k', k', ... + 'F', cat(2, F_rep, F_rep)); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_fourier_nExt_wrong_direction_warns(testCase) + swobj = copy(testCase.swobj); + nExt = [2 1 2]; + k = [1/2 0 0]; + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'fourier', ... + 'S', [1; 1i; 0], ... + 'k', k, ... + 'nExt', nExt), ... + 'spinw:genmagstr:UCExtOver') + F_rep = [1 -1; 1i -1i; 0 0]; + expected_mag_str = struct('nExt', int32(nExt), ... + 'k', k', ... + 'F', cat(2, F_rep, F_rep)); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_any_S_parallel_to_n_warns(testCase) + swobj_tri = copy(testCase.swobj_tri); + k = [1/3 0 0]; + testCase.verifyWarning(... + @() swobj_tri.genmagstr('mode', 'helical', ... + 'S', [0; 1; 1], ... + 'k', k), ... + 'spinw:genmagstr:SnParallel') + expected_mag_str = struct( ... + 'nExt', int32([1 1 1]), ... + 'k', k', ... + 'F', [-sqrt(9/8)*1i; sqrt(9/8); sqrt(9/8)]); + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_helical_S2_norm(testCase, input_norm_output_F) + swobj = spinw(); + swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]); + swobj.addatom('r', [0 0 0], 'S', 2); + k = [1/3 0 0]; + swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], 'k', k, ... + 'n', [0 1 0], 'norm', input_norm_output_F{1}); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k'; + expected_mag_str.F = input_norm_output_F{2}; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_direct_fm_chain(testCase) + swobj = copy(testCase.swobj); + swobj.genmagstr('mode', 'direct', 'k', [0 0 0], ... + 'S', [0; 1; 0]); + testCase.verify_obj(testCase.default_mag_str, swobj.mag_str); + + end + function test_direct_fm_chain_nok(testCase) + swobj = copy(testCase.swobj); + swobj.genmagstr('mode', 'direct', 'S', [0; 1; 0]); + testCase.verify_obj(testCase.default_mag_str, swobj.mag_str); + end + function test_direct_multiatom_nExt(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.5], 'S', 2); + S = [1 0 1 -1; ... + 1 1 0 0; ... + 0 -1 0 0]; + nExt = [2 1 1]; + k = [0 1/3 0]; + swobj.genmagstr('mode', 'direct', 'S', S, 'nExt', nExt, 'k', k); + expected_mag_str = struct('nExt', int32(nExt), ... + 'k', k', ... + 'F', [sqrt(2)/2 0 1 -2; ... + sqrt(2)/2 sqrt(2) 0 0; ... + 0 -sqrt(2) 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_direct_multiatom_multik(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.5], 'S', 2); + S_k = [1 0; 1 1; 0 -1]; + S = cat(3, S_k, S_k); + k = [0 1/3 0; 1/2 0 0]; + swobj.genmagstr('mode', 'direct', 'S', S, 'k', k); + F_k = [sqrt(2)/2 0; ... + sqrt(2)/2 sqrt(2); ... + 0 -sqrt(2)]; + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k'; + expected_mag_str.F = cat(3, F_k, F_k); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_direct_multik_scalar_nExt(testCase) + % Test if a scalar is used for nExt it is treated as a + % tolerance to automatically determine nExt + swobj = copy(testCase.swobj); + S_k = [1 0 1 0 1 0; 0 0 1 1 0 0; 0 0 0 1 1 1]; + S = cat(3, S_k, S_k); + nExt = 0.01; + k = [0 1/3+0.5*nExt 0; 1/2+0.5*nExt 0 0]; + swobj.genmagstr('mode', 'direct', 'S', S, 'k', k, 'nExt', nExt); + F_k = [1 0 sqrt(2)/2 0 sqrt(2)/2 0; ... + 0 0 sqrt(2)/2 sqrt(2)/2 0 0; ... + 0 0 0 sqrt(2)/2 sqrt(2)/2 1]; + expected_mag_str = struct('nExt', int32([2 3 1]), ... + 'k', k', ... + 'F', cat(3, F_k, F_k)); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_tri(testCase) + swobj_tri = copy(testCase.swobj_tri); + k = [1/3 1/3 0]; + swobj_tri.genmagstr('mode', 'helical', 'S', [1; 0; 0], ... + 'k', k); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k'; + expected_mag_str.F = [1.5; 1.5i; 0]; + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_helical_tri_n(testCase) + swobj_tri = copy(testCase.swobj_tri); + k = [1/3 1/3 0]; + swobj_tri.genmagstr('mode', 'helical', 'S', [1; 0; 0], ... + 'k', k, 'n', [0 1 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k'; + expected_mag_str.F = [1.5; 0; -1.5i]; + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_helical_tri_lu_unit(testCase) + swobj_tri = copy(testCase.swobj_tri); + swobj_tri.genmagstr('mode', 'helical', 'S', [0; 1; 0], ... + 'k', [1/3 1/3 0], 'unit', 'lu'); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = [1/3; 1/3; 0]; + expected_mag_str.F = [-0.75-1.299038105676658i; ... + 1.299038105676658-0.75i; ... + 0]; + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_fourier_tri(testCase) + swobj_tri = copy(testCase.swobj_tri); + swobj_tri.genmagstr('mode', 'fourier', 'S', [1; 0; 0], ... + 'k', [1/3 1/3 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = [1/3; 1/3; 0]; + expected_mag_str.F = [1.5; 0; 0]; + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_helical_multiatom_nExt_1spin(testCase) + % Test where only 1 spin is provided + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.0], 'S', 2); + k = [0 0 1/2]; + nExt = int32([1 1 2]); + swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], ... + 'nExt', nExt, 'k', k); + expected_mag_str = struct(... + 'k', k', ... + 'nExt', nExt, ... + 'F', [1 2 -1 -2; 1i 2i -1i -2i; 0 0 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_multiatom_nExt_1spin_r0(testCase) + swobj = spinw(); + swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]); + % Need to have nonzero r for first atom for r0 to have an effect + swobj.addatom('r', [0.5 0.5 0.5],'S', 1); + swobj.addatom('r', [0 0 0], 'S', 2); + k = [0 0 1/2]; + nExt = int32([1 1 2]); + swobj.genmagstr('mode', 'helical', 'S', [1; 0; 0], ... + 'nExt', nExt, 'k', k, 'r0', false); + expected_mag_str = struct(... + 'k', k', ... + 'nExt', nExt, ... + 'F', [1 2i -1 -2i; 1i -2 -1i 2; 0 0 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_fourier_multiatom_nExt_nMagAtom_spins(testCase) + % Test where there are the same number of spins provided as in + % the unit cell. Note result is the same as helical with + % complex spins or S=[0 1; 1 0; 0 0] + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.0], 'S', 2); + k = [0 1/2 0]; + nExt = int32([1 2 1]); + swobj.genmagstr('mode', 'fourier', 'S', [-1i 1; 1 1i; 0 0], ... + 'nExt', nExt, 'k', k); + expected_mag_str = struct( ... + 'k', k', ... + 'nExt', nExt, ... + 'F', [-1i 2 1i -2; 1 2i -1 -2i; 0 0 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_multiatom_nExt_nMagAtom_spins(testCase) + % Test where there are the same number of spins provided as in + % the unit cell + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.0], 'S', 2); + k = [0 1/2 0]; + nExt = int32([1 2 1]); + swobj.genmagstr('mode', 'helical', 'S', [0 1; 1 0; 0 0], ... + 'nExt', nExt, 'k', k); + expected_mag_str = struct( ... + 'k', k', ... + 'nExt', nExt, ... + 'F', [-1i 2 1i -2; 1 2i -1 -2i; 0 0 0 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_multiatom_nExt_nMagExt_spins(testCase) + % Test where there are the same number of spins provided as in + % the supercell + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.0], 'S', 2); + k = [0 1/2 0]; + nExt = int32([1 2 1]); + testCase.verifyWarning(... + @() swobj.genmagstr( ... + 'mode', 'helical', ... + 'S', [0 1 0 -1; 1 0 0 0; 0 0 1 0], ... + 'nExt', nExt, 'k', k), ... + 'spinw:genmagstr:SnParallel'); + expected_mag_str = struct(... + 'k', k', ... + 'nExt', nExt, ... + 'F', [-1i 2 0 -2; 1 2i 0 -2i; 0 0 1 0]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_helical_multiatom_multik_multin(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0.5], 'S', 2); + S = cat(3, [1; 0; 0], [0; 0; 1]); + k = [0 1/3 0; 1/2 0 0]; + n = [0 0 1; 0 1 0]; + % Ensure warning is not emitted as there are no S parallel to + % n within a single k + testCase.verifyWarningFree(... + @() swobj.genmagstr('mode', 'helical', ... + 'S', S, ... + 'k', k, ... + 'n', n), ... + 'spinw:genmagstr:SnParallel') + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k'; + expected_mag_str.F = cat(3, ... + [1 1-sqrt(3)*1i; 1i sqrt(3)+1i; 0 0], ... + [1i 2; 0 0; 1 -2i]); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_random_structure(testCase) + swobj = copy(testCase.swobj); + swobj.genmagstr('mode','random'); + mag_str1 = swobj.mag_str; + swobj.genmagstr('mode','random'); + mag_str2 = swobj.mag_str; + % Check structure is random each time - F is different + testCase.verifyNotEqual(mag_str1.F, mag_str2.F); + testCase.verifyEqual(size(mag_str1.F), size(mag_str2.F)); + % Check F size and magnitude + testCase.verifySize(mag_str1.F, [3 1]); + testCase.verify_val(vecnorm(real(swobj.mag_str.F), 2), 1); + % Check imaginary component of F is perpendicular to default n + testCase.verify_val(dot(imag(swobj.mag_str.F), [0 0 1]), 0); + % Check other fields + expected_mag_str = testCase.default_mag_str; + testCase.verify_obj(rmfield(expected_mag_str, 'F'), ... + rmfield(mag_str1, 'F')); + end + function test_random_structure_k_and_n(testCase) + swobj = copy(testCase.swobj); + k = [0; 0; 1/4]; + n = [1 1 0]; + swobj.genmagstr('mode','random', 'k', k', 'n', n); + mag_str1 = swobj.mag_str; + % Check F size and magnitude + testCase.verifySize(mag_str1.F, [3 1]); + testCase.verify_val(vecnorm(real(swobj.mag_str.F), 2), 1); + % Check imaginary component of F is perpendicular to n + testCase.verify_val(dot(imag(swobj.mag_str.F), n), 0); + % Check other fields + expected_mag_str = testCase.default_mag_str; + expected_mag_str.k = k; + testCase.verify_obj(rmfield(expected_mag_str, 'F'), ... + rmfield(mag_str1, 'F')); + end + function test_random_structure_multiatom_and_nExt(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 2); + nExt = int32([2 2 2]); + swobj.genmagstr('mode', 'random', 'nExt', nExt); + mag_str1 = swobj.mag_str; + swobj.genmagstr('mode', 'random', 'nExt', nExt); + mag_str2 = swobj.mag_str; + % Check structure is random each time - F is different + testCase.verifyNotEqual(mag_str1.F, mag_str2.F); + % Check F size and magnitude + testCase.verifySize(mag_str1.F, [3 16]); + testCase.verify_val( ... + vecnorm(real(swobj.mag_str.F(:, 1:2:end)), 2), ones(1, 8)); + testCase.verify_val( ... + vecnorm(real(swobj.mag_str.F(:, 2:2:end)), 2), 2*ones(1, 8)); + % Check imaginary component of F is perpendicular to default n + testCase.verifyEqual( ... + dot(imag(swobj.mag_str.F), repmat([0; 0; 1], 1, 16)), zeros(1, 16)); + % Check other fields + expected_mag_str = testCase.default_mag_str; + expected_mag_str.nExt = nExt; + testCase.verify_obj(rmfield(expected_mag_str, 'F'), ... + rmfield(mag_str1, 'F')); + end + function test_tile_existing_struct_extend_cell(testCase) + % Test that tile and increasing nExt will correctly tile + % initialised structure + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + nExt = int32([1 2 1]); + % Also test if we input 'k' it is set to 0 in final struct + swobj.genmagstr('mode', 'direct', 'S', [1 0; 0 1; 0 0], 'k', [1/2 0 0]); + swobj.genmagstr('mode', 'tile', 'nExt', nExt); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.nExt = nExt; + expected_mag_str.F = [1 0 1 0; 0 1 0 1; 0 0 0 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_tile_existing_struct_same_size(testCase) + % Test that tile with nExt same as initialised structure + % does nothing + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + nExt = int32([1 2 1]); + S = [1 0 0 -1; 0 1 0 0; 0 0 1 0]; + swobj.genmagstr('mode', 'direct', 'S', S, 'nExt', nExt); + swobj.genmagstr('mode', 'tile', 'nExt', nExt); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.nExt = nExt; + expected_mag_str.F = S; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_tile_input_S_extend_cell(testCase) + % Test that tile and input S less than nExt will correctly tile + % input S + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + nExt = int32([3 1 1]); + swobj.genmagstr('mode', 'tile', 'nExt', nExt, ... + 'S', [1 0; 0 1; 0 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.nExt = nExt; + expected_mag_str.F = [1 0 1 0 1 0; 0 1 0 1 0 1; 0 0 0 0 0 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_tile_multik(testCase) + % Test that S is summed over third dimension with tile, and k + % is not needed (is this the behaviour we want?) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + S = cat(3, [1 0; 0 1; 0 0], [0 1; 0 0; 1 0]); + swobj.genmagstr('mode', 'tile', ... + 'S', S); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = sqrt(2)/2*[1 1; 0 1; 1 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_tile_multik_provided_k_set_to_zero(testCase) + % Test that S is summed over third dimension with tile, and if + % k is provided, it is set to zero anyway + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + S = cat(3, [1 0; 0 1; 0 0], [0 1; 0 0; 1 0]); + k = 0.5*ones(size(S, 3), 3); + swobj.genmagstr('mode', 'tile', ... + 'S', S, 'k', k); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = sqrt(2)/2*[1 1; 0 1; 1 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_extend_mode_input_S_extend_cell_and_warns(testCase) + % Test undocumented 'extend' mode does same as tile + % Test that tile and input S less than nExt will correctly tile + % input S + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + nExt = int32([3 1 1]); + testCase.verifyWarning(... + @() swobj.genmagstr('mode', 'extend', 'nExt', nExt, ... + 'S', [1 0; 0 1; 0 0]), ... + 'spinw:genmagstr:DeprecationWarning'); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.nExt = nExt; + expected_mag_str.F = [1 0 1 0 1 0; 0 1 0 1 0 1; 0 0 0 0 0 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_rotate_phi(testCase, rotate_phi_inputs) + swobj = copy(testCase.swobj); + k = [1/2 0 0]; + % Need to initialise structure before rotating it + swobj.genmagstr('mode', 'direct', 'S', [1; 0; 0], 'k', k); + swobj.genmagstr('mode', 'rotate', rotate_phi_inputs{:}); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = sqrt(2)/2*[1; 1; 0]; + expected_mag_str.k = k'; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_rotate_multiatom_n(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + k = [1/2 0 0]; + swobj.genmagstr('mode', 'direct', 'S', [1 0; 0 1; 0 0], 'k', k); + swobj.genmagstr('mode', 'rotate', 'phi', pi/2, 'n', [1 1 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [0.5 0.5; 0.5 0.5; -sqrt(2)/2 sqrt(2)/2]; + expected_mag_str.k = k'; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_rotate_no_phi_collinear(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + swobj.genmagstr('mode', 'direct', 'S', [1 -1; 0 0; 0 0]); + swobj.genmagstr('mode', 'rotate', 'n', [0 1 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [0 0; 1 -1; 0 0]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_rotate_no_phi_coplanar(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + swobj.genmagstr('mode', 'direct', 'S', [1 0; 0 1; 0 0]); + swobj.genmagstr('mode', 'rotate', 'n', [0 1 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [1 0; 0 0; 0 -1]; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_rotate_no_phi_incomm(testCase) + swobj_tri = copy(testCase.swobj_tri); + k = [1/3 1/3 0]; + swobj_tri.genmagstr('mode', 'helical', 'S', [1; 0; 0], ... + 'k', k); + swobj_tri.genmagstr('mode', 'rotate', 'n', [1 0 0]); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [0; 1.5i; -1.5]; + expected_mag_str.k = k'; + testCase.verify_obj(expected_mag_str, swobj_tri.mag_str); + end + function test_func_multiatom_default(testCase) + swobj = copy(testCase.swobj); + swobj.addatom('r', [0.5 0.5 0], 'S', 1); + k = [1/3 0 0]; + x0 = [pi/2 -pi/4 0 0 k pi pi/2]; + swobj.genmagstr('mode', 'func', 'x0', x0); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [sqrt(2)/2*(1-i) 0; -sqrt(2)/2*(1+i) 0; 0 1]; + expected_mag_str.k = k'; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + function test_func_custom(testCase) + function [S, k, n] = func(S0, x0) + S = [-S0; 0; 0]; + k = x0; + n = x0; + end + swobj = copy(testCase.swobj); + x0 = [1/3 0 0]; + swobj.genmagstr('mode', 'func', 'func', @func, 'x0', x0); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = [-1; 0; 0]; + expected_mag_str.k = x0'; + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + end + methods (Test, TestTags = {'Symbolic'}) + function test_func_custom_symbolic(testCase) + function [S, k, n] = func(S0, x0) + S = [-S0; 0; 0]; + k = x0; + n = x0; + end + swobj = copy(testCase.swobj); + swobj.symbolic(true); + x0 = [1/3 0 0]; + swobj.genmagstr('mode', 'func', 'func', @func, 'x0', x0); + expected_mag_str = testCase.default_mag_str; + expected_mag_str.F = sym([-1; 0; 0]); + expected_mag_str.k = sym(x0'); + testCase.verify_obj(expected_mag_str, swobj.mag_str); + end + end +end \ No newline at end of file diff --git a/CHANGELOG.rst b/CHANGELOG.rst index be97f4ed..9fdae6f0 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -37,4 +37,31 @@ Bug Fixes basis vectors supplied to ``genlattice``, see issue `#57 `_ - Validation added for ``perm` and ``origin`` arguments supplied to ``genlattice`` (and warn users that these will be ignored if no symmetry/spacegroup is supplied in the same function call). -- Deprecated ``spgr`` argument to ``genlattice`` (users should use ``sym`` instead). \ No newline at end of file +- Deprecated ``spgr`` argument to ``genlattice`` (users should use ``sym`` instead). +- Fix ``MATLAB:nonLogicalConditional`` error raised when using multiple + k in ``genmagstr`` with ``helical`` mode +- Raise error if invalid shape ``S`` or ``k`` is provided to ``genmagstr``, + previously they would be silently set to zero +- Raise error if wrong number of spins ``S`` is provided to ``genmagstr`` in + ``helical`` mode. Previously the structure would be silently initialised + to a random structure. +- Raise error if a complex spin ``S`` is provided to ``genmagstr`` in + ``helical`` mode. Previously this meant it would silently ignore the + ``n`` option, and behave exactly like ``fourier`` mode. +- Raise error if ``rotate`` mode is used without first initialising + a magnetic structure +- Emit deprecation warning if the undocumented ``extend`` mode is used + in ``genmagstr`` +- Raise error if the first spin is parallel to ``n`` and no rotation + angle is provided in ``rotate`` mode in ``genmagstr``. Previously + this would silently result in ``NaN`` +- Raise error if ``phi`` or ``phid`` is not real in ``rotate`` mode in + ``genmagstr``. This was an undocumented feature which has been removed. +- Emit warning that the spin amplitude will be moderated if components + of ``S`` are parallel to ``n`` in ``helical`` mode in ``genmagstr`` +- Emit warning if ``nExt`` is unnecessarily large compared to ``k`` in + ``helical`` and ``fourier`` modes in ``genmagstr`` +- Emit warning if arguments that will be ignored are passed to a particular + mode in ``genmagstr`` (e.g. ``S`` is passed to ``random``) +- Raise error if complex values is provided for ``n`` in ``genmagstr``. + Previously this would've caused a crash. diff --git a/swfiles/@spinw/genmagstr.m b/swfiles/@spinw/genmagstr.m index fdc1af68..02eec46c 100644 --- a/swfiles/@spinw/genmagstr.m +++ b/swfiles/@spinw/genmagstr.m @@ -50,12 +50,18 @@ function genmagstr(obj, varargin) % % `'mode'` % : Mode that determines how the magnetic structure is generated: -% * `'random'` (reads -) -% generates random zero-k magnetic structure. -% * `'direct'` (reads `S`, `n`, `k`) +% * `'random'` (optionally reads `k`, `n`, `nExt`) +% generates a random structure in the structural cell if no other +% arguments are specified here or previously in this spinw +% object. If `nExt` is specified all spins in the supercell are +% randomised. If `k` is specified a random helical structure with +% moments perpendicular to `n` (default value: `[0 0 1]`) with +% the specified `k` propagation vector is generated. (`n` is not +% otherwise used). +% * `'direct'` (reads `S`, optionally reads `k`, `nExt`) % direct input of the magnetic structure using the % parameters of the single-k magnetic structure. -% * `'tile'` (reads `S`, `n`, `k`) +% * `'tile'` (reads `S`, optionally reads `nExt`) % Simply extends the existing or input structure % (`S`) into a magnetic supercell by replicating it. % If no structure is stored in the [spinw] object a random @@ -67,7 +73,7 @@ function genmagstr(obj, varargin) % Magnetic ordering wavevector `k` will be set to zero. To % generate structure with non-zero k, use the `'helical'` or % `'direct'` option. -% * `'helical'` (reads `S`, `n`, `k`, `r0`) +% * `'helical'` (reads `S`, optionally reads `n`, `k`, `r0`, `nExt`, `epsilon`) % generates helical structure in a single cell or in a % supercell. In contrary to the `'extend'` option the % magnetic structure is not generated by replication but @@ -80,7 +86,7 @@ function genmagstr(obj, varargin) % cell. In the first case $r$ denotes the atomic % positions, while for the second case $r$ denotes the % position of the origin of the cell. -% * `'rotate'` (reads `S`, `n`, `k`) +% * `'rotate'` (optionally reads `S`, `phi`, `phid`, `n`, `nExt`) % uniform rotation of all magnetic moments with a % `phi` angle around the given `n` vector. If % `phi=0`, all moments are rotated so, that the first @@ -102,7 +108,7 @@ function genmagstr(obj, varargin) % vector. The default value for `func` is `@gm_spherical3d`. For planar % magnetic structure use `@gm_planar`. Only `func` and `x` % have to be defined for this mode. -% * `'fourier'` (reads `S`, `n`, `k`, `r0`) +% * `'fourier'` (reads `S`, optionally reads `k`, `r0`, `nExt`, `epsilon`) % same as `'helical'`, except the `S` option is taken as the % Fourier components, thus if it contains real numbers, it will % generate a sinusoidally modulated structure instead of @@ -132,8 +138,9 @@ function genmagstr(obj, varargin) % % `'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`. +% structures, stored in a 3-element row vector, is automatically +% normalised to a unit vector. Default value `[0 0 1]`. The coordinate +% system of the vector is determined by `unit`. % % `'S'` % : Vector values of the spins (expectation value), dimensions are $[3\times n_{spin} n_K]$. @@ -220,6 +227,36 @@ function genmagstr(obj, varargin) param = sw_readparam(inpForm, varargin{:}); +% Error if S or k is provided but is empty. This means it has failed the +% input validation, but hasn't caused an error because soft=True +err_str = []; +if isstruct(varargin{1}) + varargin_struct = varargin{1}; +else + varargin_struct = struct(varargin{:}); +end +varargin_names = fields(varargin_struct); +if any(strcmpi(varargin_names, 'S')) && isempty(param.S) + err_str = ["S"]; +end +if any(strcmpi(varargin_names, 'k')) && isempty(param.k) + err_str = [err_str "k"]; +end +if length(err_str) > 0 + error('spinw:genmagstr:WrongInput', 'Incorrect input size for ' + ... + join(err_str, ', ')); +end + +if strcmpi(param.mode, 'rotate') && isempty(obj.mag_str.F) + error('spinw:genmagstr:WrongInput', ['rotate mode requires a magnetic ' ... + 'structure to be defined with another mode first']) +end + +% Complex n causes a crash, error out instead +if ~isreal(param.n) + error('spinw:genmagstr:WrongInput', 'n should be real') +end + if isempty(param.k) noK = true; param.k = k0; @@ -231,12 +268,39 @@ function genmagstr(obj, varargin) error('spinw:genmagstr:WrongInput','''nExt'' has to be larger than 0!'); end +% input type for S, check whether it is complex type +cmplxS = ~isreal(param.S); +if strcmpi(param.mode, 'helical') && cmplxS + error('spinw:genmagstr:WrongInput', ... + ['S must be real for helical mode. To specify complex basis ' ... + 'vectors directly use fourier mode.']) +end + if strcmp(param.mode,'extend') + warning('spinw:genmagstr:DeprecationWarning',... + 'extend mode is deprecated, please use tile mode instead'); param.mode = 'tile'; end -% input type for S, check whether it is complex type -cmplxS = ~isreal(param.S); +mode_args = struct("random", ["k", "n", "nExt", "unit"], ... + "direct", ["S", "k", "nExt", "unit"], ... + "tile", ["S", "nExt", "unit"], ... + "helical", ["S", "n", "k", "r0", "nExt", "epsilon", "unit"], ... + "rotate", ["S", "phi", "phid", "n", "nExt", "unit"], ... + "func", ["func", "x0"], ... + "fourier", ["S", "k", "r0", "nExt", "epsilon", "unit"]); +if ~any(strcmp(fields(mode_args), param.mode)) + error('spinw:genmagstr:WrongMode','Wrong param.mode value!'); +else + unused_args = {varargin_names{:}}; + unused_args(ismember(lower(unused_args), ["mode" "norm"])) = []; + unused_args(ismember(lower(unused_args), lower(mode_args.(lower(param.mode))))) = []; + if ~isempty(unused_args) + warning('spinw:genmagstr:UnreadInput', ... + 'Input(s) %s will be ignored in %s mode', ... + string(join(unused_args', ', ')), param.mode) + end +end switch lower(param.unit) case 'lu' @@ -284,10 +348,6 @@ function genmagstr(obj, varargin) % number of magnetic atoms in the supercell nMagExt = nMagAtom*prod(nExt); -if nMagAtom==0 - error('spinw:genmagstr:NoMagAtom','There is no magnetic atom in the unit cell with S>0!'); -end - % Create mAtom.Sext matrix. mAtom = sw_extendlattice(nExt, mAtom); @@ -298,16 +358,6 @@ function genmagstr(obj, varargin) end n = bsxfunsym(@rdivide,param.n,sqrt(sum(param.n.^2,2))); -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 - -% if the magnetic structure is not initialized start with a random real one -if strcmp(param.mode,'tile') && (nMagAtom > size(param.S,2)) - param.mode = 'random'; -end - % convert input into symbolic variables if obj.symb k = sym(k); @@ -315,7 +365,7 @@ function genmagstr(obj, varargin) n = sym(n); end -if ~cmplxS && ~strcmpi(param.mode,'fourier') && ~strcmpi(param.mode,'direct') && any(k) +if ~cmplxS && ~strcmpi(param.mode,'fourier') && ~strcmpi(param.mode,'direct') && any(k(:)) param.S = param.S + 1i*cross(repmat(permute(n,[2 3 1]),[1 size(param.S,2) 1]),param.S); end @@ -327,9 +377,11 @@ function genmagstr(obj, varargin) % cells defined in obj % -the number of spins stored in obj is not equal to the number % of spins in the final structure - if any(double(obj.mag_str.nExt) - double(param.nExt)) || (size(param.S,2) ~= nMagExt) - S = param.S(:,1:nMagAtom,:); - S = repmat(S,[1 prod(nExt) 1]); + if nMagAtom ~= size(param.S,2) && nMagExt ~= size(param.S,2) + error('spinw:genmagstr:WrongInput', ['Incorrect input size for S, ' ... + 'S must be provided for each magnetic atom']); + elseif any(double(obj.mag_str.nExt) - double(param.nExt)) || (size(param.S,2) ~= nMagExt) + S = repmat(param.S,[1 prod(nExt) 1]); else S = param.S; end @@ -352,11 +404,18 @@ function genmagstr(obj, varargin) S0 = param.S; % Magnetic ordering wavevector in the extended unit cell. kExt = k.*nExt; - % Warns about the non sufficient extension of the unit cell. - % we substitute random values for symbolic km + % Substitute random values for symbolic km skExt = sw_sub1(kExt,'rand'); if any(abs(skExt(:)-round(skExt(:)))>param.epsilon) && prod(nExt) > 1 + % Warns about the non sufficient extension of the unit cell. warning('spinw:genmagstr:UCExtNonSuff','In the extended unit cell k is still larger than epsilon!'); + else + idx = find(sum(k, 1) == 0); + if any(round(skExt(:)) > 1) || any(nExt(idx) > 1) + % Warns about overextension of the unit cell. + warning('spinw:genmagstr:UCExtOver', ... + 'Unit cell has been extended beyond what is required for the given k') + end end % number of spins in the input nSpin = size(param.S,2); @@ -395,6 +454,17 @@ function genmagstr(obj, varargin) else S = S0; end + if strcmpi(param.mode, 'helical') + for ik=1:size(k, 1) + Sk = real(param.S(:, :, ik)); + if any(dot(repmat(n(ik, :), size(Sk, 2), 1)', Sk)) + warning('spinw:genmagstr:SnParallel', ... + ['There are spin components parallel to n, the ' ... + 'amplitude of these components will be modulated']); + break; + end + end + end case 'direct' % direct input of real magnetic moments S = param.S; @@ -417,10 +487,10 @@ function genmagstr(obj, varargin) S = param.S; if ~isreal(param.phi) - % rotate the first spin along [100] - S1 = S(:,1)-sum(n*S(:,1))*n'; - S1 = S1/norm(S1); - param.phi = -atan2(cross(n,[1 0 0])*S1,[1 0 0]*S1); + error('spinw:genmagstr:ComplexPhi', ... + ['Phi should be real! If you really intended phi to have ' ... + 'a complex component, this undocumented feature has been ' ... + 'removed. Please contact the spinw developers.']) end if param.phi == 0 @@ -435,6 +505,11 @@ function genmagstr(obj, varargin) % Axis of rotation defined by the spin direction nRot = cross(n,S1); + if all(nRot == 0) + error('spinw:genmagstr:InvalidRotation', ... + ['Cannot automatically determine rotation as the ' ... + 'first spin is parallel to n']); + end % Angle of rotation. phi = -atan2(norm(cross(S1,n)),dot(S1,n)); else @@ -459,9 +534,6 @@ function genmagstr(obj, varargin) if any(k) S = S + 1i*cross(repmat(permute(n,[2 3 1]),[1 size(S,2) 1]),S); end - - otherwise - error('spinw:genmagstr:WrongMode','Wrong param.mode value!'); end % normalize the magnetic moments