Skip to content

Commit 2434933

Browse files
tsdevtsdev
committed
Implementing the rlu transformation
- minor bug fixes - newcell() and spinwave() are now compatible with the rlu transformation
1 parent c3d48a7 commit 2434933

File tree

5 files changed

+163
-120
lines changed

5 files changed

+163
-120
lines changed

‎swfiles/+swplot/plotmag.m‎

Lines changed: 93 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
% 'circle' Plot only the rotation plane of incommensurate
2929
% magnetic structures.
3030
% 'arrow' Plots only the moment directions.
31+
% 'none' Don't plot anything.
3132
% figure Handle of the swplot figure. Default is the selected figure.
3233
% legend Whether to add the plot to the legend, default is true.
3334
% label Whether to plot labels for atoms, default is true.
@@ -207,98 +208,105 @@
207208
% do nothing
208209
case 'circle'
209210
case 'arrow'
211+
case 'none'
210212
otherwise
211213
error('plotmag:WrongInput','The given mode string is invalid!');
212214
end
213215

214-
% number of unit cells
215-
nCell = prod(nExtPlot);
216-
217-
% keep track of types of atoms
218-
aIdx = repmat(mAtom.idx,[nCell 1])';
219-
a2Idx = repmat(1:numel(mAtom.idx),[nCell 1])';
220-
pos = reshape(pos,3,[]);
221-
222-
% cut out the atoms that are out of range
223-
switch param.unit
224-
case 'lu'
225-
% L>= lower range, L<= upper range
226-
pIdx = all(bsxfun(@ge,pos,range(:,1)-10*eps) & bsxfun(@le,pos,range(:,2)+10*eps),1);
227-
case 'xyz'
228-
% convert to xyz
229-
posxyz = BV*pos;
230-
pIdx = all(bsxfun(@ge,posxyz,range(:,1)-10*eps) & bsxfun(@le,posxyz,range(:,2)+10*eps),1);
231-
end
232-
233-
if ~any(pIdx)
234-
warning('plotatom:EmptyPlot','There are no magnetic moments in the plotting range!')
235-
return
236-
end
237-
238-
M = M(:,pIdx);
239-
pos = pos(:,pIdx);
240-
aIdx = aIdx(pIdx);
241-
a2Idx = a2Idx(pIdx);
242-
243-
% normalization
244-
if param.normalize
245-
% normalize moments
246-
M = bsxfun(@rdivide,M,sqrt(sum(M.^2,1)));
247-
end
248-
249-
% save magnetic moment vector values into data before rescaling
250-
MDat = mat2cell([M;pos;a2Idx],7,ones(1,size(M,2)));
251-
252-
% scale moments
253-
% get the length of the shortest bond
254-
if ~isempty(obj.coupling.atom2)
255-
apos1 = obj.matom.r(:,obj.coupling.atom1(1));
256-
apos2 = obj.matom.r(:,obj.coupling.atom2(1))+double(obj.coupling.dl(:,1));
257-
lBond = norm(BV*(apos2-apos1));
258-
else
259-
lBond = 3;
260-
end
261-
262-
% normalize the longest moment vector to scale*(shortest bond length)
263-
M = M/sqrt(max(sum(M.^2,1)))*param.scale*lBond;
264-
265-
if param.centered
266-
% double the length for centered moments
267-
M = 2*M;
268-
end
269-
270-
% convert to lu units
271-
Mlu = BV\M;
272-
273-
if param.centered
274-
% center on atoms
275-
vpos = cat(3,pos-Mlu/2,pos+Mlu/2);
276-
else
277-
vpos = cat(3,pos,pos+Mlu);
278-
end
279-
280-
% color
281-
if strcmp(param.color,'auto')
282-
color = double(obj.unit_cell.color(:,aIdx));
216+
if ~strcmp(param.mode,'none')
217+
218+
% number of unit cells
219+
nCell = prod(nExtPlot);
220+
221+
% keep track of types of atoms
222+
aIdx = repmat(mAtom.idx,[nCell 1])';
223+
a2Idx = repmat(1:numel(mAtom.idx),[nCell 1])';
224+
pos = reshape(pos,3,[]);
225+
226+
% cut out the atoms that are out of range
227+
switch param.unit
228+
case 'lu'
229+
% L>= lower range, L<= upper range
230+
pIdx = all(bsxfun(@ge,pos,range(:,1)-10*eps) & bsxfun(@le,pos,range(:,2)+10*eps),1);
231+
case 'xyz'
232+
% convert to xyz
233+
posxyz = BV*pos;
234+
pIdx = all(bsxfun(@ge,posxyz,range(:,1)-10*eps) & bsxfun(@le,posxyz,range(:,2)+10*eps),1);
235+
end
236+
237+
if ~any(pIdx)
238+
warning('plotatom:EmptyPlot','There are no magnetic moments in the plotting range!')
239+
return
240+
end
241+
242+
M = M(:,pIdx);
243+
pos = pos(:,pIdx);
244+
aIdx = aIdx(pIdx);
245+
a2Idx = a2Idx(pIdx);
246+
247+
% normalization
248+
if param.normalize
249+
% normalize moments
250+
M = bsxfun(@rdivide,M,sqrt(sum(M.^2,1)));
251+
end
252+
253+
% save magnetic moment vector values into data before rescaling
254+
MDat = mat2cell([M;pos;a2Idx],7,ones(1,size(M,2)));
255+
256+
% scale moments
257+
% get the length of the shortest bond
258+
if ~isempty(obj.coupling.atom2)
259+
apos1 = obj.matom.r(:,obj.coupling.atom1(1));
260+
apos2 = obj.matom.r(:,obj.coupling.atom2(1))+double(obj.coupling.dl(:,1));
261+
lBond = norm(BV*(apos2-apos1));
262+
else
263+
lBond = 3;
264+
end
265+
266+
% normalize the longest moment vector to scale*(shortest bond length)
267+
M = M/sqrt(max(sum(M.^2,1)))*param.scale*lBond;
268+
269+
if param.centered
270+
% double the length for centered moments
271+
M = 2*M;
272+
end
273+
274+
% convert to lu units
275+
Mlu = BV\M;
276+
277+
if param.centered
278+
% center on atoms
279+
vpos = cat(3,pos-Mlu/2,pos+Mlu/2);
280+
else
281+
vpos = cat(3,pos,pos+Mlu);
282+
end
283+
284+
% color
285+
if strcmp(param.color,'auto')
286+
color = double(obj.unit_cell.color(:,aIdx));
287+
else
288+
color = swplot.color(param.color);
289+
end
290+
291+
% shift positions
292+
vpos = bsxfun(@plus,vpos,BV\param.shift);
293+
294+
% prepare legend labels
295+
mAtom.name = obj.unit_cell.label(mAtom.idx);
296+
lLabel = repmat(mAtom.name,[nCell 1]);
297+
lLabel = lLabel(pIdx);
298+
299+
% plot moment vectors
300+
swplot.plot('type','arrow','name','mag','position',vpos,'text','',...
301+
'figure',hFigure,'legend',false,'color',color,'R',param.radius0,...
302+
'fontsize',param.fontsize,'tooltip',false,'replace',param.replace,...
303+
'data',MDat,'label',lLabel,'nmesh',param.nmesh,'ang',param.ang,...
304+
'lHead',param.lHead,'translate',param.translate,'zoom',param.zoom);
305+
283306
else
284-
color = swplot.color(param.color);
307+
% don't plot anything just remove previous plot
285308
end
286309

