diff --git a/swfiles/@spinw/brille_init.m b/swfiles/@spinw/brille_init.m new file mode 100644 index 000000000..1a9ec909e --- /dev/null +++ b/swfiles/@spinw/brille_init.m @@ -0,0 +1,239 @@ +function brille_init(obj, varargin) +% Initialises and constructs a brille grid for interpolation +% +% ### Syntax +% +% `brille_init(obj, varargin)` +% +% ### Description +% +% This function sets up a Brille BZTrellisQ grid for interpolating the magnon +% energies (eigenvalues) and either the spin-spin correlation function Sab or +% the eigenvectors of the spin Hamiltonian. It then uses `spinwave` to +% calculate these quantities at specified points in the grid (`fill` the grid). +% +% This function is meant to be called from `spinwave` when the option +% `use_brille` is set. +% +% There is a check that if the grid has been filled and the state of the spinw +% object has not changed then it will not be refilled. +% +% ### Input Arguments +% +% `k` +% : A propagation vector - if specified will override the propagation vector +% stored in the .mag_str field of the spinW object. +% `nExt` +% : The super cell size - if specified will override the size in .mag_str. +% +% In order to generate the correct first (irreducible) _magnetic_ Brillouin zone, +% we need to use either the magnetic propagation vector or supercell size. +% Normally, the values in the .mag_str field is used but this can be overriden +% by the user if it turns out that the interpolation is incorrect. +% +% `use_primitive` +% : If set to true, brille will find the (irreducible) first Brillouin zone +% of the primitive rather than conventional unit cell. default is true +% +% `wedge_search` +% : If set to true, brille will find the irreducible first Brillouin zone, +% otherwise just the first BZ will be used. default is true +% +% `search_length` +% : An integer to control how-far the vertex-finding algorithm should +% search in τ-index. The default (==1) indicates that (1̄1̄1̄), (1̄1̄0), (1̄1̄1), +% (1̄0̄1), ..., (111) are included. +% +% `node_volume_fraction` +% : The fractional volume of a single tetrahedron in the grid. Smaller numbers +% will result in more accurate interpolation at the cost of greate computing +% time. default is 1e-5 +% +% `use_vectors` +% : If true will interpolate the eigenvectors of the Hamiltonian rather than +% the spin-spin correlation functions Sab. default is false +% +% In addition, all parameters accepted by spinwave is accepted here and will +% be passed to `spinwave` to `fill` the grid. +% +% ### Output Arguments +% +% None. The function sets the hidden .brille field of the spinW object which +% is then used by `spinwave` to interpolate. +% +% ### See Also +% +% [spinw.spinwave] +% +% (C) 2020 Greg Tucker and Duc Le + +% TODO: + +inpForm.fname = {'k' 'nExt' 'use_primitive' 'search_length'}; +inpForm.defval = {NaN*[0;0;0] NaN*[0;0;0] true 1 }; +inpForm.size = {[1 -1] [1 -1] [1 1] [1 1] }; + +inpForm.fname = [inpForm.fname {'wedge_search' 'node_volume_fraction'}]; +inpForm.defval = [inpForm.defval {true 1e-5 }]; +inpForm.size = [inpForm.size {[1 1] [1 1] }]; + +inpForm.fname = [inpForm.fname {'use_vectors' 'fid' 'hermit' 'Qtrans'}]; +inpForm.defval = [inpForm.defval {false -1 true NaN }]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] [-1 -2] }]; + +[kwds, passthrough] = sw_readparam(inpForm, varargin{:}); +passthrough = [passthrough {'hermit' kwds.hermit}]; + +pref = swpref; +if kwds.fid == -1 + kwds.fid = pref.fid; +end +fid = kwds.fid; + +% Check if we have done this before +obj_to_hash = struct(obj); +obj_to_hash.kwds = kwds; +obj_to_hash.passthrough = passthrough; +new_hash = get_hash(obj_to_hash); +if isfield(obj.brille, 'hash') && all(obj.brille.hash == new_hash) + return; +end + +% In SpinW notation, nExt and k define separate functionalities and need not be related. +% SpinW uses nExt to define the supercell and k to define a single-k magnetic structure +% which can be represented by a single rotating coordinate frame (so not all single-k +% structures can be represented). It is possible in SpinW to have both a supercell +% and a single-k structure by specifying non-unity values for nExt and non-zero +% values for k. In this case the spins of adjacent supercells will have their +% angles rotated by the amount determined by k. + +user_specified_nExt_or_k = true; +if all(isnan(kwds.k)) && all(isnan(kwds.nExt)) + % For incommensurate structure, should use the full (structural) BZ + % rather than the reduced BZ corresponding to the magnetic k-vector + % So here we ignore the SpinW k-vector and just use the supercell size. + nExt = double(obj.mag_str.nExt); + k = 1 ./ nExt; + user_specified_nExt_or_k = false; + % But for incommensurate structures we should use the full rather than + % irreducible BZ. + [~, kd] = rat(obj.mag_str.k(:), 1e-5); + if any(kd ~= 1) + kwds.wedge_search = false; + end +elseif ~all(isnan(kwds.k)) && ~all(isnan(kwds.nExt)) + % User specified nExt and k + k = kwds.k; + nExt = kwds.nExt; +elseif ~all(isnan(kwds.k)) + % User specified k only + warning('Using user specified k, ignoring SpinW nExt and k values'); + k = kwds.k; + [~, nExt] = rat(k(:), 1e-5); +else + % User specified nExt only + warning('Using user specified nExt, ignoring SpinW nExt and k values'); + nExt = kwds.nExt; + k = 1 ./ nExt; +end + +assert(numel(nExt)==3, 'the number of unit cell extensions, nExt, must be (1,3) or (3,1)') +nExt = nExt(:); + +if numel(k)==3 + % nExt and k should be compatible, with nExt a direct lattice "vector" and + % k a reciprocal lattice vector + [kn, kd] = rat(k(:), 1e-5); + % the rationalized denominator of k should (normally) be nExt + if sum(abs(kd - nExt)) > sum(abs(kd + nExt))*eps() + warning('k=(%d/%d,%d/%d,%d/%d) and nExt=(%d,%d,%d) are not compatible',... + kn(1), kd(1), kn(2), kd(2), kn(3), kd(3), ... + nExt(1), nExt(2), nExt(3)) + nExt = max( kd, nExt); + warning('Returning supercell lattice for nExt=(%d,%d,%d)',nExt(1),nExt(2),nExt(3)) + end +end + +lens = obj.lattice.lat_const(:); +angs = obj.lattice.angle(:); % SpinW stores angles in radian +spg = obj.lattice.label; + +if ~isnan(kwds.Qtrans) + obj.brille.Qtrans = kwds.Qtrans; +elseif kwds.use_vectors || user_specified_nExt_or_k + obj.brille.Qtrans = diag(nExt); + lens = lens .* nExt; +else + obj.brille.Qtrans = eye(3); + kwds.node_volume_fraction = kwds.node_volume_fraction * prod(nExt); +end + +% Construct the brille grids and fill it. +fprintf0(fid, 'Creating Brille grid\n'); +obj.brille.bz = brille.create_bz(lens, angs, spg, ... + 'use_primitive', kwds.use_primitive, ... + 'search_length', kwds.search_length, ... + 'time_reversal_symmetry', false, ... + 'wedge_search', kwds.wedge_search); +obj.brille.grid = brille.create_grid(obj.brille.bz, ... + 'node_volume_fraction', kwds.node_volume_fraction, ... + 'complex_vectors', kwds.use_vectors, ... + 'complex_values', ~kwds.hermit); +hkl = permute(light_python_wrapper.p2m(obj.brille.grid.rlu), [2 1]); +if sum(sum(abs(obj.brille.Qtrans - eye(3)))) > 1e-6 + hkl = obj.brille.Qtrans \ hkl; + kwds.node_volume_fraction = kwds.node_volume_fraction * prod(diag(obj.brille.Qtrans)); +end +fprintf0(fid, 'Filling Brille grid\n'); +if kwds.use_vectors + spec = obj.spinwave(hkl, passthrough{:}, 'saveV', true, 'sortMode', false, 'optmem', 0); + [omega, V] = parse_twin(spec); + if (size(omega, 1) / size(V, 1)) == 3 && (size(V, 3) / size(omega, 2)) == 3 + % Incommensurate + kmIdx = repmat(sort(repmat([1 2 3],1,size(omega, 2))),1,1); + eigvec = permute(cat(1, V(:,:,kmIdx==1), V(:,:,kmIdx==2), V(:,:,kmIdx==3)), [3 1 2]); + else + eigvec = permute(V, [3 1 2]); + end +else + spec = obj.spinwave(hkl, passthrough{:}, 'sortMode', false, 'formfact', false, 'optmem', 0); + [omega, Sab] = parse_twin(spec, 'Sab'); + eigvec = permute(real(Sab), [4 3 1 2]); +end +eigval = permute(real(omega), [2 1]); +eigvec_sz = size(eigvec); +eigvec_span = prod(eigvec_sz(3:end)); +obj.brille.grid.fill(eigval, 1, eigvec, eigvec_span); + +% Update hash so we don't refill it unecessarily +obj.brille.hash = new_hash; + +end + +function out = get_hash(obj) + % Calculates a hash for an object or struct using undocumented built-ins + % Based on DataHash (https://uk.mathworks.com/matlabcentral/fileexchange/31272-datahash) + Engine = java.security.MessageDigest.getInstance('MD5'); + Engine.update(getByteStreamFromArray(obj)); + out = typecast(Engine.digest, 'uint8'); +end + +function [omega, V] = parse_twin(spec, use_Sab) + if iscell(spec.omega) + % Has twins + omega = spec.omega{1}; + if nargin > 1 && any(use_Sab) + V = spec.Sab{1}; + else + V = spec.V{1}; + end + else + omega = spec.omega; + if nargin > 1 && all(logical(use_Sab)) + V = spec.Sab; + else + V = spec.V; + end + end +end + diff --git a/swfiles/@spinw/spinw.m b/swfiles/@spinw/spinw.m index beedd7630..f8022adf3 100644 --- a/swfiles/@spinw/spinw.m +++ b/swfiles/@spinw/spinw.m @@ -427,6 +427,8 @@ % use the version property as contant, this will be executed only % once ver = sw_version; + % stores brille objects + brille = struct(); end methods (Static) @@ -709,4 +711,4 @@ function notify(varargin) end end % classdef -end \ No newline at end of file +end diff --git a/swfiles/@spinw/spinwave.m b/swfiles/@spinw/spinwave.m index 4824b1657..b5a541b91 100644 --- a/swfiles/@spinw/spinwave.m +++ b/swfiles/@spinw/spinwave.m @@ -298,7 +298,11 @@ inpForm.defval = [inpForm.defval {false -1 -1 }]; inpForm.size = [inpForm.size {[1 1] [1 1] [1 1] }]; -param = sw_readparam(inpForm, varargin{:}); +inpForm.fname = [inpForm.fname {'use_brille' 'use_vectors'}]; +inpForm.defval = [inpForm.defval {false false}]; +inpForm.size = [inpForm.size {[1 1] [1 -1]}]; + +[param, ~] = sw_readparam(inpForm, varargin{:}); if ~param.fitmode % save the time of the beginning of the calculation @@ -365,7 +369,7 @@ hkl = {hkl}; end -if incomm +if incomm && ~param.use_brille % TODO if ~helical && ~param.fitmode warning('spinw:spinwave:Twokm',['The two times the magnetic ordering '... @@ -516,7 +520,7 @@ % magnetic couplings, 3x3xnJ JJ = cat(3,reshape(SS.all(6:14,:),3,3,[]),SI.aniso); -if incomm +if incomm && ~param.use_brille % transform JJ due to the incommensurate wavevector [~, K] = sw_rot(n,km*dR*2*pi); % multiply JJ with K matrices for every interaction @@ -603,6 +607,15 @@ end +% If we're using Brille, initialise it +if param.use_brille + idx = 2 * find(cellfun(@(c) strcmp(c, 'use_brille'), varargin(1:2:end))); + pars = varargin([1:(idx-2) (idx+1):end]); + obj.brille_init(pars{:}); + if incomm + nHkl0 = nHkl0 * 3; + end +end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -699,6 +712,22 @@ twinIdxMEM = twinIdx(hklIdxMEM); nHklMEM = size(hklExtMEM,2); + if param.use_brille && ~param.use_vectors + [eigval, eigvec] = obj.brille.grid.ir_interpolate_at(obj.brille.Qtrans * hkl(:,hklIdxMEM)); + if incomm + kmIdx = repmat(sort(repmat([1 2 3],1,size(eigval,2)/3)),1,1); + eigval = permute(eigval, [2 1]); + omega = cat(2, omega, eigval(kmIdx==1,:), eigval(kmIdx==2,:), eigval(kmIdx==3,:)); + eigvec = permute(eigvec, [3 4 2 1]); + Sab = cat(4, Sab, eigvec(:,:,kmIdx==1,:), eigvec(:,:,kmIdx==2,:), eigvec(:,:,kmIdx==3,:)); + else + omega = cat(2, omega, permute(eigval, [2 1])); + Sab = cat(4, Sab, permute(eigvec, [3 4 2 1])); + end + sw_timeit(jj/nSlice*100,0,param.tid); + continue; + end + % Creates the matrix of exponential factors nCoupling x nHkl size. % Extends dR into 3 x 3 x nCoupling x nHkl % ExpF = exp(1i*permute(sum(repmat(dR,[1 1 nHklMEM]).*repmat(... @@ -774,7 +803,21 @@ gComm = diag(gCommd); %gd = diag(g); - if param.hermit + if param.use_brille && param.use_vectors + [eigval, V] = obj.brille.grid.ir_interpolate_at(obj.brille.Qtrans * hkl(:,hklIdxMEM)); + V = permute(V, [2 3 1]); + if incomm + kmIdx = repmat(sort(repmat([1 2 3],1,size(V,2))),1,1); + V = cat(3, V(kmIdx==1,:,:), V(kmIdx==2,:,:), V(kmIdx==3,:,:)); + eigval = permute(eigval, [2 1]); + omega = cat(2, omega, eigval(kmIdx==1,:), eigval(kmIdx==2,:), eigval(kmIdx==3,:)); + hklExt0MEM = repmat(hklExt0MEM, [1 3]); + nHklMEM = nHklMEM * 3; + hklIdxMEM = repmat(hklIdxMEM, [1 3]); + else + omega = cat(2, omega, permute(eigval, [2 1])); + end + elseif param.hermit % All the matrix calculations are according to Colpa's paper % J.H.P. Colpa, Physica 93A (1978) 327-353 @@ -904,7 +947,7 @@ % Sab(alpha,beta,iMode,iHkl), size: 3 x 3 x 2*nMagExt x nHkl. % Normalizes the intensity to single unit cell. Sab = cat(4,Sab,squeeze(sum(zeda.*ExpFL.*VExtL,4)).*squeeze(sum(zedb.*ExpFR.*VExtR,3))/prod(nExt)); - + sw_timeit(jj/nSlice*100,0,param.tid); end @@ -939,7 +982,15 @@ if incomm % resize matrices due to the incommensurability (k-km,k,k+km) multiplicity - kmIdx = repmat(sort(repmat([1 2 3],1,nHkl0/3)),1,nTwin); + if param.use_brille + kmIdx = []; + for jj = 1:nSlice + hklIdxMEM = hklIdx(jj):(hklIdx(jj+1)-1); + kmIdx = cat(2, kmIdx, repmat(sort(repmat([1 2 3],1,size(hklIdxMEM,2))),1,nTwin)); + end + else + kmIdx = repmat(sort(repmat([1 2 3],1,nHkl0/3)),1,nTwin); + end % Rodrigues' rotation formula. nx = [0 -n(3) n(2); n(3) 0 -n(1); -n(2) n(1) 0]; nxn = n'*n; @@ -973,12 +1024,19 @@ Sab = cat(3,mmat(Sab(:,:,:,kmIdx==1),K1), mmat(Sab(:,:,:,kmIdx==2),K2), ... mmat(Sab(:,:,:,kmIdx==3),conj(K1))); - hkl = hkl(:,kmIdx==2); + if ~param.use_brille + hkl = hkl(:,kmIdx==2); + end nHkl0 = nHkl0/3; else helical = false; end +if param.use_brille && ~param.use_vectors && param.formfact + % Assume that all ions are the same so just use first row of FF + Sab = Sab .* repmat(permute(FF(1,:), [1 3 4 2]), [3 3 size(Sab,3)]); +end + if ~param.notwin % Rotate the calculated correlation function into the twin coordinate % system using rotC @@ -1072,4 +1130,4 @@ %fprintf(repmat('\b',[1 30])); end -end \ No newline at end of file +end diff --git a/swfiles/sw_readparam.m b/swfiles/sw_readparam.m index 67b4ed60e..a661ba5bb 100644 --- a/swfiles/sw_readparam.m +++ b/swfiles/sw_readparam.m @@ -1,4 +1,4 @@ -function input = sw_readparam(format, varargin) +function [input, remainder] = sw_readparam(format, varargin) % parse input arguments % % ### Syntax @@ -132,15 +132,23 @@ if ~all(usedField) mName = rName(~usedField); - if numel(mName) == 1 - wName = [' "' mName{1} '"']; - elseif numel(mName) == 2 - wName = ['s "' mName{1} '" and "' mName{2} '"']; + if nargout > 1 + remainder = {}; + for ii = 1:numel(mName) + remainder = [remainder {mName{ii} raw.(mName{ii})}]; + end else - wName = ['s ' sprintf('"%s", "',mName{1:(end-2)}) mName{end-1} '" and "' mName{end} '"']; + if numel(mName) == 1 + wName = [' "' mName{1} '"']; + elseif numel(mName) == 2 + wName = ['s "' mName{1} '" and "' mName{2} '"']; + else + wName = ['s ' sprintf('"%s", "',mName{1:(end-2)}) mName{end-1} '" and "' mName{end} '"']; + end + warning('sw_readparam:UnreadInput','Unregognised input parameter%s!',wName); end - - warning('sw_readparam:UnreadInput','Unregognised input parameter%s!',wName); +elseif nargout > 1 + remainder = {}; end end