Skip to content

Commit

Permalink
bug fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@stestoll
stestoll committed Mar 26, 2024
1 parent f9fecc9 commit 0434f90
Show file tree
Hide file tree
Showing 4 changed files with 62 additions and 39 deletions.
2 changes: 1 addition & 1 deletion easyspin/resfields.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@
error('Options.Threshold must be a number >=0 and <1.');
end
preSelectionThreshold = Opt.Threshold(1);
if numel(Opt.Threshold)==1
if isscalar(Opt.Threshold)
postSelectionThreshold = preSelectionThreshold;
else
postSelectionThreshold = Opt.Threshold(2);
Expand Down
38 changes: 25 additions & 13 deletions easyspin/resfields_perturb.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@
end

% Determine excitation mode
[xi1,xik,nB1,nk,nB0,mwmode] = p_excitationgeometry(Exp.mwMode);
[xi1,xik,nB1,nk,nB0_L,mwmode] = p_excitationgeometry(Exp.mwMode);

% Temperature, non-equilibrium populations
if isfield(Exp,'Temperature')
Expand Down Expand Up @@ -306,8 +306,8 @@

% Loop over all orientations
for iOri = nOrientations:-1:1
R = erot(Orientations(iOri,:));
n0 = R.'*nB0; % transform to molecular frame representation
R_L2M = erot(Orientations(iOri,:)).'; % lab frame -> molecular frame
n0 = R_L2M*nB0_L; % transform to molecular frame representation
vecs(:,iOri) = n0;