287-
% shift positions
288-
vpos = bsxfun(@plus,vpos,BV\param.shift);
289-
290-
% prepare legend labels
291-
mAtom.name = obj.unit_cell.label(mAtom.idx);
292-
lLabel = repmat(mAtom.name,[nCell 1]);
293-
lLabel = lLabel(pIdx);
294-
295-
296-
% plot moment vectors
297-
swplot.plot('type','arrow','name','mag','position',vpos,'text','',...
298-
'figure',hFigure,'legend',false,'color',color,'R',param.radius0,...
299-
'fontsize',param.fontsize,'tooltip',false,'replace',param.replace,...
300-
'data',MDat,'label',lLabel,'nmesh',param.nmesh,'ang',param.ang,...
301-
'lHead',param.lHead,'translate',param.translate,'zoom',param.zoom);
302310

303311
% save range
304312
setappdata(hFigure,'range',struct('range',range,'unit',param.unit));

‎swfiles/@spinw/energy.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,8 @@
113113
dRIdx = find(any(dR));
114114
for ii = 1:numel(dRIdx)
115115
%M2(:,dRIdx(ii)) = sw_rot(n, kExt*dR(:,dRIdx(ii))*2*pi, M2(:,dRIdx(ii)));
116-
M2(:,dRIdx(ii)) = real(bsxfunsym(@times,M2cmplx(:,dRIdx(ii)),exp(1i*kExt*dR(:,dRIdx(ii))*2*pi)));
116+
% TODO check - sign in front of phase
117+
M2(:,dRIdx(ii)) = real(bsxfunsym(@times,M2cmplx(:,dRIdx(ii)),exp(-1i*kExt*dR(:,dRIdx(ii))*2*pi)));
117118
end
118119
end
119120

