From 243493372c658d91097bb5bb28122027dc40958d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A1ndor=20T=C3=B3th?= Date: Wed, 12 Apr 2017 15:34:59 +0200 Subject: [PATCH] Implementing the rlu transformation - minor bug fixes - newcell() and spinwave() are now compatible with the rlu transformation --- swfiles/+swplot/plotmag.m | 178 ++++++++++++++++++----------------- swfiles/@spinw/energy.m | 3 +- swfiles/@spinw/newcell.m | 94 +++++++++++------- swfiles/@spinw/optmagsteep.m | 2 +- swfiles/@spinw/spinwave.m | 6 +- 5 files changed, 163 insertions(+), 120 deletions(-) diff --git a/swfiles/+swplot/plotmag.m b/swfiles/+swplot/plotmag.m index 6c710f3b1..e339b4f46 100644 --- a/swfiles/+swplot/plotmag.m +++ b/swfiles/+swplot/plotmag.m @@ -28,6 +28,7 @@ % 'circle' Plot only the rotation plane of incommensurate % magnetic structures. % 'arrow' Plots only the moment directions. +% 'none' Don't plot anything. % figure Handle of the swplot figure. Default is the selected figure. % legend Whether to add the plot to the legend, default is true. % label Whether to plot labels for atoms, default is true. @@ -207,98 +208,105 @@ % do nothing case 'circle' case 'arrow' + case 'none' otherwise error('plotmag:WrongInput','The given mode string is invalid!'); end -% number of unit cells -nCell = prod(nExtPlot); - -% keep track of types of atoms -aIdx = repmat(mAtom.idx,[nCell 1])'; -a2Idx = repmat(1:numel(mAtom.idx),[nCell 1])'; -pos = reshape(pos,3,[]); - -% cut out the atoms that are out of range -switch param.unit - case 'lu' - % L>= lower range, L<= upper range - pIdx = all(bsxfun(@ge,pos,range(:,1)-10*eps) & bsxfun(@le,pos,range(:,2)+10*eps),1); - case 'xyz' - % convert to xyz - posxyz = BV*pos; - pIdx = all(bsxfun(@ge,posxyz,range(:,1)-10*eps) & bsxfun(@le,posxyz,range(:,2)+10*eps),1); -end - -if ~any(pIdx) - warning('plotatom:EmptyPlot','There are no magnetic moments in the plotting range!') - return -end - -M = M(:,pIdx); -pos = pos(:,pIdx); -aIdx = aIdx(pIdx); -a2Idx = a2Idx(pIdx); - -% normalization -if param.normalize - % normalize moments - M = bsxfun(@rdivide,M,sqrt(sum(M.^2,1))); -end - -% save magnetic moment vector values into data before rescaling -MDat = mat2cell([M;pos;a2Idx],7,ones(1,size(M,2))); - -% scale moments -% get the length of the shortest bond -if ~isempty(obj.coupling.atom2) - apos1 = obj.matom.r(:,obj.coupling.atom1(1)); - apos2 = obj.matom.r(:,obj.coupling.atom2(1))+double(obj.coupling.dl(:,1)); - lBond = norm(BV*(apos2-apos1)); -else - lBond = 3; -end - -% normalize the longest moment vector to scale*(shortest bond length) -M = M/sqrt(max(sum(M.^2,1)))*param.scale*lBond; - -if param.centered - % double the length for centered moments - M = 2*M; -end - -% convert to lu units -Mlu = BV\M; - -if param.centered - % center on atoms - vpos = cat(3,pos-Mlu/2,pos+Mlu/2); -else - vpos = cat(3,pos,pos+Mlu); -end - -% color -if strcmp(param.color,'auto') - color = double(obj.unit_cell.color(:,aIdx)); +if ~strcmp(param.mode,'none') + + % number of unit cells + nCell = prod(nExtPlot); + + % keep track of types of atoms + aIdx = repmat(mAtom.idx,[nCell 1])'; + a2Idx = repmat(1:numel(mAtom.idx),[nCell 1])'; + pos = reshape(pos,3,[]); + + % cut out the atoms that are out of range + switch param.unit + case 'lu' + % L>= lower range, L<= upper range + pIdx = all(bsxfun(@ge,pos,range(:,1)-10*eps) & bsxfun(@le,pos,range(:,2)+10*eps),1); + case 'xyz' + % convert to xyz + posxyz = BV*pos; + pIdx = all(bsxfun(@ge,posxyz,range(:,1)-10*eps) & bsxfun(@le,posxyz,range(:,2)+10*eps),1); + end + + if ~any(pIdx) + warning('plotatom:EmptyPlot','There are no magnetic moments in the plotting range!') + return + end + + M = M(:,pIdx); + pos = pos(:,pIdx); + aIdx = aIdx(pIdx); + a2Idx = a2Idx(pIdx); + + % normalization + if param.normalize + % normalize moments + M = bsxfun(@rdivide,M,sqrt(sum(M.^2,1))); + end + + % save magnetic moment vector values into data before rescaling + MDat = mat2cell([M;pos;a2Idx],7,ones(1,size(M,2))); + + % scale moments + % get the length of the shortest bond + if ~isempty(obj.coupling.atom2) + apos1 = obj.matom.r(:,obj.coupling.atom1(1)); + apos2 = obj.matom.r(:,obj.coupling.atom2(1))+double(obj.coupling.dl(:,1)); + lBond = norm(BV*(apos2-apos1)); + else + lBond = 3; + end + + % normalize the longest moment vector to scale*(shortest bond length) + M = M/sqrt(max(sum(M.^2,1)))*param.scale*lBond; + + if param.centered + % double the length for centered moments + M = 2*M; + end + + % convert to lu units + Mlu = BV\M; + + if param.centered + % center on atoms + vpos = cat(3,pos-Mlu/2,pos+Mlu/2); + else + vpos = cat(3,pos,pos+Mlu); + end + + % color + if strcmp(param.color,'auto') + color = double(obj.unit_cell.color(:,aIdx)); + else + color = swplot.color(param.color); + end + + % shift positions + vpos = bsxfun(@plus,vpos,BV\param.shift); + + % prepare legend labels + mAtom.name = obj.unit_cell.label(mAtom.idx); + lLabel = repmat(mAtom.name,[nCell 1]); + lLabel = lLabel(pIdx); + + % plot moment vectors + swplot.plot('type','arrow','name','mag','position',vpos,'text','',... + 'figure',hFigure,'legend',false,'color',color,'R',param.radius0,... + 'fontsize',param.fontsize,'tooltip',false,'replace',param.replace,... + 'data',MDat,'label',lLabel,'nmesh',param.nmesh,'ang',param.ang,... + 'lHead',param.lHead,'translate',param.translate,'zoom',param.zoom); + else - color = swplot.color(param.color); + % don't plot anything just remove previous plot end -% shift positions -vpos = bsxfun(@plus,vpos,BV\param.shift); - -% prepare legend labels -mAtom.name = obj.unit_cell.label(mAtom.idx); -lLabel = repmat(mAtom.name,[nCell 1]); -lLabel = lLabel(pIdx); - - -% plot moment vectors -swplot.plot('type','arrow','name','mag','position',vpos,'text','',... - 'figure',hFigure,'legend',false,'color',color,'R',param.radius0,... - 'fontsize',param.fontsize,'tooltip',false,'replace',param.replace,... - 'data',MDat,'label',lLabel,'nmesh',param.nmesh,'ang',param.ang,... - 'lHead',param.lHead,'translate',param.translate,'zoom',param.zoom); % save range setappdata(hFigure,'range',struct('range',range,'unit',param.unit)); diff --git a/swfiles/@spinw/energy.m b/swfiles/@spinw/energy.m index 0fecc22d0..b9d2fa01e 100644 --- a/swfiles/@spinw/energy.m +++ b/swfiles/@spinw/energy.m @@ -113,7 +113,8 @@ dRIdx = find(any(dR)); for ii = 1:numel(dRIdx) %M2(:,dRIdx(ii)) = sw_rot(n, kExt*dR(:,dRIdx(ii))*2*pi, M2(:,dRIdx(ii))); - M2(:,dRIdx(ii)) = real(bsxfunsym(@times,M2cmplx(:,dRIdx(ii)),exp(1i*kExt*dR(:,dRIdx(ii))*2*pi))); + % TODO check - sign in front of phase + M2(:,dRIdx(ii)) = real(bsxfunsym(@times,M2cmplx(:,dRIdx(ii)),exp(-1i*kExt*dR(:,dRIdx(ii))*2*pi))); end end diff --git a/swfiles/@spinw/newcell.m b/swfiles/@spinw/newcell.m index c80878384..e0d0dc060 100644 --- a/swfiles/@spinw/newcell.m +++ b/swfiles/@spinw/newcell.m @@ -1,22 +1,35 @@ -function varargout = newcell(obj,bvect, bshift) +function varargout = newcell(obj,varargin) % changes lattice vectors while keeping atoms % -% {T} = NEWCELL(obj, bvect, {bshift}) +% {T} = NEWCELL(obj, 'option1', value1, ...) % % The function defines new unit cell using the 3 vectors contained in % bvect. The three vectors in lattice units define a parallelepiped. This % will be the new unit cell. The atoms from the original unit cell will % fill the new unit cell. Also the magnetic structure and bond and single -% ion property definitions will be erased from the structure. +% ion property definitions will be erased from the structure. The new cell +% will naturally have different reciprocal lattice, however the original +% reciprocal lattice units can be retained using the 'keepq' option. In +% this case the spinw.spinwave() function will calculate spin wave +% dispersion at reciprocal lattice points of the original lattice. The +% transformation between the two lattices is stored in spinw.unit.qmat. % % Input: % % obj spinw class object. +% +% Options: +% % bvect Defines the new lattice vectors in the original lattice % coordinate system. Cell with the following elements -% {v1 v2 v3}. -% bshift Vector defines a shift of the position of the unit cell. +% {v1 v2 v3} or a 3x3 matrix with v1, v2 and v3 as column +% vectors: [v1 v2 v3]. +% bshift Row vector defines a shift of the position of the unit cell. % Optional. +% keepq If true, the reciprocal lattice units of the new model will be +% the same as in the old model. This is achieved by storing the +% transformation matrix between the new and the old coordinate +% system in spinw.unit.qmat. % % Output: % @@ -28,14 +41,31 @@ % % Example: % -% tri = spinw; -% tri.genlattice('lat_const',[3 3 5],'angled',[90 90 120]) -% tri.addatom('r',[0 0 0]) -% tri.newcell({[1 0 0] [1 2 0] [0 0 1]}) -% plot(tri) +% In this example we generate the triangular lattice antiferromagnet and +% convert the hexagonal cell to orthorhombic. This doubles the number of +% magnetic atoms in the cell and changes the reciprocal lattice. However we +% use 'keepq' to able to index the reciprocal lattice of the orthorhombic +% cell with the reciprocal lattice of the original hexagonal cell. To show +% that the two models are equivalent, we calculate the spin wave spectrum +% on both model using the same rlu. On the orthorhombic cell, the q value +% will be converted automatically and the calculated spectrum will be the +% same for both cases. +% +% tri = sw_model('triAF',1); +% tri_orth = copy(tri); +% tri_orth.newcell('bvect',{[1 0 0] [1 2 0] [0 0 1]},'keepq',true); +% tri_orth.gencoupling +% tri_orth.addcoupling('bond',1,'mat','J1') +% newk = ((tri_orth.unit.qmat)*tri.magstr.k')'; +% tri_orth.genmagstr('mode','helical','k',newk,'S',[1 0 0]') +% plot(tri_orth) +% +% figure +% subplot(2,1,1) +% sw_plotspec(sw_egrid(tri.spinwave({[0 0 0] [1 1 0] 501})),'mode','color','dE',0.2) +% subplot(2,1,2) +% sw_plotspec(sw_egrid(tri_orth.spinwave({[0 0 0] [1 1 0] 501})),'mode','color','dE',0.2) % -% The example show how to convert a triangular lattice into orthorhombic -% lattice vectors and plots the new unit cell. % % See also SPINW.GENLATTICE, SPINW.GENCOUPLING, SPINW.NOSYM. % @@ -45,18 +75,18 @@ return end -%warning('spinw:newcell:PossibleBug','There might be an error, if there are atoms at the faces of the original cell!') +inpForm.fname = {'bvect' 'bshift' 'keepq'}; +inpForm.defval = {eye(3) [0 0 0] true }; +inpForm.size = {[-1 3] [1 3] [1 1] }; + +param = sw_readparam(inpForm, varargin{:}); -if ~iscell(bvect) || numel(bvect)~=3 - error('spinw:newcell:WrongInput','Input has to be cell type with 3 vectors inside!'); +if ~iscell(param.bvect) && size(param.bvect,1)~=3 + error('spinw:newcell:WrongInput','Input has to be 1x3 cell or 3x3 matrix!'); end % shift -if nargin == 2 - bshift = [0;0;0]; -else - bshift = bshift(:); -end +bshift = param.bshift(:); % here 3 coordinate systems are used: % - xyz real space, Cartesian @@ -66,8 +96,11 @@ % transformation matrix from the new lattice units into the original % lattice -% v_orig = Tn_o*v_new -Tn_o = [bvect{1}(:) bvect{2}(:) bvect{3}(:)]; +if iscell(param.bvect) + Tn_o = [param.bvect{1}(:) param.bvect{2}(:) param.bvect{3}(:)]; +else + Tn_o = param.bvect; +end % transformation from the original lattice into xyz real space % xyz = To_xyz * v_orig @@ -87,8 +120,6 @@ % number of cells needed for the extension nExt = ceil(max(pp,[],2) - min(pp,[],2))'+2; -%obj.mag_str.N_ext = int32(nExt); - % generated atoms atomList = obj.atom; @@ -100,7 +131,6 @@ idxExt = atomList.idxext; - rExt = bsxfun(@plus,rExt,bshift); % atomic positions in the new unit cell rNew = inv(Tn_o)*rExt; %#ok @@ -111,7 +141,7 @@ idxCut = any((rNew<-epsilon) | (rNew>(1-epsilon)),1); rNew(:,idxCut) = []; idxExt(idxCut) = []; -atomList.Sext(idxCut) = []; +atomList.Sext(idxCut) = []; % atoms are close to the face or origin --> put them exactly rNew(rNew0 - varargout{1} = Tn_o; + varargout{1} = Tn_o'; +end + +if param.keepq + % keep the new coordinate system + obj.unit.qmat = Tn_o'; end end \ No newline at end of file diff --git a/swfiles/@spinw/optmagsteep.m b/swfiles/@spinw/optmagsteep.m index 2b0df094a..e9ec76f98 100644 --- a/swfiles/@spinw/optmagsteep.m +++ b/swfiles/@spinw/optmagsteep.m @@ -280,7 +280,7 @@ end if fid == 1 - sw_status(0,1); + sw_status(0,1,'Magnetic structure optimization'); end if nargout == 1 diff --git a/swfiles/@spinw/spinwave.m b/swfiles/@spinw/spinwave.m index 8809b2abe..9909218b2 100644 --- a/swfiles/@spinw/spinwave.m +++ b/swfiles/@spinw/spinwave.m @@ -725,6 +725,10 @@ K = chol(ham(:,:,ii)+eye(2*nMagExt)*param.omega_tol); warn1 = true; catch PD + if param.tid == 2 + % close timer window + sw_status(100,2,param.tid); + end error('spinw:spinwave:NonPosDefHamiltonian',... ['Hamiltonian matrix is not positive definite, probably'... ' the magnetic structure is wrong! For approximate'... @@ -929,7 +933,7 @@ % Creates output structure with the calculated values. spectra.omega = omega; spectra.Sab = Sab; -spectra.hkl = hkl(:,1:nHkl0); +spectra.hkl = obj.unit.qmat\hkl(:,1:nHkl0); spectra.hklA = hklA; spectra.incomm = incomm; spectra.helical = helical;