geff(iOri) = norm(g.'*n0);
Expand Down Expand Up @@ -346,13 +346,13 @@
end
else
if mwmode.linearpolarizedMode
nB1_ = R.'*nB1; % transform to molecular frame representation
nB1_ = R_L2M*nB1; % transform to molecular frame representation
TransitionRate(:,iOri) = c2*norm(cross(g.'*nB1_,u))^2;
elseif mwmode.unpolarizedMode
nk_ = R.'*nk; % transform to molecular frame representation
nk_ = R_L2M*nk; % transform to molecular frame representation
TransitionRate(:,iOri) = c2/2*(trgg-norm(g*u)^2-norm(cross(g.'*nk_,u))^2);
elseif mwmode.circpolarizedMode
nk_ = R.'*nk; % transform to molecular frame representation
nk_ = R_L2M*nk; % transform to molecular frame representation
TransitionRate(:,iOri) = c2*(trgg-norm(g*u)^2-norm(cross(g.'*nk_,u))^2) + ...
mwmode.circSense*2*c2*det(g)*xik/norm(g.'*n0);
end
Expand All @@ -376,7 +376,7 @@
for mS = S:-1:-S+1
imS = imS + 1;

% first-order
% first-order correction
if highSpin
E1D = -uDu/2*(3-6*mS);
else
Expand All @@ -398,7 +398,7 @@
E1A = 0;
end

% second-order
% second-order correction
E2D = 0;
if secondOrder
if highSpin
Expand Down Expand Up @@ -469,6 +469,7 @@
B = [];
Int = [];
Wid = [];
Transitions = [];
spec = spec/dB/prod(nNucStates);
spec = spec*(2*pi); % powder chi integral
else
Expand All @@ -481,6 +482,7 @@
%-------------------------------------------------------------------
nNucSublevels = prod(nNucStates);
Int = repelem(Intensity,nNucSublevels,1)/nNucSublevels;
Int = flipud(Int);

% Widths
%-------------------------------------------------------------------
Expand All @@ -495,9 +497,9 @@
if any(Sys.gStrain(:)) || any(Sys.AStrain(:))

if any(Sys.gStrain(:))
gStrainMatrix = diag(Sys.gStrain(1,:)./Sys.g(1,:))*Exp.mwFreq*1e3; % -> MHz
gStrainMatrix = diag(Sys.gStrain(1,:)./Sys.g(1,:))*Exp.mwFreq*1e3; % -> MHz
if any(Sys.gFrame(:))
R_g2M = erot(Sys.gFrame(1,:)).'; % g frame -> molecular frame
R_g2M = erot(Sys.gFrame(1,:)).'; % g frame -> molecular frame
gStrainMatrix = R_g2M*gStrainMatrix*R_g2M.';
end
else
Expand Down Expand Up @@ -571,8 +573,20 @@
Wid2_DE = repmat(lwD.^2+lwE.^2,nNucTrans,1);
Wid = sqrt(Wid2_DE + Wid.^2);
end


% Transitions
%-------------------------------------------------------------------
Transitions = [];
lowerLevels = (1:nNucSublevels).';
for mSidx = 1:2*S
upperLevels = lowerLevels + nNucSublevels; % only correct for weak HFC
newTransitions = [lowerLevels upperLevels];
Transitions = [Transitions; newTransitions]; %#ok
lowerLevels = upperLevels;
end

spec = 0;

end

% Reshape arrays in the case of crystals with site splitting
Expand All @@ -585,8 +599,6 @@
if ~isempty(Wid), Wid = reshape(Wid,siz); end
end

Transitions = NaN; % not implemented

% Arrange output
%---------------------------------------------------------------
Output = {B,Int,Wid,Transitions,spec};
Expand Down
47 changes: 29 additions & 18 deletions easyspin/resfreqs_perturb.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@
end
end

for iNuc = 1:nNuclei
for iNuc = nNuclei:-1:1
if Sys.fullA
A{iNuc} = Sys.A((iNuc-1)*3+(1:3),:);
else
Expand Down Expand Up @@ -163,7 +163,7 @@
err = '';
if ~isfield(Exp,'Field'), err = 'Exp.Field is missing.'; end

[xi1,xik,nB1,nk,nB0,mwmode] = p_excitationgeometry(Exp.mwMode);
[xi1,xik,nB1,nk,nB0_L,mwmode] = p_excitationgeometry(Exp.mwMode);

if isfield(Exp,'Temperature')
if numel(Exp.Temperature)>1
Expand Down Expand Up @@ -242,7 +242,7 @@

B0 = Exp.Field*1e-3; % mT -> T

for iNuc = 1:nNuclei
for iNuc = nNuclei:-1:1
A_ = A{iNuc};
detA(iNuc) = det(A_);
invA{iNuc} = inv(A_); % gives an error with zero hf couplings
Expand Down Expand Up @@ -285,16 +285,16 @@
gg = g*g.';
trgg = trace(gg);

% prefactor for transition rate
% Prefactor for transition rate
mS = S:-1:-S+1;
c = bmagn/2 * sqrt(S*(S+1)-mS.*(mS-1));
c = c/planck/1e9;
c2 = c.^2;

% Loop over all orientations
for iOri = nOrientations:-1:1
R = erot(Orientations(iOri,:));
n0 = R.'*nB0; % transform to molecular frame representation
R_L2M = erot(Orientations(iOri,:)).'; % lab frame -> molecular frame
n0 = R_L2M*nB0_L; % transform to molecular frame representation
vecs(:,iOri) = n0;

geff(iOri) = norm(g.'*n0);
Expand Down Expand Up @@ -332,13 +332,13 @@
end
else
if mwmode.linearpolarizedMode
nB1_ = R.'*nB1; % transform to molecular frame representation
nB1_ = R_L2M*nB1; % transform to molecular frame representation
TransitionRate(:,iOri) = c2*norm(cross(g.'*nB1_,u))^2;
elseif mwmode.unpolarizedMode
nk_ = R.'*nk; % transform to molecular frame representation
nk_ = R_L2M*nk; % transform to molecular frame representation
TransitionRate(:,iOri) = c2/2*(trgg-norm(g*u)^2-norm(cross(g.'*nk_,u))^2);
elseif mwmode.circpolarizedMode
nk_ = R.'*nk; % transform to molecular frame representation
nk_ = R_L2M*nk; % transform to molecular frame representation
TransitionRate(:,iOri) = c2*(trgg-norm(g*u)^2-norm(cross(g.'*nk_,u))^2) + ...
circSense*2*c2*det(g)*xik/norm(g.'*n0);
end
Expand Down Expand Up @@ -451,6 +451,7 @@
nu = [];
Int = [];
Wid = [];
Transitions = [];
spec = spec/dnu/prod(nNucStates);
spec = spec*(2*pi); % powder chi integral
else
Expand All @@ -462,6 +463,7 @@
%-------------------------------------------------------------------
nNucSublevels = prod(nNucStates);
Int = repelem(Intensity,nNucSublevels,1)/nNucSublevels;
Int = flipud(Int);

% Widths
%-------------------------------------------------------------------
Expand All @@ -476,9 +478,9 @@
if any(Sys.gStrain(:)) || any(Sys.AStrain(:))

if any(Sys.gStrain(:))
gStrainMatrix = diag(Sys.gStrain(1,:)./Sys.g(1,:))*E0*1e3; % -> MHz
gStrainMatrix = diag(Sys.gStrain(1,:)./Sys.g(1,:))*E0*1e3; % -> MHz
if any(Sys.gFrame(:))
R_g2M = erot(Sys.gFrame(1,:)).'; % g frame -> molecular frame
R_g2M = erot(Sys.gFrame(1,:)).'; % g frame -> molecular frame
gStrainMatrix = R_g2M*gStrainMatrix*R_g2M.';
end
else
Expand Down Expand Up @@ -551,26 +553,35 @@
else
Wid = [];
end

% Transitions
%-------------------------------------------------------------------
Transitions = [];
lowerLevels = (1:nNucSublevels).';
for mSidx = 1:2*S
upperLevels = lowerLevels + nNucSublevels; % only correct for weak HFC
newTransitions = [lowerLevels upperLevels];
Transitions = [Transitions; newTransitions]; %#ok
lowerLevels = upperLevels;
end

spec = 0;

end

% Reshape arrays in the case of crystals with site splitting
d = dbstack;
pepperCall = numel(d)>2 && strcmp(d(2).name,'pepper');
if (nSites>1) && ~pepperCall
db = dbstack;
pepperCall = numel(db)>2 && strcmp(db(2).name,'pepper');
if nSites>1 && ~pepperCall
siz = [nTransitions*nSites, numel(nu)/nTransitions/nSites];
nu = reshape(nu,siz);
if ~isempty(Int), Int = reshape(Int,siz); end
if ~isempty(Wid), Wid = reshape(Wid,siz); end
end

Transitions = NaN; % not implemented

% Arrange output
%---------------------------------------------------------------
Output = {nu,Int,Wid,Transitions,spec};
varargout = Output(1:max(nargout,1));

return

end
14 changes: 7 additions & 7 deletions tests/pepper_pt_temp.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,24 +3,24 @@
% temperature effects perturbation <-> matrix diagonalization

Sys.S = 1;
Sys.g = [2];
Sys.g = 2;
Sys.D = 400;
%Sys.A = [50 80];
%Sys.Nucs = '1H';
Sys.lwpp = 0.5;

Exp.mwFreq = 9.5;
Exp.CenterSweep = [330 60];
Exp.Temperature = 0.5;
Exp.Harmonic = 0;

Opt.Method = 'matrix';
[x,y0]=pepper(Sys,Exp,Opt);
[B,spc0] = pepper(Sys,Exp,Opt);
Opt.Method = 'perturb';
[x,y1]=pepper(Sys,Exp,Opt);
[~,spc1] = pepper(Sys,Exp,Opt);

if opt.Display
plot(x,y0,x,y1);
plot(B,spc0,B,spc1);
legend('matrix','perturb');
legend boxoff
end

ok = areequal(y0,y1,2e-2,'rel');
ok = areequal(spc0,spc1,2e-2,'rel');

0 comments on commit 0434f90

Please sign in to comment.