‎swfiles/@spinw/newcell.m‎

Lines changed: 62 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,35 @@
1-
function varargout = newcell(obj,bvect, bshift)
1+
function varargout = newcell(obj,varargin)
22
% changes lattice vectors while keeping atoms
33
%
4-
% {T} = NEWCELL(obj, bvect, {bshift})
4+
% {T} = NEWCELL(obj, 'option1', value1, ...)
55
%
66
% The function defines new unit cell using the 3 vectors contained in
77
% bvect. The three vectors in lattice units define a parallelepiped. This
88
% will be the new unit cell. The atoms from the original unit cell will
99
% fill the new unit cell. Also the magnetic structure and bond and single
10-
% ion property definitions will be erased from the structure.
10+
% ion property definitions will be erased from the structure. The new cell
11+
% will naturally have different reciprocal lattice, however the original
12+
% reciprocal lattice units can be retained using the 'keepq' option. In
13+
% this case the spinw.spinwave() function will calculate spin wave
14+
% dispersion at reciprocal lattice points of the original lattice. The
15+
% transformation between the two lattices is stored in spinw.unit.qmat.
1116
%
1217
% Input:
1318
%
1419
% obj spinw class object.
20+
%
21+
% Options:
22+
%
1523
% bvect Defines the new lattice vectors in the original lattice
1624
% coordinate system. Cell with the following elements
17-
% {v1 v2 v3}.
18-
% bshift Vector defines a shift of the position of the unit cell.
25+
% {v1 v2 v3} or a 3x3 matrix with v1, v2 and v3 as column
26+
% vectors: [v1 v2 v3].
27+
% bshift Row vector defines a shift of the position of the unit cell.
1928
% Optional.
29+
% keepq If true, the reciprocal lattice units of the new model will be
30+
% the same as in the old model. This is achieved by storing the
31+
% transformation matrix between the new and the old coordinate
32+
% system in spinw.unit.qmat.
2033
%
2134
% Output:
2235
%
@@ -28,14 +41,31 @@
2841
%
2942
% Example:
3043
%
31-
% tri = spinw;
32-
% tri.genlattice('lat_const',[3 3 5],'angled',[90 90 120])
33-
% tri.addatom('r',[0 0 0])
34-
% tri.newcell({[1 0 0] [1 2 0] [0 0 1]})
35-
% plot(tri)
44+
% In this example we generate the triangular lattice antiferromagnet and
45+
% convert the hexagonal cell to orthorhombic. This doubles the number of
46+
% magnetic atoms in the cell and changes the reciprocal lattice. However we
47+
% use 'keepq' to able to index the reciprocal lattice of the orthorhombic
48+
% cell with the reciprocal lattice of the original hexagonal cell. To show
49+
% that the two models are equivalent, we calculate the spin wave spectrum
50+
% on both model using the same rlu. On the orthorhombic cell, the q value
51+
% will be converted automatically and the calculated spectrum will be the
52+
% same for both cases.
53+
%
54+
% tri = sw_model('triAF',1);
55+
% tri_orth = copy(tri);
56+
% tri_orth.newcell('bvect',{[1 0 0] [1 2 0] [0 0 1]},'keepq',true);
57+
% tri_orth.gencoupling
58+
% tri_orth.addcoupling('bond',1,'mat','J1')
59+
% newk = ((tri_orth.unit.qmat)*tri.magstr.k')';
60+
% tri_orth.genmagstr('mode','helical','k',newk,'S',[1 0 0]')
61+
% plot(tri_orth)
62+
%
63+
% figure
64+
% subplot(2,1,1)
65+
% sw_plotspec(sw_egrid(tri.spinwave({[0 0 0] [1 1 0] 501})),'mode','color','dE',0.2)
66+
% subplot(2,1,2)
67+
% sw_plotspec(sw_egrid(tri_orth.spinwave({[0 0 0] [1 1 0] 501})),'mode','color','dE',0.2)
3668
%
37-
% The example show how to convert a triangular lattice into orthorhombic
38-
% lattice vectors and plots the new unit cell.
3969
%
4070
% See also SPINW.GENLATTICE, SPINW.GENCOUPLING, SPINW.NOSYM.
4171
%
@@ -45,18 +75,18 @@
4575
return
4676
end
4777

48-
%warning('spinw:newcell:PossibleBug','There might be an error, if there are atoms at the faces of the original cell!')
78+
inpForm.fname = {'bvect' 'bshift' 'keepq'};
79+
inpForm.defval = {eye(3) [0 0 0] true };
80+
inpForm.size = {[-1 3] [1 3] [1 1] };
81+
82+
param = sw_readparam(inpForm, varargin{:});
4983

50-
if ~iscell(bvect) || numel(bvect)~=3
51-
error('spinw:newcell:WrongInput','Input has to be cell type with 3 vectors inside!');
84+
if ~iscell(param.bvect) && size(param.bvect,1)~=3
85+
error('spinw:newcell:WrongInput','Input has to be 1x3 cell or 3x3 matrix!');
5286
end
5387

5488
% shift
55-
if nargin == 2
56-
bshift = [0;0;0];
57-
else
58-
bshift = bshift(:);
59-
end
89+
bshift = param.bshift(:);
6090

6191
% here 3 coordinate systems are used:
6292
% - xyz real space, Cartesian
@@ -66,8 +96,11 @@
6696

6797
% transformation matrix from the new lattice units into the original
6898
% lattice
69-
% v_orig = Tn_o*v_new
70-
Tn_o = [bvect{1}(:) bvect{2}(:) bvect{3}(:)];
99+
if iscell(param.bvect)
100+
Tn_o = [param.bvect{1}(:) param.bvect{2}(:) param.bvect{3}(:)];
101+
else
102+
Tn_o = param.bvect;
103+
end
71104

72105
% transformation from the original lattice into xyz real space
73106
% xyz = To_xyz * v_orig
@@ -87,8 +120,6 @@
87120

88121
% number of cells needed for the extension
89122
nExt = ceil(max(pp,[],2) - min(pp,[],2))'+2;
90-
%obj.mag_str.N_ext = int32(nExt);
91-
92123

93124
% generated atoms
94125
atomList = obj.atom;
@@ -100,7 +131,6 @@
100131
idxExt = atomList.idxext;
101132

102133

103-
104134
rExt = bsxfun(@plus,rExt,bshift);
105135
% atomic positions in the new unit cell
106136
rNew = inv(Tn_o)*rExt; %#ok<MINV>
@@ -111,7 +141,7 @@
111141
idxCut = any((rNew<-epsilon) | (rNew>(1-epsilon)),1);
112142
rNew(:,idxCut) = [];
113143
idxExt(idxCut) = [];
114-
atomList.Sext(idxCut) = [];
144+
atomList.Sext(idxCut) = [];
115145

116146
% atoms are close to the face or origin --> put them exactly
117147
rNew(rNew<epsilon) = 0;
@@ -135,11 +165,6 @@
135165
obj.unit_cell.(fNames{ii}) = obj.unit_cell.(fNames{ii})(:,idxExt);
136166
end
137167

138-
% reset the magnetic structure
139-
% obj.mag_str.F = zeros(3,0,0);
140-
% obj.mag_str.k = zeros(3,0);
141-
% obj.mag_str.N_ext = int32([1 1 1]);
142-
143168
% reset the magnetic structure and the bonds
144169
obj.initfield({'coupling' 'mag_str'});
145170

@@ -149,7 +174,12 @@
149174
% transformation from the original reciprocal lattice into the new
150175
% reciprocal lattice
151176
if nargout>0
152-
varargout{1} = Tn_o;
177+
varargout{1} = Tn_o';
178+
end
179+
180+
if param.keepq
181+
% keep the new coordinate system
182+
obj.unit.qmat = Tn_o';
153183
end
154184

155185
end

‎swfiles/@spinw/optmagsteep.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -280,7 +280,7 @@
280280
end
281281

282282
if fid == 1
283-
sw_status(0,1);
283+
sw_status(0,1,'Magnetic structure optimization');
284284
end
285285

286286
if nargout == 1

0 commit comments

Comments
 (0)