From 5ef8bf70d5688127ece2c325e93cd164ca578890 Mon Sep 17 00:00:00 2001 From: Stefan Stoll Date: Thu, 16 May 2024 21:55:42 -0700 Subject: [PATCH 1/6] work on strains --- easyspin/private/getdstrainops.m | 72 --------------- easyspin/private/strains_de.m | 90 +++++++++++++++++++ easyspin/private/strains_ga.m | 83 ++++++++++++++++++ easyspin/resfields.m | 140 ++++++++---------------------- easyspin/resfields_perturb.m | 50 ++++++----- easyspin/resfreqs_matrix.m | 54 ++++++------ tests/pepper_dstraincorrelated.m | 33 +++---- tests/resfields_gstrain.m | 16 ++++ tests/resfields_gstrainframe.m | 24 +++++ tests/resfields_pt_gstrain.m | 19 ++++ tests/resfields_pt_gstrainframe.m | 24 +++++ tests/resfreqs_matrix_gstrain.m | 27 +++--- 12 files changed, 385 insertions(+), 247 deletions(-) delete mode 100644 easyspin/private/getdstrainops.m create mode 100644 easyspin/private/strains_de.m create mode 100644 easyspin/private/strains_ga.m create mode 100644 tests/resfields_gstrain.m create mode 100644 tests/resfields_gstrainframe.m create mode 100644 tests/resfields_pt_gstrain.m create mode 100644 tests/resfields_pt_gstrainframe.m diff --git a/easyspin/private/getdstrainops.m b/easyspin/private/getdstrainops.m deleted file mode 100644 index 9df7e00a..00000000 --- a/easyspin/private/getdstrainops.m +++ /dev/null @@ -1,72 +0,0 @@ -% Calculate Hamiltonian derivatives with respect to D and E, taking into -% account correlation if present), premultipy with strains. - -% Sys.DStrain: nElectrons x 2 or nElectrons x 3 array -% with [FWHM_D FWHM_E rDE] on each row. rDE is the correlation coefficient -% between D and E. - -function [useDStrain,dHdD,dHdE] = getdstrainops(Sys) - -useDStrain = any(Sys.DStrain(:)); -dHdD = cell(1,Sys.nElectrons); -dHdE = cell(1,Sys.nElectrons); - -if ~useDStrain - return -end - -% Loop over all electron spins -for iEl = 1:Sys.nElectrons - - % Skip if no D strain is given for this electron spin - if ~any(Sys.DStrain(iEl,1:2)) - dHdD{iEl} = 0; - dHdE{iEl} = 0; - continue - end - - % Construct Zeeman operators in molecular frame - SxM_ = sop(Sys,[iEl,1]); - SyM_ = sop(Sys,[iEl,2]); - SzM_ = sop(Sys,[iEl,3]); - - % Construct Zeeman operators in D frame, and square - if any(Sys.DFrame(iEl,:)) - R_M2D = erot(Sys.DFrame(iEl,:)); % molecular frame -> D frame - SxD2_ = (R_M2D(1,1)*SxM_ + R_M2D(1,2)*SyM_ + R_M2D(1,3)*SzM_)^2; - SyD2_ = (R_M2D(2,1)*SxM_ + R_M2D(2,2)*SyM_ + R_M2D(2,3)*SzM_)^2; - SzD2_ = (R_M2D(3,1)*SxM_ + R_M2D(3,2)*SyM_ + R_M2D(3,3)*SzM_)^2; - else - % D frame aligns with molecular frame - SxD2_ = SxM_^2; - SyD2_ = SyM_^2; - SzD2_ = SzM_^2; - end - - % Calculate derivatives of Hamiltonian w.r.t. D and E - % Spin Hamiltonian terms (in D tensor frame xD, yD, zD) - % H = D*(SzD^2-S(S+1)/3) + E*(SxD^2-SyD^2) - % = D*(2*SzD^2-SxD^2-SyD^2)/3 + E*(SxD^2-SyD^2) - dHdD_ = (2*SzD2_-SxD2_-SyD2_)/3; - dHdE_ = SxD2_-SyD2_; - - % Compute Hamiltonian derivatives, pre-multiply with strain FWHMs. - DeltaD = Sys.DStrain(iEl,1); - DeltaE = Sys.DStrain(iEl,2); - rDE = Sys.DStrainCorr(iEl); % correlation coefficient between D and E - if (rDE~=0) - % Transform correlated D-E strain to uncorrelated coordinates - logmsg(1,' correlated D strain for electron spin %d (r = %f)',iEl,rDE); - % Construct and diagonalize covariance matrix - R12 = rDE*DeltaD*DeltaE; - CovMatrix = [DeltaD^2 R12; R12 DeltaE^2]; - [V,L] = eig(CovMatrix); - L = sqrt(diag(L)); - dHdD{iEl} = L(1)*(V(1,1)*dHdD_ + V(1,2)*dHdE_); - dHdE{iEl} = L(2)*(V(2,1)*dHdD_ + V(2,2)*dHdE_); - else - dHdD{iEl} = DeltaD*dHdD_; - dHdE{iEl} = DeltaE*dHdE_; - end - -end diff --git a/easyspin/private/strains_de.m b/easyspin/private/strains_de.m new file mode 100644 index 00000000..ef4bebd6 --- /dev/null +++ b/easyspin/private/strains_de.m @@ -0,0 +1,90 @@ +% Calculate Hamiltonian derivatives with respect to D and E, taking into +% account correlation if present; premultipy with linewidth. + +% Required input fields +% Sys.DStrain: +% nElectrons x 2 array with [FWHM_D FWHM_E] on each row. +% Sys.DStrainCorr: +% Array with nElectron correlation coefficient between D and E (between -1 and +1). +% Sys.DFrame: +% Euler angles describing D tensor orientation in molecular frame + +function dHdDE = strains_de(Sys) + +% Return if no D/E strain present +if ~any(Sys.DStrain(:)) + dHdDE = {}; + return +end + +% Loop over all electron spins +nElectrons = size(Sys.DStrain,1); +dHdDE = cell(1,2*nElectrons); % preallocate maximal possible size +idx = 0; +for iEl = 1:nElectrons + + % Skip if no D strain is given for this electron spin + if ~any(Sys.DStrain(iEl,1:2)) + continue + end + + % Strain information + Delta = Sys.DStrain(iEl,:); % FWHMs for D and E + rDE = Sys.DStrainCorr(iEl); % correlation coefficient between D and E + if rDE<-1 || rDE>1 + error('Electron spin %d: D/E correlation coefficient must be between -1 and 1',iEl); + end + + % Construct spin operators in molecular frame + SxM = sop(Sys,[iEl,1]); + SyM = sop(Sys,[iEl,2]); + SzM = sop(Sys,[iEl,3]); + + % Transform spin operators to D frame + if any(Sys.DFrame(iEl,:)) + R_M2D = erot(Sys.DFrame(iEl,:)); % molecular frame -> D frame + SxD = R_M2D(1,1)*SxM + R_M2D(1,2)*SyM + R_M2D(1,3)*SzM; + SyD = R_M2D(2,1)*SxM + R_M2D(2,2)*SyM + R_M2D(2,3)*SzM; + SzD = R_M2D(3,1)*SxM + R_M2D(3,2)*SyM + R_M2D(3,3)*SzM; + else + % D frame aligns with molecular frame + SxD = SxM; + SyD = SyM; + SzD = SzM; + end + + % Calculate derivatives of Hamiltonian w.r.t. D and E + % Spin Hamiltonian terms (in D tensor frame xD, yD, zD) + % H = D*(SzD^2-S(S+1)/3) + E*(SxD^2-SyD^2) + % = D*(2*SzD^2-SxD^2-SyD^2)/3 + E*(SxD^2-SyD^2) + dHdD = (2*SzD^2 - SxD^2 - SyD^2)/3; + dHdE = SxD^2 - SyD^2; + + % Diagonalize covariance matrix in the case of non-zero correlation + if rDE~=0 + % Transform correlated D-E strain to uncorrelated coordinates + logmsg(1,' correlated D strain for electron spin %d (r = %f)',iEl,rDE); + % Construct covariance matrix + R12 = rDE*Delta(1)*Delta(2); % covariance + covMatrix = [Delta(1)^2 R12; R12 Delta(2)^2]; + % Diagonalize covariance matrix, get derivatives along principal directions + [V,sigma2] = eig(covMatrix); + Delta = sqrt(diag(sigma2)); + dHdDE1 = Delta(1)*(V(1,1)*dHdD + V(1,2)*dHdE); + dHdDE2 = Delta(2)*(V(2,1)*dHdD + V(2,2)*dHdE); + else + dHdDE1 = Delta(1)*dHdD; + dHdDE2 = Delta(2)*dHdE; + end + + % Store in list + idx = idx+1; + dHdDE{idx} = dHdDE1; + idx = idx+1; + dHdDE{idx} = dHdDE2; + +end + +dHdDE = dHdDE(1:idx); + +end diff --git a/easyspin/private/strains_ga.m b/easyspin/private/strains_ga.m new file mode 100644 index 00000000..ccde0ba5 --- /dev/null +++ b/easyspin/private/strains_ga.m @@ -0,0 +1,83 @@ +% Calculates quantities needed by resfields for g/A strain broadening + +function [usegStrain,useAStrain,simplegStrain,lw2_gA,ops] = strains_ga(CoreSys,Exp,mwFreq,Transitions,nTransitions,kH0,kmuzM) + +usegStrain = any(CoreSys.gStrain(:)); +useAStrain = CoreSys.nNuclei>0 && any(CoreSys.AStrain); + +ops = struct; +lw2_gA = []; + +% g-A strain +%------------------------------------------------- +% g strain tensor is taken to be aligned with the g tensor +% A strain tensor is taken to be aligned with the A tensor +% g strain can be specified for each electron spin +% A strain is limited to the first electron and first nuclear spin +simplegStrain = CoreSys.nElectrons==1; +if usegStrain + logmsg(1,' g strain present'); + lw_g = mwFreq*CoreSys.gStrain./CoreSys.g; % MHz + for e = CoreSys.nElectrons:-1:1 + if any(CoreSys.gFrame(e,:)) + R_g2M{e} = erot(CoreSys.gFrame(e,:)).'; % g frame -> molecular frame + else + R_g2M{e} = eye(3); + end + end + if ~simplegStrain + logmsg(1,' multiple g strains present'); + for e = CoreSys.nElectrons:-1:1 + kSxM{e} = sop(CoreSys,[e,1]); + kSyM{e} = sop(CoreSys,[e,2]); + kSzM{e} = sop(CoreSys,[e,3]); + end + ops.kSxM = kSxM; + ops.kSyM = kSyM; + ops.kSzM = kSzM; + end +else + lw_g = zeros(CoreSys.nElectrons,3); +end + +if useAStrain + if usegStrain + if any(CoreSys.gFrame(1,:)~=CoreSys.AFrame(1,1:3)) + error('For g/A strain, Sys.gFrame and Sys.AFrame must be collinear.'); + end + R_strain2M = R_g2M; + else + R_A2M = erot(CoreSys.AFrame).'; % strain tensor frame -> mol. frame + R_strain2M = R_A2M; + end + lw_A = diag(CoreSys.AStrain); + % Diagonalize Hamiltonian at center field. + centerB = mean(Exp.Range); + [Vecs,E] = eig(kH0 - centerB*kmuzM); + [~,idx] = sort(real(diag(E))); + Vecs = Vecs(:,idx); + % Calculate effective mI of nucleus 1 for all eigenstates. + mI1 = real(diag(Vecs'*sop(CoreSys,[2,3])*Vecs)); + mI1Tr = mean(mI1(Transitions),2); + % compute A strain array + lw_A = reshape(mI1Tr(:,ones(1,9)).',[3,3,nTransitions]).*... + repmat(lw_A,[1,1,nTransitions]); + + % Combine g and A strain + rho = CoreSys.gAStrainCorr; % correlation coefficient + for e = CoreSys.nElectrons:-1:1 + lw2_gA{e} = repmat(diag(lw_g(e,:).^2),[1,1,nTransitions]) + lw_A.^2 + ... + 2*rho*repmat(diag(lw_g(e,:)),[1,1,nTransitions]).*lw_A; + for tr = 1:nTransitions + lw2_gA{e}(:,:,tr) = R_strain2M{e}*lw2_gA{e}(:,:,tr)*R_strain2M{e}.'; + end + end +else + if usegStrain + for e = CoreSys.nElectrons:-1:1 + lw2_gA{e} = repmat(R_g2M{e}*diag(lw_g(e,:).^2)*R_g2M{e}.',[1,1,nTransitions]); + end + end +end +% lw2_gA = a (cell array of) 3D array with 3x3 strain line-width matrices +% for each transition stacked along the third dimension. diff --git a/easyspin/resfields.m b/easyspin/resfields.m index 5ab39998..b08412c3 100644 --- a/easyspin/resfields.m +++ b/easyspin/resfields.m @@ -81,10 +81,9 @@ Sys = adddefaults(Sys,DefaultSys); if numel(Sys.gAStrainCorr)~=1 || ~isnumeric(Sys.gAStrainCorr) || ... - Sys.gAStrainCorr==0 || ~isfinite(Sys.gAStrainCorr) - error('Sys.gAStrainCorr must be a single number, either +1 or -1.'); + Sys.gAStrainCorr<-1 || Sys.gAStrainCorr>1 + error('Sys.gAStrainCorr must be a single number between -1 and +1.'); end -Sys.gAStrainCorr = sign(Sys.gAStrainCorr); if Sys.nElectrons>1 if any(Sys.AStrain(:)) @@ -93,9 +92,7 @@ end if any(Sys.gStrain(:)) || any(Sys.AStrain(:)) - gFull = size(Sys.g,1)==3*numel(Sys.S); - %aFull = size(System.A,1)==3*(1+sum(System.Nucs==',')); - if gFull + if Sys.fullg error('gStrain and AStrain are not supported when full g matrices are given!'); end if any(Sys.DStrain) @@ -267,8 +264,8 @@ end computeFreq2Field = Opt.Freq2Field; -StrainsPresent = any([Sys.HStrain(:); Sys.DStrain(:); Sys.gStrain(:); Sys.AStrain(:)]); -computeStrains = StrainsPresent && (nargout>2); +strainsPresent = any([Sys.HStrain(:); Sys.DStrain(:); Sys.gStrain(:); Sys.AStrain(:)]); +computeStrains = strainsPresent && nargout>2; computeGradient = (computeStrains || (nargout>4)) && GradientSwitch; computeIntensities = (nargout>1 && IntensitySwitch) || computeGradient; @@ -381,10 +378,10 @@ end % Components of S vectors for computing - for iEl = Sys.nElectrons:-1:1 - S(iEl).x = sop(CoreSys,[iEl,1]); - S(iEl).y = sop(CoreSys,[iEl,2]); - S(iEl).z = sop(CoreSys,[iEl,3]); + for e = Sys.nElectrons:-1:1 + S(e).x = sop(CoreSys,[e,1]); + S(e).y = sop(CoreSys,[e,2]); + S(e).z = sop(CoreSys,[e,3]); end else @@ -723,90 +720,30 @@ % Line width preparations %======================================================================= logmsg(1,'- Broadenings'); -simplegStrain = true; -usegStrain = false; -useAStrain = false; + if computeStrains logmsg(1,' using strains'); % Frequency-domain residual width tensor - %----------------------------------------------- HStrain2 = CoreSys.HStrain.^2; - % D strain - %----------------------------------------------- - [useDStrain,dHdD,dHdE] = getdstrainops(CoreSys); - - % g-A strain - %------------------------------------------------- - % g strain tensor is taken to be aligned with the g tensor - % A strain tensor is taken to be aligned with the A tensor - % g strain can be specified for each electron spin - % A strain is limited to the first electron and first nuclear spin - usegStrain = any(CoreSys.gStrain(:)); - if usegStrain - logmsg(1,' g strain present'); - simplegStrain = CoreSys.nElectrons==1; - for iEl = CoreSys.nElectrons:-1:1 - gStrainMatrix{iEl} = diag(CoreSys.gStrain(iEl,:)./CoreSys.g(iEl,:))*mwFreq; % MHz - if any(CoreSys.gFrame(iEl,:)) - R_g2M = erot(CoreSys.gFrame(iEl,:)).'; % g frame -> molecular frame - gStrainMatrix{iEl} = R_g2M*gStrainMatrix{iEl}*R_g2M.'; - end - end - if ~simplegStrain - logmsg(1,' multiple g strains present'); - for iEl = CoreSys.nElectrons:-1:1 - kSxM{iEl} = sop(CoreSys,[iEl,1]); - kSyM{iEl} = sop(CoreSys,[iEl,2]); - kSzM{iEl} = sop(CoreSys,[iEl,3]); - end - end - else - for e = CoreSys.nElectrons:-1:1 - gStrainMatrix{e} = zeros(3); - end - end - - useAStrain = (CoreSys.nNuclei>0) && any(CoreSys.AStrain); - if useAStrain - % Transform A strain matrix to molecular frame. - AStrainMatrix = diag(CoreSys.AStrain); - if isfield(CoreSys,'AFrame') - R_A2M = erot(CoreSys.AFrame(1,:)).'; % A frame -> molecular frame - AStrainMatrix = R_A2M*AStrainMatrix*R_A2M.'; - end - % Diagonalize Hamiltonian at center field. - centerB = mean(Exp.Range); - [Vecs,E] = eig(kH0 - centerB*kmuzM); - [~,idx] = sort(real(diag(E))); - Vecs = Vecs(:,idx); - % Calculate effective mI of nucleus 1 for all eigenstates. - mI = real(diag(Vecs'*sop(CoreSys,[2,3])*Vecs)); - mITr = mean(mI(Transitions),2); - % compute A strain array - AStrainMatrix = reshape(mITr(:,ones(1,9)).',[3,3,nTransitions]).*... - repmat(AStrainMatrix,[1,1,nTransitions]); - corr = Sys.gAStrainCorr; - for e = Sys.nElectrons:-1:1 - gAslw2{e} = (repmat(gStrainMatrix{e},[1,1,nTransitions])+corr*AStrainMatrix).^2; - end - clear AStrainMatrix Vecs E idx mI mITr - else - for e = Sys.nElectrons:-1:1 - gAslw2{e} = repmat(gStrainMatrix{e}.^2,[1,1,nTransitions]); - end - end - clear gslw - % gAslw2 = a (cell array of) 3D array with 3x3 strain line-width matrices - % for each transition piled up along the third dimension. + % D/E strain + sigmadHdDE = strains_de(CoreSys); + useDStrain = ~isempty(sigmadHdDE); + % g/A strain + [usegStrain,useAStrain,simplegStrain,lw2_gA,ops] = strains_ga(CoreSys,Exp,mwFreq,Transitions,nTransitions,kH0,kmuzM); + if any(HStrain2), logmsg(2,' ## using H strain'); end if usegStrain || useAStrain, logmsg(2,' ## using g/A strain'); end if useDStrain, logmsg(2,' ## using D strain'); end else - logmsg(1,' no strains specified',nTransitions); + logmsg(1,' no strains specified'); + + usegStrain = false; + useAStrain = false; + useDStrain = false; end @@ -929,7 +866,7 @@ kmuyL = yLab_M(1)*kmuxM + yLab_M(2)*kmuyM + yLab_M(3)*kmuzM; if usegStrain && ~simplegStrain for e = Sys.nElectrons:-1:1 - kSzL{e} = zLab_M(1)*kSxM{e} + zLab_M(2)*kSyM{e} + zLab_M(3)*kSzM{e}; + kSzL{e} = zLab_M(1)*ops.kSxM{e} + zLab_M(2)*ops.kSyM{e} + zLab_M(3)*ops.kSzM{e}; end end end @@ -1287,31 +1224,30 @@ % Calculate width if requested %-------------------------------------------------- if computeStrains - %m = @(Op) real(V'*Op*V) - real(U'*Op*U); - m = @(Op) real((V'-U')*Op*(V+U)); % equivalent to prev. line + deltaVU = @(Op) real(V'*Op*V - U'*Op*U); % H strain LineWidth2 = LineWidthSquared; - % D strain + % D/E strain if useDStrain - for iEl = 1:CoreSys.nElectrons - LineWidth2 = LineWidth2 + abs(m(dHdD{iEl}))^2; - LineWidth2 = LineWidth2 + abs(m(dHdE{iEl}))^2; + for i = 1:numel(sigmadHdDE) + LineWidth2 = LineWidth2 + deltaVU(sigmadHdDE{i})^2; end end % g and A strain if usegStrain || useAStrain if simplegStrain - gA2 = gAslw2{1}(:,:,iTrans); + LineWidth2_gA = lw2_gA{1}(:,:,iTrans); else - gA2 = 0; - for iEl = 1:Sys.nElectrons - gA2 = gA2 + abs(m(kSzL{iEl}))*gAslw2{iEl}(:,:,iTrans); + LineWidth2_gA = 0; + for e = 1:Sys.nElectrons + delta_mS = deltaVU(kSzL{e}); + LineWidth2_gA = LineWidth2_gA + abs(delta_mS)*lw2_gA{e}(:,:,iTrans); end end - LineWidth2 = LineWidth2 + zLab_M.'*gA2*zLab_M; + LineWidth2 = LineWidth2 + zLab_M.'*LineWidth2_gA*zLab_M; end % Convert to field value and save @@ -1325,18 +1261,18 @@ %------------------------------------------------------- if nPerturbNuclei>0 % Compute S vector expectation values for all electron spins - for iEl = Sys.nElectrons:-1:1 - Su(:,iEl) = [U'*S(iEl).x*U; U'*S(iEl).y*U; U'*S(iEl).z*U]; - Sv(:,iEl) = [V'*S(iEl).x*V; V'*S(iEl).y*V; V'*S(iEl).z*V]; + for e = Sys.nElectrons:-1:1 + Su(:,e) = [U'*S(e).x*U; U'*S(e).y*U; U'*S(e).z*U]; + Sv(:,e) = [V'*S(e).x*V; V'*S(e).y*V; V'*S(e).z*V]; end % Build and diagonalize nuclear sub-Hamiltonians for iiNuc = 1:nPerturbNuclei Hu = 0; Hv = 0; % Hyperfine (dependent on S) - for iEl = 1:Sys.nElectrons - Hu = Hu + Su(1,iEl)*Hhfi(iEl,iiNuc).x + Su(2,iEl)*Hhfi(iEl,iiNuc).y + Su(3,iEl)*Hhfi(iEl,iiNuc).z; - Hv = Hv + Sv(1,iEl)*Hhfi(iEl,iiNuc).x + Sv(2,iEl)*Hhfi(iEl,iiNuc).y + Sv(3,iEl)*Hhfi(iEl,iiNuc).z; + for e = 1:Sys.nElectrons + Hu = Hu + Su(1,e)*Hhfi(e,iiNuc).x + Su(2,e)*Hhfi(e,iiNuc).y + Su(3,e)*Hhfi(e,iiNuc).z; + Hv = Hv + Sv(1,e)*Hhfi(e,iiNuc).x + Sv(2,e)*Hhfi(e,iiNuc).y + Sv(3,e)*Hhfi(e,iiNuc).z; end % Nuclear Zeeman and quadrupole (independent of S) if ~Opt.HybridOnlyHFI diff --git a/easyspin/resfields_perturb.m b/easyspin/resfields_perturb.m index 269a8d6a..625fa265 100644 --- a/easyspin/resfields_perturb.m +++ b/easyspin/resfields_perturb.m @@ -471,7 +471,7 @@ if immediateBinning B = []; Int = []; - Wid = []; + lw = []; Transitions = []; spec = spec/dB/prod(nNucStates); spec = spec*(2*pi); % powder chi integral @@ -490,49 +490,53 @@ % Widths %------------------------------------------------------------------- if any(Sys.HStrain) - lw2 = sum(Sys.HStrain.^2*vecs.^2,1); % MHz^2 - lw = sqrt(lw2)*1e6*planck./geff/bmagn*1e3; % mT - Wid = repmat(lw,nNucTrans*2*S,1); + lw2 = sum(Sys.HStrain.^2*vecs.^2,1); % MHz^2 + lw = sqrt(lw2)*1e6*planck./geff/bmagn*1e3; % MHz -> mT + lw = repmat(lw,nNucTrans*2*S,1); else - Wid = 0; + lw = 0; end if any(Sys.gStrain(:)) || any(Sys.AStrain(:)) if any(Sys.gStrain(:)) - gStrainMatrix = diag(Sys.gStrain(1,:)./Sys.g(1,:))*Exp.mwFreq*1e3; % -> MHz + lw_g = 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 - gStrainMatrix = R_g2M*gStrainMatrix*R_g2M.'; + else + R_g2M = eye(3); end else - gStrainMatrix = zeros(3); + lw_g = zeros(3); end if any(Sys.AStrain) && Sys.nNuclei>0 - AStrainMatrix = diag(Sys.AStrain); + lw_A = diag(Sys.AStrain); if isfield(Sys,'AFrame') - Rp = erot(Sys.AFrame(1,:)).'; % A frame -> molecular frame - AStrainMatrix = Rp*AStrainMatrix*Rp.'; + if any(Sys.gFrame~=Sys.AFrame(1,:)) + error('For g/A strain, the g and A tensors must be collinear.'); + end end - corr = Sys.gAStrainCorr; + rho = Sys.gAStrainCorr; % correlation coefficient mI1 = -Sys.I(1):+Sys.I(1); for idx = 1:numel(mI1) - StrainMatrix = gStrainMatrix + corr*(mI1(idx))*AStrainMatrix; + lw_A_ = mI1(idx).*lw_A; + lw2_gA_ = lw_g.^2 + lw_A_.^2 + 2*rho*lw_g.*lw_A_; + lw2_gA_ = R_g2M*lw2_gA_*R_g2M.'; for iOri = 1:nOrientations - lw2(idx,iOri) = vecs(:,iOri).'*StrainMatrix.^2*vecs(:,iOri); + lw2_gA(idx,iOri) = vecs(:,iOri).'*lw2_gA_*vecs(:,iOri); end end - Wid_gA = sqrt(lw2)*planck*1e6./repmat(geff,numel(mI1),1)/bmagn*1e3; % MHz -> mT + whos lw2 geff mI1 + lw_gA = sqrt(lw2_gA)*planck*1e6./repmat(geff,numel(mI1),1)/bmagn*1e3; % MHz -> mT idx = repmat(1:numel(mI1),2*S*nNucTrans/numel(mI1),1); - Wid = sqrt(Wid_gA(idx(:),:).^2 + Wid.^2); + lw = sqrt(lw_gA(idx(:),:).^2 + lw.^2); else - StrainMatrix = gStrainMatrix; for iOri = 1:nOrientations - lw2(1,iOri) = vecs(:,iOri).'*StrainMatrix.^2*vecs(:,iOri); + lw2(1,iOri) = vecs(:,iOri).'*R_g2M*lw_g.^2*R_g2M.'*vecs(:,iOri); end - Wid_gA = sqrt(lw2)*planck*1e6./geff/bmagn*1e3; % MHz -> mT - Wid = sqrt(repmat(Wid_gA.^2,2*S*nNucTrans,1)+Wid.^2); + lw_gA = sqrt(lw2)*planck*1e6./geff/bmagn*1e3; % MHz -> mT + lw = sqrt(repmat(lw_gA.^2,2*S*nNucTrans,1)+lw.^2); end elseif any(Sys.DStrain(:)) @@ -574,7 +578,7 @@ lwD = bsxfun(@times,lwD,MHz2mT); lwE = bsxfun(@times,lwE,MHz2mT); Wid2_DE = repmat(lwD.^2+lwE.^2,nNucTrans,1); - Wid = sqrt(Wid2_DE + Wid.^2); + lw = sqrt(Wid2_DE + lw.^2); end % Transitions @@ -599,12 +603,12 @@ siz = [nTransitions*nSites, numel(B)/nTransitions/nSites]; B = reshape(B,siz); if ~isempty(Int), Int = reshape(Int,siz); end - if ~isempty(Wid), Wid = reshape(Wid,siz); end + if ~isempty(lw), lw = reshape(lw,siz); end end % Arrange output %--------------------------------------------------------------- -Output = {B,Int,Wid,Transitions,spec}; +Output = {B,Int,lw,Transitions,spec}; varargout = Output(1:max(nargout,1)); end diff --git a/easyspin/resfreqs_matrix.m b/easyspin/resfreqs_matrix.m index b1a61e12..b54bf254 100644 --- a/easyspin/resfreqs_matrix.m +++ b/easyspin/resfreqs_matrix.m @@ -252,7 +252,7 @@ StrainsPresent = any([Sys.HStrain(:); Sys.DStrain(:); Sys.gStrain(:); Sys.AStrain(:)]); computeStrains = StrainsPresent && (nargout>2); -computeIntensities = ((nargout>1) & Opt.Intensity); +computeIntensities = nargout>1 && Opt.Intensity; % Preparing kernel and perturbing system Hamiltonians. @@ -534,7 +534,8 @@ % D strain %----------------------------------------------- - [useDStrain,dHdD,dHdE] = getdstrainops(CoreSys); + sigmadHdDE = strains_de(CoreSys); + useDStrain = ~isempty(sigmadHdDE); % g-A strain %------------------------------------------------- @@ -547,7 +548,8 @@ gStrainMatrix{iEl} = diag(CoreSys.gStrain(iEl,:)./CoreSys.g(iEl,:)); if any(CoreSys.gFrame(iEl,:)) R_g2M = erot(CoreSys.gFrame(iEl,:)).'; % g frame -> molecular frame - gStrainMatrix{iEl} = R_g2M*gStrainMatrix{iEl}*R_g2M.'; + else + R_g2M = eye(3); end end if ~simplegStrain @@ -560,7 +562,7 @@ end else for iEl = CoreSys.nElectrons:-1:1 - gStrainMatrix{iEl} = 0; + gStrainMatrix{iEl} = zeros(3); end end @@ -644,9 +646,9 @@ % Set up Hamiltonians for 3 lab principal directions %----------------------------------------------------- - [xLab,yLab,zLab] = erot(Orientations(iOri,:),'rows'); + [xLab_M,yLab_M,zLab_M] = erot(Orientations(iOri,:),'rows'); if higherOrder - [Vs,E] = gethamdata_hO(Exp.Field,zLab, CoreSys,Opt.Sparse, [], nLevels); + [Vs,E] = gethamdata_hO(Exp.Field,zLab_M, CoreSys,Opt.Sparse, [], nLevels); if Opt.Sparse, sp = 'sparse'; else, sp = ''; end [g0e{1},g0e{2},g0e{3}] = ham_ez(CoreSys,[],sp); [g0n{1},g0n{2},g0n{3}] = ham_nz(CoreSys,[],sp); @@ -659,18 +661,18 @@ kGM{n} = g0{n}-g1{1}{n}; end % z laboratoy axis: external static field - kmuzL = zLab(1)*kGM{1} + zLab(2)*kGM{2} + zLab(3)*kGM{3}; + kmuzL = zLab_M(1)*kGM{1} + zLab_M(2)*kGM{2} + zLab_M(3)*kGM{3}; % x laboratory axis: B1 excitation field - kmuxL = xLab(1)*kGM{1} + xLab(2)*kGM{2} + xLab(3)*kGM{3}; + kmuxL = xLab_M(1)*kGM{1} + xLab_M(2)*kGM{2} + xLab_M(3)*kGM{3}; % y laboratory vector: needed for integration over all B1 field orientations. - kmuyL = yLab(1)*kGM{1} + yLab(2)*kGM{2} + yLab(3)*kGM{3}; + kmuyL = yLab_M(1)*kGM{1} + yLab_M(2)*kGM{2} + yLab_M(3)*kGM{3}; else % z laboratoy axis: external static field - kmuzL = zLab(1)*kmuxM + zLab(2)*kmuyM + zLab(3)*kmuzM; + kmuzL = zLab_M(1)*kmuxM + zLab_M(2)*kmuyM + zLab_M(3)*kmuzM; % x laboratory axis: B1 excitation field - kmuxL = xLab(1)*kmuxM + xLab(2)*kmuyM + xLab(3)*kmuzM; + kmuxL = xLab_M(1)*kmuxM + xLab_M(2)*kmuyM + xLab_M(3)*kmuzM; % y laboratory vector: needed for integration over all B1 field orientations. - kmuyL = yLab(1)*kmuxM + yLab(2)*kmuyM + yLab(3)*kmuzM; + kmuyL = yLab_M(1)*kmuxM + yLab_M(2)*kmuyM + yLab_M(3)*kmuzM; [Vs,E] = gethamdata(Exp.Field, kH0, kmuzL, [], nLevels); end @@ -777,32 +779,34 @@ % Calculate width if requested. %-------------------------------------------------- if computeStrains - LineWidthSquared = CoreSys.HStrain.^2*zLab.^2; + LineWidthSquared = CoreSys.HStrain.^2*zLab_M.^2; for iTrans = 1:nTransitions - m = @(Op)Vs(:,v(iTrans))'*Op*Vs(:,v(iTrans)) - Vs(:,u(iTrans))'*Op*Vs(:,u(iTrans)); - + V = Vs(:,v(iTrans)); + U = Vs(:,u(iTrans)); + deltaVU = @(Op) real(V'*Op*V - U'*Op*U); + % H strain: Frequency-domain residual width tensor LineWidth2 = LineWidthSquared; - % D strain + % D/E strain if useDStrain - for iEl = 1:CoreSys.nElectrons - LineWidth2 = LineWidth2 + abs(m(dHdD{iEl}))^2; - LineWidth2 = LineWidth2 + abs(m(dHdE{iEl}))^2; + for i = 1:numel(sigmadHdDE) + LineWidth2 = LineWidth2 + deltaVU(sigmadHdDE{i})^2; end end % A strain if useAStrain - LineWidth2 = LineWidth2 + abs(m(dHdAx))^2; - LineWidth2 = LineWidth2 + abs(m(dHdAy))^2; - LineWidth2 = LineWidth2 + abs(m(dHdAz))^2; + LineWidth2 = LineWidth2 + abs(deltaVU(dHdAx))^2; + LineWidth2 = LineWidth2 + abs(deltaVU(dHdAy))^2; + LineWidth2 = LineWidth2 + abs(deltaVU(dHdAz))^2; end % g strain if usegStrain - dg2 = (m(kmuzL)*Exp.Field*zLab.'*gStrainMatrix{1}*zLab)^2; - LineWidth2 = LineWidth2 + abs(dg2); + dg2 = R_g2M*gStrainMatrix{1}.^2*R_g2M.'; + dg2 = abs(deltaVU(kmuzL))*Exp.Field*1e3*dg2; + LineWidth2 = LineWidth2 + zLab_M.'*dg2*zLab_M; end Wdat(iTrans,iOri) = sqrt(LineWidth2); end @@ -831,7 +835,7 @@ % Nuclear Zeeman and quadrupole (independent of S) if ~Opt.HybridOnlyHFI Hc = Hquad{iiNuc} + Exp.Field*... - (zLab(1)*Hzeem(iiNuc).x + zLab(2)*Hzeem(iiNuc).y + zLab(3)*Hzeem(iiNuc).z); + (zLab_M(1)*Hzeem(iiNuc).x + zLab_M(2)*Hzeem(iiNuc).y + zLab_M(3)*Hzeem(iiNuc).z); Hu = Hu + Hc; Hv = Hv + Hc; end diff --git a/tests/pepper_dstraincorrelated.m b/tests/pepper_dstraincorrelated.m index c2385013..e3988f58 100644 --- a/tests/pepper_dstraincorrelated.m +++ b/tests/pepper_dstraincorrelated.m @@ -12,33 +12,36 @@ Opt.GridSize = 15; -Sys.DStrainCorr = 0; [x,y0] = pepper(Sys,Exp,Opt); -Sys.DStrainCorr = +1; [x,yp] = pepper(Sys,Exp,Opt); -Sys.DStrainCorr = -1; [x,ym] = pepper(Sys,Exp,Opt); +Sys.DStrainCorr = 0; [B,spc0] = pepper(Sys,Exp,Opt); +Sys.DStrainCorr = +1; [B,spcp] = pepper(Sys,Exp,Opt); +Sys.DStrainCorr = -1; [B,spcm] = pepper(Sys,Exp,Opt); -if (opt.Display) +if opt.Display if isempty(olddata) - plot(x,y0,'b',x,yp,'r',x,ym,'g'); + plot(B,spc0,'b',B,spcp,'r',B,spcm,'g'); legend('corr = 0','corr = +1','corr = -1'); legend boxoff else - subplot(3,1,1); plot(x,olddata.y0,'b',x,y0,'r'); - subplot(3,1,2); plot(x,olddata.yp,'b',x,yp,'r'); - subplot(3,1,3); plot(x,olddata.ym,'b',x,ym,'r'); + subplot(3,1,1); plot(B,olddata.y0,'b',B,spc0,'r'); + title('corr = 0'); + subplot(3,1,2); plot(B,olddata.yp,'b',B,spcp,'r'); + title('corr = +1'); + subplot(3,1,3); plot(B,olddata.ym,'b',B,spcm,'r'); + title('corr = -1'); legend('old','new'); legend boxoff; end end -data.y0 = y0; -data.yp = yp; -data.ym = ym; +data.y0 = spc0; +data.yp = spcp; +data.ym = spcm; if ~isempty(olddata) - m = max([y0 yp ym]); - ok = areequal(y0,olddata.y0,m*0.01,'abs'); - ok = ok && areequal(yp,olddata.yp,m*0.01,'abs'); - ok = ok && areequal(ym,olddata.ym,m*0.01,'abs'); + m = max([spc0 spcp spcm]); + ok = areequal(spc0,olddata.y0,m*0.01,'abs'); + ok = ok && areequal(spcp,olddata.yp,m*0.01,'abs'); + ok = ok && areequal(spcm,olddata.ym,m*0.01,'abs'); else ok = []; end diff --git a/tests/resfields_gstrain.m b/tests/resfields_gstrain.m new file mode 100644 index 00000000..66682da3 --- /dev/null +++ b/tests/resfields_gstrain.m @@ -0,0 +1,16 @@ +function ok = test() + +% Test whether g-strain-based broadening gives correct line width + +Sys.g = 2; +Sys.gStrain = 0.01; + +Exp.mwFreq = 9.5; +Exp.Range = [300 400]; + +[B,I,W] = resfields(Sys,Exp); + +ok = areequal(W/B,Sys.gStrain/Sys.g,1e-8,'rel'); + +end + diff --git a/tests/resfields_gstrainframe.m b/tests/resfields_gstrainframe.m new file mode 100644 index 00000000..14bc7578 --- /dev/null +++ b/tests/resfields_gstrainframe.m @@ -0,0 +1,24 @@ +function ok = test + +Sys.S = 1/2; +Sys.g = [2.000 2.005 2.010]; +Sys.gStrain = [0.0001 0.0003 0.005]; + +Exp.mwFreq = 34; % GHz +Exp.Range = [1204 1217]; % mT +Exp.Harmonic = 0; + +gFrame = [0 0 0; -20 90 0; -20 50 80; -90 30 0]*pi/180; +nFrames = size(gFrame,1); + +for i = 1:nFrames + Sys.gFrame = gFrame(i,:); + Exp.SampleFrame = -fliplr(gFrame(i,:)); + [B(i),I(i),W(i)] = resfields(Sys,Exp); +end + +for k = 1:nFrames-1 + ok(k) = areequal(W(k),W(end),1e-8,'rel'); +end + +end diff --git a/tests/resfields_pt_gstrain.m b/tests/resfields_pt_gstrain.m new file mode 100644 index 00000000..d72fd3c0 --- /dev/null +++ b/tests/resfields_pt_gstrain.m @@ -0,0 +1,19 @@ +function ok = test() + +% Test whether g strain linewidths are correct along principal axes + +Sys.g = [2 2 2]; +Sys.gStrain = [0.01 0.02 0.04]; +Exp.mwFreq = 9.5; +Exp.Range = [300 380]; + +Exp.SampleFrame = [0 0 0]; +[Bz,~,Wz] = resfields_perturb(Sys,Exp); +Exp.SampleFrame = [0 pi/2 0]; +[Bx,~,Wx] = resfields_perturb(Sys,Exp); +Exp.SampleFrame = [0 pi/2 pi/2]; +[By,~,Wy] = resfields_perturb(Sys,Exp); + +ok = areequal([Wx/Bx Wy/By Wz/Bz],Sys.gStrain./Sys.g,1e-10,'rel'); + +end \ No newline at end of file diff --git a/tests/resfields_pt_gstrainframe.m b/tests/resfields_pt_gstrainframe.m new file mode 100644 index 00000000..ceeb800d --- /dev/null +++ b/tests/resfields_pt_gstrainframe.m @@ -0,0 +1,24 @@ +function ok = test() + +Sys.S = 1/2; +Sys.g = [2.000 2.005 2.010]; +Sys.gStrain = [0.0001 0.0003 0.005]; + +Exp.mwFreq = 34; % GHz +Exp.Range = [1204 1217]; % mT +Exp.Harmonic = 0; + +gFrame = [0 0 0; -20 90 0; -20 50 80; -90 30 0]*pi/180; +nFrames = size(gFrame,1); + +for i = 1:nFrames + Sys.gFrame = gFrame(i,:); + Exp.SampleFrame = -fliplr(gFrame(i,:)); + [~,~,W(i)] = resfields_perturb(Sys,Exp); +end + +for k = 1:nFrames-1 + ok(k) = areequal(W(k),W(end),1e-10,'rel'); +end + +end diff --git a/tests/resfreqs_matrix_gstrain.m b/tests/resfreqs_matrix_gstrain.m index 629b8871..3e975ce5 100644 --- a/tests/resfreqs_matrix_gstrain.m +++ b/tests/resfreqs_matrix_gstrain.m @@ -1,27 +1,34 @@ function [ok,data] = test(opt,olddata) -% Check whether resfreqs_matrix handlex gStrain correctly. - +% Check whether resfreqs_matrix handles gStrain correctly. Sys.S = 1/2; Sys.g = 2; Sys.gStrain = 0.01; -Exp.Field = 350; -Exp.mwRange = [9.4 10]; +mw = 9.8; % GHz +Exp.Field = mw*1e9*planck/bmagn/Sys.g/1e-3; % mT +Exp.mwRange = mw + [-1 1]*0.1; % GHz -[x,y] = pepper(Sys,Exp); -y = y/max(y); +[B,spc] = pepper(Sys,Exp); +spc = spc/max(spc); if opt.Display - plot(x,y,x,olddata.y); + plot(B,spc,B,olddata.y); + yline(0.5); + xline(mw); + g_up = Sys.g + Sys.gStrain/2; + g_lo = Sys.g - Sys.gStrain/2; + dmw_up = g_up*bmagn*Exp.Field*1e-3/planck/1e9; + dmw_lo = g_lo*bmagn*Exp.Field*1e-3/planck/1e9; + xline(dmw_up); + xline(dmw_lo); legend('new','old'); end -data.y = y; +data.y = spc; if ~isempty(olddata) - ok = areequal(y,olddata.y,1e-3,'abs'); + ok = areequal(spc,olddata.y,1e-3,'abs'); else ok = []; end - From 2d4b452b5816f22ece1655f166efa3a8f8db1bb6 Mon Sep 17 00:00:00 2001 From: fmentink Date: Tue, 1 Apr 2025 16:09:33 -0400 Subject: [PATCH 2/6] updated functions ee, ez, hf, zf with derivatives --- easyspin/ham_eeD.m | 146 +++++++++++++++++++++++++++++++++++++ easyspin/ham_ezD.m | 173 ++++++++++++++++++++++++++++++++++++++++++++ easyspin/ham_hfD.m | 139 +++++++++++++++++++++++++++++++++++ easyspin/ham_zfD.m | 176 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 634 insertions(+) create mode 100644 easyspin/ham_eeD.m create mode 100644 easyspin/ham_ezD.m create mode 100644 easyspin/ham_hfD.m create mode 100644 easyspin/ham_zfD.m diff --git a/easyspin/ham_eeD.m b/easyspin/ham_eeD.m new file mode 100644 index 00000000..37dc4faa --- /dev/null +++ b/easyspin/ham_eeD.m @@ -0,0 +1,146 @@ +% ham_ee Electron-electron spin interaction Hamiltonian +% +% F = ham_ee(SpinSystem) +% F = ham_ee(SpinSystem,eSpins) +% F = ham_ee(SpinSystem,eSpins,'sparse') +% +% Returns the electron-electron spin interaction (EEI) +% Hamiltonian, in MHz. +% +% Input: +% - SpinSystem: Spin system structure. EEI +% parameters are in the ee, eeFrame, and ee2 fields. +% - eSpins: If given, specifies electron spins +% for which the EEI should be computed. If +% absent, all electrons are included. +% - 'sparse': If given, the matrix is returned in sparse format. + +% Output: +% - F: Hamiltonian matrix containing the EEI for +% electron spins specified in eSpins. + +function [F,dF] = ham_ee(System,Spins,opt) + +if nargin==0, help(mfilename); return; end + +if nargin<1 || nargin>3, error('Wrong number of input arguments!'); end +if nargout<0, error('Not enough output arguments.'); end +if nargout>2, error('Too many output arguments.'); end + +if nargin<3, opt = ''; end +if nargin<2, Spins = []; end +if ~ischar(opt) + error('Third input must be a string, ''sparse''.'); +end +useSparseMatrices = strcmp(opt,'sparse'); + +[System,err] = validatespinsys(System); +error(err); + +sys = spinvec(System); +n = prod(2*sys+1); + +% Special cases: only one spins, ee not given or all zero +F = sparse(n,n); +if (System.nElectrons==1), return; end +if ~any(System.ee(:)) && ~any(System.ee2(:)), return; end + +if isempty(Spins), Spins = 1:System.nElectrons; end + +% Some error checking on the second input argument +if numel(Spins)<2 + error('Spins (2nd argument) must contain at least 2 values!'); +end +if any(Spins<1) || any(Spins>System.nElectrons) + error('Spins (2nd argument) must contain values between 1 and %d!',System.nElectrons); +end +if numel(unique(Spins))~=numel(Spins) + error('Spins (2nd argument) contains double entries!'); +end + +F = sparse(n,n); + +% Compile list of wanted interactions +Spins = sort(Spins); +[idx1,idx2] = find(tril(ones(numel(Spins)),-1)); +idx = [idx2,idx1]; + +Pairs = Spins(idx); +nPairs = size(Pairs,1); +Coupl = Pairs(:,1) + (Pairs(:,2)-1)*System.nElectrons; + +% Compile list of all spin pairs +[e2,e1] = find(tril(ones(System.nElectrons),-1)); +allPairsIdx = e1 + (e2-1)*System.nElectrons; + +ee = System.ee; +if isfield(System,'eeFrame') + eeFrame = System.eeFrame; +else + eeFrame = []; +end + +%for conistency accross methods +nStates=System.nStates; +% Bilinear coupling term S1*ee*S2 +%---------------------------------------------------------------- +for iPair = 1:nPairs + iCoupling = find(Coupl(iPair)==allPairsIdx); + + % Construct matrix representing coupling tensor + if System.fullee + J = ee(3*(iCoupling-1)+(1:3),:); + else + J = diag(ee(iCoupling,:)); + end + if ~isempty(eeFrame) && any(eeFrame(iCoupling,:)) + R_M2ee = erot(eeFrame(iCoupling,:)); % mol frame -> ee frame + R_ee2M = R_M2ee.'; % ee frame -> mol frame + J = R_ee2M*J*R_ee2M.'; + else + R_ee2M = eye(3); + end + + % preparing the derivatives (specific for each electron) + dFdeex = sparse(nStates,nStates); + dFdeey = sparse(nStates,nStates); + dFdeez = sparse(nStates,nStates); + + % preparing the derivatives coefficient + deexM = R_ee2M(:,1)*R_ee2M(:,1).'; % rotate derivative wrt eex to molecular frame + deeyM = R_ee2M(:,2)*R_ee2M(:,2).'; % rotate derivative wrt eey to molecular frame + deezM = R_ee2M(:,3)*R_ee2M(:,3).'; % rotate derivative wrt eez to molecular frame + + % Sum up Hamiltonian terms + for c1 = 1:3 + so1 = sop(sys,[Pairs(iPair,1),c1],'sparse'); + for c2 = 1:3 + so2 = sop(sys,[Pairs(iPair,2),c2],'sparse'); + tempProduct=so1*so2; + F = F + J(c1,c2)*tempProduct; + dFdeex = dFdeex + deexM(c1,c2)*tempProduct; + dFdeey = dFdeey + deeyM(c1,c2)*tempProduct; + dFdeez = dFdeez + deezM(c1,c2)*tempProduct; + end + end + + dF{iPair} = {dFdeex,dFdeey,dFdeez}; % derivatives + + % Isotropic biquadratic exchange coupling term +ee2*(S1.S2)^2 + %----------------------------------------------------------------- + if System.ee2(iCoupling)==0, continue; end + F2 = 0; + for c = 1:3 + F2 = F2 + sop(sys,[Pairs(iPair,1),c;Pairs(iPair,2),c],'sparse'); + end + dFdJ=F2^2; + F = F + System.ee2(iCoupling)*dFdJ; + dF{iPair} = {dFdeex,dFdeey,dFdeez,dFdJ}; % derivatives +end + +F = (F+F')/2; % Hermitianise +if ~useSparseMatrices + F = full(F); % sparse -> full +end + +end diff --git a/easyspin/ham_ezD.m b/easyspin/ham_ezD.m new file mode 100644 index 00000000..4764deba --- /dev/null +++ b/easyspin/ham_ezD.m @@ -0,0 +1,173 @@ +% ham_ez Electron Zeeman interaction Hamiltonian +% +% H = ham_ez(SpinSystem, B) +% H = ham_ez(SpinSystem, B, eSpins) +% H = ham_ez(SpinSystem, B, eSpins, 'sparse') +% [mux, muy, muz] = ham_ez(SpinSystem) +% [mux, muy, muz] = ham_ez(SpinSystem, eSpins) +% [mux, muy, muz] = ham_ez(SpinSystem, eSpins, 'sparse') +% +% Returns the electron Zeeman interaction Hamiltonian matrix for +% the electron spins 'eSpins' of the spin system 'SpinSystem'. +% +% Input: +% - SpinSystem: Spin system structure. +% - B: Magnetic field vector, in millitesla. +% - eSpins: Vector of indices for electron spins to include. If eSpins is +% omitted, all electron spins are included. +% - 'sparse': If given, results returned in sparse format. +% +% Output: +% - mux, muy, muz: Components of the magnetic dipole moment operator +% for the selected electron spins as defined by mui=-d(H)/d(B_i) +% i=x,y,z where B_i are the cartesian components of +% the external field in the molecular frame. Units are MHz/mT = GHz/T. +% To get the full electron Zeeman Hamiltonian, use +% H = -(mux*B(1)+muy*B(2)+muz*B(3)), where B is the magnetic field vector +% in mT. +% - H: Electron Zeeman Hamiltonian matrix. + +function varargout = ham_ez(SpinSystem,varargin) + +if nargin==0, help(mfilename); return; end + +if nargout~=1 && nargout~=6, error('Wrong number of output arguments!'); end +singleOutput = nargout==1; + +if singleOutput + if nargin<1 || nargin>4 + error('Wrong number of input arguments!'); + end + if nargin<2 + error('Field vector (second input, in mT) is missing.') + else + B0 = varargin{1}; + end + if nargin<3, eSpins = []; else, eSpins = varargin{2}; end + if nargin<4, opt = ''; else, opt = varargin{3}; end +else + B0 = []; + if nargin<2, eSpins = []; else, eSpins = varargin{1}; end + if nargin<3, opt = ''; else, opt = varargin{2}; end +end + +if ~ischar(opt) + error('Last input must be a string, ''sparse''.'); +end +useSparseMatrices = strcmp(opt,'sparse'); + +% Validate spin system +[Sys,err] = validatespinsys(SpinSystem); +error(err); + +% Vector of all spin quantum numbers +spins = Sys.Spins; +nElectrons = Sys.nElectrons; + +% No 'eSpins' specified -> use all +if isempty(eSpins), eSpins = 1:nElectrons; end + +% Validate magnetic field if given +if ~isempty(B0) + if numel(B0)~=3 + error('Magnetic field vector (2nd input) must be a 3-element array.'); + end +end + +% Validate third argument (eSpins) +if any(eSpins<1) || any(eSpins>nElectrons) + error('Electron spin indices (2nd input argument) invalid!'); +end + +% Initialize Zeeman interaction component matrices to zero +nStates = Sys.nStates; +muxM = sparse(nStates,nStates); +muyM = sparse(nStates,nStates); +muzM = sparse(nStates,nStates); + +% Calculate prefactors +pre = -bmagn/planck*Sys.g; % Hz/T +pre = pre/1e9; % Hz/T -> GHz/T = MHz/mT + +% Loop over all electron spins +for i = eSpins + + % Get full g matrix + if Sys.fullg + g = pre((i-1)*3+(1:3),:); + else + g = diag(pre(i,:)); + end + + % Transform g matrix to molecular frame + ang = Sys.gFrame(i,:); + if any(ang) + R_M2g = erot(ang); % mol frame -> g frame + R_g2M = R_M2g.'; % g frame -> mol frame + g = R_g2M*g*R_g2M.'; + else + R_g2M=eye(3); + end + + % preparing the derivatives (specific for each electron) + dMuxxM = sparse(nStates,nStates); + dMuxyM = sparse(nStates,nStates); + dMuxzM = sparse(nStates,nStates); + + dMuyxM = sparse(nStates,nStates); + dMuyyM = sparse(nStates,nStates); + dMuyzM = sparse(nStates,nStates); + + dMuzxM = sparse(nStates,nStates); + dMuzyM = sparse(nStates,nStates); + dMuzzM = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dgxM = R_g2M(:,1)*R_g2M(:,1).'; % rotate derivative wrt gx to molecular frame + dgyM = R_g2M(:,2)*R_g2M(:,2).'; % rotate derivative wrt gy to molecular frame + dgzM = R_g2M(:,3)*R_g2M(:,3).'; % rotate derivative wrt gz to molecular frame + + % Build magnetic dipole moment components in MHz/mT + for k = 1:3 + Sk = sop(spins,[i,k],'sparse'); + muxM = muxM + g(1,k)*Sk; + muyM = muyM + g(2,k)*Sk; + muzM = muzM + g(3,k)*Sk; + + dMuxxM = dMuxxM + dgxM(1,k)*Sk; + dMuyxM = dMuyxM + dgxM(2,k)*Sk; + dMuzxM = dMuzxM + dgxM(3,k)*Sk; + + dMuxyM = dMuxyM + dgyM(1,k)*Sk; + dMuyyM = dMuyyM + dgyM(2,k)*Sk; + dMuzyM = dMuzyM + dgyM(3,k)*Sk; + + dMuxzM = dMuxzM + dgzM(1,k)*Sk; + dMuyzM = dMuyzM + dgzM(2,k)*Sk; + dMuzzM = dMuzzM + dgzM(3,k)*Sk; + end + +end + +if isempty(B0) + % Return magnetic dipole moment components + if ~useSparseMatrices + muxM = full(muxM); + muyM = full(muyM); + muzM = full(muzM); + + dMuxM = {dMuxxM, dMuyxM, dMuzxM}; + dMuyM = {dMuxyM, dMuyyM, dMuzyM}; + dMuzM = {dMuxzM, dMuyzM, dMuzzM}; + end + varargout = {muxM, muyM, muzM, dMuxM, dMuyM, dMuzM}; +else + % Return Zeeman Hamiltonian + H = -(muxM*B0(1) + muyM*B0(2) + muzM*B0(3)); + if ~useSparseMatrices + H = full(H); + end + varargout = {H}; +end + +end diff --git a/easyspin/ham_hfD.m b/easyspin/ham_hfD.m new file mode 100644 index 00000000..45917045 --- /dev/null +++ b/easyspin/ham_hfD.m @@ -0,0 +1,139 @@ +% ham_hf Hyperfine interaction Hamiltonian +% +% Hhf = ham_hf(System) +% Hhf = ham_hf(System,eSpins) +% Hhf = ham_hf(System,eSpins,nSpins) +% Hhf = ham_hf(System,eSpins,nSpins,'sparse') +% +% Returns the hyperfine interaction Hamiltonian (in units of MHz) between +% electron spins 'eSpins' and nuclear spins 'nSpins' of the system +% 'System'. eSpins=1 is the first electron spins, nSpins=1 is the first +% nuclear spin. +% +% If 'sparse' is given, the matrix is returned in sparse format. + +function [Hhf,dHhf] = ham_hf(System,elSpins,nucSpins,opt) + +if nargin==0 + help(mfilename); + return +end + +switch nargin + case 1 + elSpins = []; + nucSpins = []; + opt = ''; + case 2 + nucSpins = []; + opt = ''; + case 3 + opt = ''; + case 4 + otherwise + error('Incorrect number of input arguments.') +end + +if ~ischar(opt) + error('Fourth input must be a string, ''sparse''.'); +end +useSparseMatrices = strcmp(opt,'sparse'); + +% Validate spin system. +[Sys,err] = validatespinsys(System); +error(err); + +% Get spin vector and space dimension +SpinVec = Sys.Spins; +nStates = Sys.nStates; +nElectrons = Sys.nElectrons; +nNuclei = Sys.nNuclei; + +Hhf = sparse(nStates,nStates); % sparse zero matrix + +% Special case: no nuclei present, so no hyperfine +if nNuclei==0 + if ~useSparseMatrices + Hhf = full(Hhf); + end + return +end + +% Get electron spin list +if isempty(elSpins) + elSpins = 1:nElectrons; +else + if any(elSpins<1) || any(elSpins>nElectrons) + error('Electron spins (2nd argument) contains out-of-range values!'); + end +end + +% Get electron spin list +if isempty(nucSpins) + nucSpins = 1:nNuclei; +else + if any(nucSpins<1) || any(nucSpins>nNuclei) + error('Nuclear spins (3rd argument) contains out-of-range values!'); + end +end + +fullAMatrix = size(System.A,1)>nNuclei; + +% Generate Hamiltonian for hyperfine interaction +for eSp = elSpins + eidx = (eSp-1)*3+(1:3); + for nSp = nucSpins + if Sys.I(nSp)==0, continue; end + + % Construct full hyperfine matrix + if fullAMatrix + A = Sys.A((nSp-1)*3+(1:3),eidx); + else + A = diag(Sys.A(nSp,eidx)); + end + + if ~any(A(:)) + continue + end + + % Transform matrix into molecular frame representation + ang = Sys.AFrame(nSp,eidx); + if any(ang) + R_M2A = erot(ang); % mol frame -> A frame + R_A2M = R_M2A.'; % A frame -> mol frame + A = R_A2M*A*R_A2M.'; + else + R_A2M=eye(3); + end + + % preparing the derivatives (specific for each electron) + dHhfx = sparse(nStates,nStates); + dHhfy = sparse(nStates,nStates); + dHhfz = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dAxM = R_A2M(:,1)*R_A2M(:,1).'; % rotate derivative wrt Ax to molecular frame + dAyM = R_A2M(:,2)*R_A2M(:,2).'; % rotate derivative wrt Ay to molecular frame + dAzM = R_A2M(:,3)*R_A2M(:,3).'; % rotate derivative wrt Az to molecular frame + + % Construct hyperfine Hamiltonian + for c1 = 1:3 + for c2 = 1:3 + tempProduct=sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + Hhf = Hhf + A(c1,c2)*tempProduct; + dHhfx = dHhfx + dAxM(c1,c2)*tempProduct; + dHhfy = dHhfy + dAyM(c1,c2)*tempProduct; + dHhfz = dHhfz + dAzM(c1,c2)*tempProduct; + end + end + + end % for all specified nuclei + dHhf{eSp,nSp}={dHhfx,dHhfy,dHhfz}; % --> the indices in the derivative may not be needed? +end % for all specified electrons + +Hhf = (Hhf+Hhf')/2; % hermitianise, e.g. guards against small imaginary remainders on the diagonal +if ~useSparseMatrices + Hhf = full(Hhf); % convert sparse to full +end + +end diff --git a/easyspin/ham_zfD.m b/easyspin/ham_zfD.m new file mode 100644 index 00000000..91f3e5af --- /dev/null +++ b/easyspin/ham_zfD.m @@ -0,0 +1,176 @@ +% ham_zf Electronic zero field interaction Hamiltonian +% +% F = ham_zf(SpinSystem) +% F = ham_zf(SpinSystem,Electrons) +% F = ham_zf(SpinSystem,Electrons,'sparse') +% +% Returns the electronic zero-field interaction (ZFI) +% Hamiltonian of the system SpinSystem, in units of MHz. +% +% Can retrun the derivative of the ZFI if Sys.Dstrain has been given +% +% If the vector Electrons is given, the ZFI of only the +% specified electrons is returned (1 is the first, 2 the +% second, etc). Otherwise, all electrons are included. +% +% If 'sparse' is given, the matrix is returned in sparse format. + +function [H,dH] = ham_zf(SpinSystem,idxElectrons,opt) + +if nargin==0, help(mfilename); return; end + +[Sys,err] = validatespinsys(SpinSystem); +error(err); + +if nargin<2, idxElectrons = []; end +if nargin<3, opt = ''; end +if ~ischar(opt) + error('Third input must be a string, ''sparse''.'); +end +sparseResult = strcmp(opt,'sparse'); + +if isempty(idxElectrons) + idxElectrons = 1:Sys.nElectrons; +end + +if any(idxElectrons>Sys.nElectrons) || any(idxElectrons<1) + error('Electron spin index/indices (2nd argument) out of range!'); +end + +nStates=Sys.nStates; +H = sparse(nStates,nStates); + +Spins = Sys.Spins; + +% Quadratic term S*D*S +%------------------------------------------------------------------------------- +for iSpin = idxElectrons + + % Prepare full 3x3 D matrix + if Sys.fullD + D = Sys.D(3*(iSpin-1)+(1:3),:); + else + D = diag(Sys.D(iSpin,:)); + end + + if ~any(D(:)) + continue + end + + % Transform D tensor to molecular frame and obtain derivatives + ang = Sys.DFrame(iSpin,:); + if any(ang) + R_M2D = erot(ang); % mol frame -> D frame + R_D2M = R_M2D.'; % D frame -> mol frame + D = R_D2M*D*R_D2M.'; + else + R_D2M = eye(3); + end + + % preparing the derivatives (specific for each electron) + dHdDx = sparse(nStates,nStates); + dHdDy = sparse(nStates,nStates); + dHdDz = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dDxM = R_D2M(:,1)*R_D2M(:,1).'; % rotate derivative wrt Dx to molecular frame + dDyM = R_D2M(:,2)*R_D2M(:,2).'; % rotate derivative wrt Dy to molecular frame + dDzM = R_D2M(:,3)*R_D2M(:,3).'; % rotate derivative wrt Dz to molecular frame + + % Construct spin operator matrices + for c = 3:-1:1 + Sxyz{c} = sop(Spins,[iSpin,c],'sparse'); + end + + % Construct SDS term + for c1 = 1:3 + for c2 = 1:3 + tempProduct=Sxyz{c1}*Sxyz{c2}; % storing the matrix product temporarily + H = H + D(c1,c2)*tempProduct; + dHdDx = dHdDx + dDxM(c1,c2)*tempProduct; + dHdDy = dHdDy + dDyM(c1,c2)*tempProduct; + dHdDz = dHdDz + dDzM(c1,c2)*tempProduct; + end + end + + dH{iSpin} = {dHdDx,dHdDy,dHdDz}; +end + +% Fourth-order terms a and F +%----------------------------------------------------------------------------- +% Spin system fields: a and F + +% Abragam/Bleaney, pp. 142, 437 +% Bleaney/Trenam, Proc.Roy.Soc.A, 223(1152), 1-14, (1954) +% Doetschman/McCool, Chem.Phys. 8, 1-16 (1975) +% Scullane, J.Magn.Reson. 47, 383-397 (1982) +% Jain/Lehmann, phys.stat.sol.(b) 159, 495-544 (1990) + +% These terms are no longer supported. +if isfield(Sys,'aF') + error('Sys.aF is no longer supported. Use Sys.B4 and Sys.B4Frame instead.'); +end + +% High-order terms in extended Stevens operator format +%----------------------------------------------------------------------------- +% Spin system fields: +% B1, B2, B3, ... (k = 1...12) -> processed to Sys.B +% B1Frame, B2Frame, B3Frame, ... -> processed to Sys.BFrame + +% Run over all ranks k +for k = 1:numel(Sys.B) + Bk = Sys.B{k}; + if isempty(Bk), continue; end + + % Run over all desired electron spins + for iSpin = idxElectrons + + % If D is used, skip rank-2 Stevens operator terms + D_present = isfield(Sys,'D') && ~isempty(D) && any(D(:)~=0); + if D_present && k==2, continue; end + + % Skip if all rank-k coefficients are zero + if all(Bk(iSpin,:)==0), continue; end + + % Apply transformation if non-zero tilt angles are given + % (Sys.BFrame is processed from Sys.B?Frame by validatespinsys) + tiltang = Sys.BFrame{k}(iSpin,:); + if any(tiltang) + % Calculate transformation matrix for ISTOs + Dk = wignerd(k,tiltang(1),tiltang(2),tiltang(3)); + Ck = isto2stev(k); + % Transform to Stevens operator basis + DB = Ck*Dk.'/Ck; + DB = removenumericalnoise(DB); + % Transform Bk + Bk_M = Bk(iSpin,:)*DB; % eigenframe of Bk -> molecular frame + else + Bk_M = Bk(iSpin,:); + end + + % Build Hamiltonian + q = k:-1:-k; + for iq = find(Bk_M~=0) + H = H + Bk_M(iq)*stev(Spins,[k,q(iq),iSpin],'sparse'); + end + + end % for all electron spins specified + +end % for all tensor ranks + +H = (H+H')/2; % Hermitianize +if ~sparseResult + H = full(H); +end + +end +%=============================================================================== + +function A_ = removenumericalnoise(A) +reA = real(A); +imA = imag(A); +thr = 1e-14; +reA(abs(reA) Date: Tue, 1 Apr 2025 16:15:19 -0400 Subject: [PATCH 3/6] added eeD, ezD hfD and zfsD for derivatives --- easyspin/private/strains_de.m | 4 +- easyspin/private/strains_g.m | 175 ++++++++++++++++++++++++++++++++++ easyspin/private/strains_ga.m | 30 +++--- 3 files changed, 193 insertions(+), 16 deletions(-) create mode 100644 easyspin/private/strains_g.m diff --git a/easyspin/private/strains_de.m b/easyspin/private/strains_de.m index ef4bebd6..8e4aa357 100644 --- a/easyspin/private/strains_de.m +++ b/easyspin/private/strains_de.m @@ -32,7 +32,7 @@ Delta = Sys.DStrain(iEl,:); % FWHMs for D and E rDE = Sys.DStrainCorr(iEl); % correlation coefficient between D and E if rDE<-1 || rDE>1 - error('Electron spin %d: D/E correlation coefficient must be between -1 and 1',iEl); + error('Electron spin %d: D/E correlation coefficient must be between -1 and 1, verify your input',iEl); end % Construct spin operators in molecular frame @@ -42,7 +42,7 @@ % Transform spin operators to D frame if any(Sys.DFrame(iEl,:)) - R_M2D = erot(Sys.DFrame(iEl,:)); % molecular frame -> D frame + R_M2D = erot(Sys.DFrame(iEl,:)); % molecular frame -> D frame -----------------------M2D pr D2M? confused due to g2M in strains_ga SxD = R_M2D(1,1)*SxM + R_M2D(1,2)*SyM + R_M2D(1,3)*SzM; SyD = R_M2D(2,1)*SxM + R_M2D(2,2)*SyM + R_M2D(2,3)*SzM; SzD = R_M2D(3,1)*SxM + R_M2D(3,2)*SyM + R_M2D(3,3)*SzM; diff --git a/easyspin/private/strains_g.m b/easyspin/private/strains_g.m new file mode 100644 index 00000000..831efd22 --- /dev/null +++ b/easyspin/private/strains_g.m @@ -0,0 +1,175 @@ +% Calculate Hamiltonian derivatives with respect to the g Principle Axis values, taking into +% account correlation if present; premultipy with linewidth. + +% Required input fields +% Sys.gStrain: +% nElectrons x 3 array with [gxx gyy gzz] on each row. +% Sys.gStrainCorr: +% Array with nElectron correlation coefficient between gxx, gyy and gzz (between -1 and +1). +% Sys.gFrame: +% Euler angles describing g tensor orientation in molecular frame + +function dHdDg = strains_g(Sys) + +% Return if no g strain present +if ~any(Sys.gStrain(:)) + dHdDg = {}; + return +end + +% Loop over all electron spins +nElectrons = size(Sys.gStrain,1); +dHdDg = cell(3,nElectrons); % preallocate maximal possible size 3*number of electrons +idx = 0; +for iEl = 1:nElectrons + + % Skip if no g strain is given for this electron spin + if ~any(Sys.gStrain(iEl,:)) + continue + end + + % Strain information + Delta = Sys.gStrain(iEl,:); % FWHMs for g + rgg = Sys.gStrainCorr(:,:,iEl); % correlation coefficient between g components (3x3 matrix per electron) --> should be a 3D matrix? + if find(rgg)<-1 || find(rgg)>1 + error('Electron spin %d: gi/gj correlation coefficient must be between -1 and 1, verify your input',iEl); + end + + % Construct spin operators in molecular frame + SxM = sop(Sys,[iEl,1]); + SyM = sop(Sys,[iEl,2]); + SzM = sop(Sys,[iEl,3]); + + % Transform spin operators to g frame + if any(Sys.gFrame(iEl,:)) + R_g2M = erot(Sys.gFrame(iEl,:)).'; % molecular frame -> g Frame + + + SxD = R_M2D(1,1)*SxM + R_M2D(1,2)*SyM + R_M2D(1,3)*SzM; + SyD = R_M2D(2,1)*SxM + R_M2D(2,2)*SyM + R_M2D(2,3)*SzM; + SzD = R_M2D(3,1)*SxM + R_M2D(3,2)*SyM + R_M2D(3,3)*SzM; + else + % g frame aligns with molecular frame + SxD = SxM; + SyD = SyM; + SzD = SzM; + end + + % Calculate derivatives of Hamiltonian w.r.t. gxx, gyy and gzz + + % H = D*(SzD^2-S(S+1)/3) + E*(SxD^2-SyD^2) + % = D*(2*SzD^2-SxD^2-SyD^2)/3 + E*(SxD^2-SyD^2) + dHdD = (2*SzD^2 - SxD^2 - SyD^2)/3; + dHdE = SxD^2 - SyD^2; + + % Diagonalize covariance matrix in the case of non-zero correlation + if rDE~=0 + % Transform correlated D-E strain to uncorrelated coordinates + logmsg(1,' correlated D strain for electron spin %d (r = %f)',iEl,rDE); + % Construct covariance matrix + R12 = rDE*Delta(1)*Delta(2); % covariance + covMatrix = [Delta(1)^2 R12; R12 Delta(2)^2]; + % Diagonalize covariance matrix, get derivatives along principal directions + [V,sigma2] = eig(covMatrix); + Delta = sqrt(diag(sigma2)); + dHdDE1 = Delta(1)*(V(1,1)*dHdD + V(1,2)*dHdE); + dHdDE2 = Delta(2)*(V(2,1)*dHdD + V(2,2)*dHdE); + else + dHdDE1 = Delta(1)*dHdD; + dHdDE2 = Delta(2)*dHdE; + end + + % Store in list + idx = idx+1; + dHdDE{idx} = dHdDE1; + idx = idx+1; + dHdDE{idx} = dHdDE2; + +end + +dHdDE = dHdDE(1:idx); + +end +% Calculates quantities needed by resfields for g/A strain broadening + +function [usegStrain,useAStrain,simplegStrain,lw2_gA,ops] = strains_ga(CoreSys,Exp,mwFreq,Transitions,nTransitions,kH0,kmuzM) + +usegStrain = any(CoreSys.gStrain(:)); +useAStrain = CoreSys.nNuclei>0 && any(CoreSys.AStrain); + +ops = struct; +lw2_gA = []; + +% g-A strain +%------------------------------------------------- +% g strain tensor is taken to be aligned with the g tensor +% A strain tensor is taken to be aligned with the A tensor +% g strain can be specified for each electron spin +% A strain is limited to the first electron and first nuclear spin +simplegStrain = CoreSys.nElectrons==1; +if usegStrain + logmsg(1,' g strain present'); + lw_g = mwFreq*CoreSys.gStrain./CoreSys.g; % MHz + for Ne = CoreSys.nElectrons:-1:1 + if any(CoreSys.gFrame(Ne,:)) + R_g2M{Ne} = erot(CoreSys.gFrame(Ne,:)).'; % g frame -> molecular frame [question about why .' not consistent with strains_de] + else + R_g2M{Ne} = eye(3); + end + end + if ~simplegStrain + logmsg(1,' multiple g strains present'); + for Ne = CoreSys.nElectrons:-1:1 + kSxM{Ne} = sop(CoreSys,[Ne,1]); + kSyM{Ne} = sop(CoreSys,[Ne,2]); + kSzM{Ne} = sop(CoreSys,[Ne,3]); + end + ops.kSxM = kSxM; + ops.kSyM = kSyM; + ops.kSzM = kSzM; + end +else + lw_g = zeros(CoreSys.nElectrons,3); +end + +if useAStrain + if usegStrain + if any(CoreSys.gFrame(1,:)~=CoreSys.AFrame(1,1:3)) + error('For g/A strain, Sys.gFrame and Sys.AFrame must be collinear.'); + end + R_strain2M = R_g2M; + else + R_A2M = erot(CoreSys.AFrame).'; % strain tensor frame -> mol. frame + R_strain2M = R_A2M; + end + lw_A = diag(CoreSys.AStrain); + % Diagonalize Hamiltonian at center field. + centerB = mean(Exp.Range); + [Vecs,E] = eig(kH0 - centerB*kmuzM); + [~,idx] = sort(real(diag(E))); + Vecs = Vecs(:,idx); + % Calculate effective mI of nucleus 1 for all eigenstates. + mI1 = real(diag(Vecs'*sop(CoreSys,[2,3])*Vecs)); + mI1Tr = mean(mI1(Transitions),2); + % compute A strain array + lw_A = reshape(mI1Tr(:,ones(1,9)).',[3,3,nTransitions]).*... + repmat(lw_A,[1,1,nTransitions]); + + % Combine g and A strain + rho = CoreSys.gAStrainCorr; % correlation coefficient + for Ne = CoreSys.nElectrons:-1:1 + lw2_gA{Ne} = repmat(diag(lw_g(Ne,:).^2),[1,1,nTransitions]) + lw_A.^2 + ... + 2*rho*repmat(diag(lw_g(Ne,:)),[1,1,nTransitions]).*lw_A; + for tr = 1:nTransitions + lw2_gA{Ne}(:,:,tr) = R_strain2M{Ne}*lw2_gA{Ne}(:,:,tr)*R_strain2M{Ne}.'; + end + end +else + if usegStrain + for Ne = CoreSys.nElectrons:-1:1 + lw2_gA{Ne} = repmat(R_g2M{Ne}*diag(lw_g(Ne,:).^2)*R_g2M{Ne}.',[1,1,nTransitions]); + end + end +end +% lw2_gA = a (cell array of) 3D array with 3x3 strain line-width matrices +% for each transition stacked along the third dimension. diff --git a/easyspin/private/strains_ga.m b/easyspin/private/strains_ga.m index ccde0ba5..7765ad19 100644 --- a/easyspin/private/strains_ga.m +++ b/easyspin/private/strains_ga.m @@ -18,19 +18,19 @@ if usegStrain logmsg(1,' g strain present'); lw_g = mwFreq*CoreSys.gStrain./CoreSys.g; % MHz - for e = CoreSys.nElectrons:-1:1 - if any(CoreSys.gFrame(e,:)) - R_g2M{e} = erot(CoreSys.gFrame(e,:)).'; % g frame -> molecular frame + for Ne = CoreSys.nElectrons:-1:1 + if any(CoreSys.gFrame(Ne,:)) + R_g2M{Ne} = erot(CoreSys.gFrame(Ne,:)).'; % g frame -> molecular frame [question about why .' not consistent with strains_de] else - R_g2M{e} = eye(3); + R_g2M{Ne} = eye(3); end end if ~simplegStrain logmsg(1,' multiple g strains present'); - for e = CoreSys.nElectrons:-1:1 - kSxM{e} = sop(CoreSys,[e,1]); - kSyM{e} = sop(CoreSys,[e,2]); - kSzM{e} = sop(CoreSys,[e,3]); + for Ne = CoreSys.nElectrons:-1:1 + kSxM{Ne} = sop(CoreSys,[Ne,1]); + kSyM{Ne} = sop(CoreSys,[Ne,2]); + kSzM{Ne} = sop(CoreSys,[Ne,3]); end ops.kSxM = kSxM; ops.kSyM = kSyM; @@ -41,6 +41,7 @@ end if useAStrain + %%I ignore it for now if usegStrain if any(CoreSys.gFrame(1,:)~=CoreSys.AFrame(1,1:3)) error('For g/A strain, Sys.gFrame and Sys.AFrame must be collinear.'); @@ -65,17 +66,18 @@ % Combine g and A strain rho = CoreSys.gAStrainCorr; % correlation coefficient - for e = CoreSys.nElectrons:-1:1 - lw2_gA{e} = repmat(diag(lw_g(e,:).^2),[1,1,nTransitions]) + lw_A.^2 + ... - 2*rho*repmat(diag(lw_g(e,:)),[1,1,nTransitions]).*lw_A; + for Ne = CoreSys.nElectrons:-1:1 + lw2_gA{Ne} = repmat(diag(lw_g(Ne,:).^2),[1,1,nTransitions]) + lw_A.^2 + ... + 2*rho*repmat(diag(lw_g(Ne,:)),[1,1,nTransitions]).*lw_A; for tr = 1:nTransitions - lw2_gA{e}(:,:,tr) = R_strain2M{e}*lw2_gA{e}(:,:,tr)*R_strain2M{e}.'; + lw2_gA{Ne}(:,:,tr) = R_strain2M{Ne}*lw2_gA{Ne}(:,:,tr)*R_strain2M{Ne}.'; end end + else if usegStrain - for e = CoreSys.nElectrons:-1:1 - lw2_gA{e} = repmat(R_g2M{e}*diag(lw_g(e,:).^2)*R_g2M{e}.',[1,1,nTransitions]); + for Ne = CoreSys.nElectrons:-1:1 + lw2_gA{Ne} = repmat(R_g2M{Ne}*diag(lw_g(Ne,:).^2)*R_g2M{Ne}.',[1,1,nTransitions]); end end end From bb2700d777bc7de16980efd88d4e6d6f4e7c6922 Mon Sep 17 00:00:00 2001 From: fmentink Date: Tue, 1 Apr 2025 16:31:13 -0400 Subject: [PATCH 4/6] updated functions ham_ez, ee, hf and zf --- easyspin/ham_ee.m | 42 ++++++++++++++++++++++++++++++++---------- easyspin/ham_ez.m | 40 ++++++++++++++++++++++++++++++++++++++-- easyspin/ham_hf.m | 23 ++++++++++++++++++++--- easyspin/ham_zf.m | 41 +++++++++++++++++++++++++++++++---------- 4 files changed, 121 insertions(+), 25 deletions(-) diff --git a/easyspin/ham_ee.m b/easyspin/ham_ee.m index 9a73e29e..37dc4faa 100644 --- a/easyspin/ham_ee.m +++ b/easyspin/ham_ee.m @@ -1,4 +1,4 @@ -% ham_ee Electron-electron spin interaction Hamiltonian +% ham_ee Electron-electron spin interaction Hamiltonian % % F = ham_ee(SpinSystem) % F = ham_ee(SpinSystem,eSpins) @@ -19,13 +19,13 @@ % - F: Hamiltonian matrix containing the EEI for % electron spins specified in eSpins. -function F = ham_ee(System,Spins,opt) +function [F,dF] = ham_ee(System,Spins,opt) if nargin==0, help(mfilename); return; end if nargin<1 || nargin>3, error('Wrong number of input arguments!'); end if nargout<0, error('Not enough output arguments.'); end -if nargout>1, error('Too many output arguments.'); end +if nargout>2, error('Too many output arguments.'); end if nargin<3, opt = ''; end if nargin<2, Spins = []; end @@ -49,7 +49,7 @@ % Some error checking on the second input argument if numel(Spins)<2 - error('Spins (2nd argument) must contain at least 2 values!'); + error('Spins (2nd argument) must contain at least 2 values!'); end if any(Spins<1) || any(Spins>System.nElectrons) error('Spins (2nd argument) must contain values between 1 and %d!',System.nElectrons); @@ -80,11 +80,13 @@ eeFrame = []; end +%for conistency accross methods +nStates=System.nStates; % Bilinear coupling term S1*ee*S2 %---------------------------------------------------------------- for iPair = 1:nPairs iCoupling = find(Coupl(iPair)==allPairsIdx); - + % Construct matrix representing coupling tensor if System.fullee J = ee(3*(iCoupling-1)+(1:3),:); @@ -95,25 +97,45 @@ R_M2ee = erot(eeFrame(iCoupling,:)); % mol frame -> ee frame R_ee2M = R_M2ee.'; % ee frame -> mol frame J = R_ee2M*J*R_ee2M.'; + else + R_ee2M = eye(3); end - + + % preparing the derivatives (specific for each electron) + dFdeex = sparse(nStates,nStates); + dFdeey = sparse(nStates,nStates); + dFdeez = sparse(nStates,nStates); + + % preparing the derivatives coefficient + deexM = R_ee2M(:,1)*R_ee2M(:,1).'; % rotate derivative wrt eex to molecular frame + deeyM = R_ee2M(:,2)*R_ee2M(:,2).'; % rotate derivative wrt eey to molecular frame + deezM = R_ee2M(:,3)*R_ee2M(:,3).'; % rotate derivative wrt eez to molecular frame + % Sum up Hamiltonian terms for c1 = 1:3 so1 = sop(sys,[Pairs(iPair,1),c1],'sparse'); for c2 = 1:3 so2 = sop(sys,[Pairs(iPair,2),c2],'sparse'); - F = F + so1*J(c1,c2)*so2; + tempProduct=so1*so2; + F = F + J(c1,c2)*tempProduct; + dFdeex = dFdeex + deexM(c1,c2)*tempProduct; + dFdeey = dFdeey + deeyM(c1,c2)*tempProduct; + dFdeez = dFdeez + deezM(c1,c2)*tempProduct; end end - + + dF{iPair} = {dFdeex,dFdeey,dFdeez}; % derivatives + % Isotropic biquadratic exchange coupling term +ee2*(S1.S2)^2 - %----------------------------------------------------------------- + %----------------------------------------------------------------- if System.ee2(iCoupling)==0, continue; end F2 = 0; for c = 1:3 F2 = F2 + sop(sys,[Pairs(iPair,1),c;Pairs(iPair,2),c],'sparse'); end - F = F + System.ee2(iCoupling)*F2^2; + dFdJ=F2^2; + F = F + System.ee2(iCoupling)*dFdJ; + dF{iPair} = {dFdeex,dFdeey,dFdeez,dFdJ}; % derivatives end F = (F+F')/2; % Hermitianise diff --git a/easyspin/ham_ez.m b/easyspin/ham_ez.m index 20c02be7..4764deba 100644 --- a/easyspin/ham_ez.m +++ b/easyspin/ham_ez.m @@ -31,7 +31,7 @@ if nargin==0, help(mfilename); return; end -if nargout~=1 && nargout~=3, error('Wrong number of output arguments!'); end +if nargout~=1 && nargout~=6, error('Wrong number of output arguments!'); end singleOutput = nargout==1; if singleOutput @@ -105,14 +105,46 @@ R_M2g = erot(ang); % mol frame -> g frame R_g2M = R_M2g.'; % g frame -> mol frame g = R_g2M*g*R_g2M.'; + else + R_g2M=eye(3); end + % preparing the derivatives (specific for each electron) + dMuxxM = sparse(nStates,nStates); + dMuxyM = sparse(nStates,nStates); + dMuxzM = sparse(nStates,nStates); + + dMuyxM = sparse(nStates,nStates); + dMuyyM = sparse(nStates,nStates); + dMuyzM = sparse(nStates,nStates); + + dMuzxM = sparse(nStates,nStates); + dMuzyM = sparse(nStates,nStates); + dMuzzM = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dgxM = R_g2M(:,1)*R_g2M(:,1).'; % rotate derivative wrt gx to molecular frame + dgyM = R_g2M(:,2)*R_g2M(:,2).'; % rotate derivative wrt gy to molecular frame + dgzM = R_g2M(:,3)*R_g2M(:,3).'; % rotate derivative wrt gz to molecular frame + % Build magnetic dipole moment components in MHz/mT for k = 1:3 Sk = sop(spins,[i,k],'sparse'); muxM = muxM + g(1,k)*Sk; muyM = muyM + g(2,k)*Sk; muzM = muzM + g(3,k)*Sk; + + dMuxxM = dMuxxM + dgxM(1,k)*Sk; + dMuyxM = dMuyxM + dgxM(2,k)*Sk; + dMuzxM = dMuzxM + dgxM(3,k)*Sk; + + dMuxyM = dMuxyM + dgyM(1,k)*Sk; + dMuyyM = dMuyyM + dgyM(2,k)*Sk; + dMuzyM = dMuzyM + dgyM(3,k)*Sk; + + dMuxzM = dMuxzM + dgzM(1,k)*Sk; + dMuyzM = dMuyzM + dgzM(2,k)*Sk; + dMuzzM = dMuzzM + dgzM(3,k)*Sk; end end @@ -123,8 +155,12 @@ muxM = full(muxM); muyM = full(muyM); muzM = full(muzM); + + dMuxM = {dMuxxM, dMuyxM, dMuzxM}; + dMuyM = {dMuxyM, dMuyyM, dMuzyM}; + dMuzM = {dMuxzM, dMuyzM, dMuzzM}; end - varargout = {muxM, muyM, muzM}; + varargout = {muxM, muyM, muzM, dMuxM, dMuyM, dMuzM}; else % Return Zeeman Hamiltonian H = -(muxM*B0(1) + muyM*B0(2) + muzM*B0(3)); diff --git a/easyspin/ham_hf.m b/easyspin/ham_hf.m index 2f641eed..45917045 100644 --- a/easyspin/ham_hf.m +++ b/easyspin/ham_hf.m @@ -1,4 +1,4 @@ -% ham_hf Hyperfine interaction Hamiltonian +% ham_hf Hyperfine interaction Hamiltonian % % Hhf = ham_hf(System) % Hhf = ham_hf(System,eSpins) @@ -12,7 +12,7 @@ % % If 'sparse' is given, the matrix is returned in sparse format. -function Hhf = ham_hf(System,elSpins,nucSpins,opt) +function [Hhf,dHhf] = ham_hf(System,elSpins,nucSpins,opt) if nargin==0 help(mfilename); @@ -102,16 +102,33 @@ R_M2A = erot(ang); % mol frame -> A frame R_A2M = R_M2A.'; % A frame -> mol frame A = R_A2M*A*R_A2M.'; + else + R_A2M=eye(3); end + % preparing the derivatives (specific for each electron) + dHhfx = sparse(nStates,nStates); + dHhfy = sparse(nStates,nStates); + dHhfz = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dAxM = R_A2M(:,1)*R_A2M(:,1).'; % rotate derivative wrt Ax to molecular frame + dAyM = R_A2M(:,2)*R_A2M(:,2).'; % rotate derivative wrt Ay to molecular frame + dAzM = R_A2M(:,3)*R_A2M(:,3).'; % rotate derivative wrt Az to molecular frame + % Construct hyperfine Hamiltonian for c1 = 1:3 for c2 = 1:3 - Hhf = Hhf + A(c1,c2)*sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + tempProduct=sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + Hhf = Hhf + A(c1,c2)*tempProduct; + dHhfx = dHhfx + dAxM(c1,c2)*tempProduct; + dHhfy = dHhfy + dAyM(c1,c2)*tempProduct; + dHhfz = dHhfz + dAzM(c1,c2)*tempProduct; end end end % for all specified nuclei + dHhf{eSp,nSp}={dHhfx,dHhfy,dHhfz}; % --> the indices in the derivative may not be needed? end % for all specified electrons Hhf = (Hhf+Hhf')/2; % hermitianise, e.g. guards against small imaginary remainders on the diagonal diff --git a/easyspin/ham_zf.m b/easyspin/ham_zf.m index f7640343..91f3e5af 100644 --- a/easyspin/ham_zf.m +++ b/easyspin/ham_zf.m @@ -1,4 +1,4 @@ -% ham_zf Electronic zero field interaction Hamiltonian +% ham_zf Electronic zero field interaction Hamiltonian % % F = ham_zf(SpinSystem) % F = ham_zf(SpinSystem,Electrons) @@ -7,13 +7,15 @@ % Returns the electronic zero-field interaction (ZFI) % Hamiltonian of the system SpinSystem, in units of MHz. % +% Can retrun the derivative of the ZFI if Sys.Dstrain has been given +% % If the vector Electrons is given, the ZFI of only the % specified electrons is returned (1 is the first, 2 the % second, etc). Otherwise, all electrons are included. % % If 'sparse' is given, the matrix is returned in sparse format. -function H = ham_zf(SpinSystem,idxElectrons,opt) +function [H,dH] = ham_zf(SpinSystem,idxElectrons,opt) if nargin==0, help(mfilename); return; end @@ -35,13 +37,15 @@ error('Electron spin index/indices (2nd argument) out of range!'); end -H = sparse(Sys.nStates,Sys.nStates); +nStates=Sys.nStates; +H = sparse(nStates,nStates); + Spins = Sys.Spins; % Quadratic term S*D*S %------------------------------------------------------------------------------- for iSpin = idxElectrons - + % Prepare full 3x3 D matrix if Sys.fullD D = Sys.D(3*(iSpin-1)+(1:3),:); @@ -53,14 +57,26 @@ continue end - % Transform D tensor to molecular frame + % Transform D tensor to molecular frame and obtain derivatives ang = Sys.DFrame(iSpin,:); if any(ang) R_M2D = erot(ang); % mol frame -> D frame - R_D2M = R_M2D.'; % D frame -> mol frame + R_D2M = R_M2D.'; % D frame -> mol frame D = R_D2M*D*R_D2M.'; + else + R_D2M = eye(3); end + % preparing the derivatives (specific for each electron) + dHdDx = sparse(nStates,nStates); + dHdDy = sparse(nStates,nStates); + dHdDz = sparse(nStates,nStates); + + % preparing the derivatives coefficient + dDxM = R_D2M(:,1)*R_D2M(:,1).'; % rotate derivative wrt Dx to molecular frame + dDyM = R_D2M(:,2)*R_D2M(:,2).'; % rotate derivative wrt Dy to molecular frame + dDzM = R_D2M(:,3)*R_D2M(:,3).'; % rotate derivative wrt Dz to molecular frame + % Construct spin operator matrices for c = 3:-1:1 Sxyz{c} = sop(Spins,[iSpin,c],'sparse'); @@ -69,10 +85,15 @@ % Construct SDS term for c1 = 1:3 for c2 = 1:3 - H = H + D(c1,c2)*(Sxyz{c1}*Sxyz{c2}); + tempProduct=Sxyz{c1}*Sxyz{c2}; % storing the matrix product temporarily + H = H + D(c1,c2)*tempProduct; + dHdDx = dHdDx + dDxM(c1,c2)*tempProduct; + dHdDy = dHdDy + dDyM(c1,c2)*tempProduct; + dHdDz = dHdDz + dDzM(c1,c2)*tempProduct; end end + dH{iSpin} = {dHdDx,dHdDy,dHdDz}; end % Fourth-order terms a and F @@ -110,7 +131,7 @@ % Skip if all rank-k coefficients are zero if all(Bk(iSpin,:)==0), continue; end - + % Apply transformation if non-zero tilt angles are given % (Sys.BFrame is processed from Sys.B?Frame by validatespinsys) tiltang = Sys.BFrame{k}(iSpin,:); @@ -126,13 +147,13 @@ else Bk_M = Bk(iSpin,:); end - + % Build Hamiltonian q = k:-1:-k; for iq = find(Bk_M~=0) H = H + Bk_M(iq)*stev(Spins,[k,q(iq),iSpin],'sparse'); end - + end % for all electron spins specified end % for all tensor ranks From 389b1b798525b2f29965606ebd17dbb0961bd006 Mon Sep 17 00:00:00 2001 From: Stefan Stoll Date: Wed, 2 Apr 2025 14:16:06 -0700 Subject: [PATCH 5/6] ham_hf: edits derivatives code, add tests --- easyspin/ham_hf.m | 13 ++++++++++--- tests/ham_hf_deriv_size.m | 30 ++++++++++++++++++++++++++++++ tests/ham_hf_deriv_value.m | 16 ++++++++++++++++ 3 files changed, 56 insertions(+), 3 deletions(-) create mode 100644 tests/ham_hf_deriv_size.m create mode 100644 tests/ham_hf_deriv_value.m diff --git a/easyspin/ham_hf.m b/easyspin/ham_hf.m index 45917045..57f293cf 100644 --- a/easyspin/ham_hf.m +++ b/easyspin/ham_hf.m @@ -103,7 +103,7 @@ R_A2M = R_M2A.'; % A frame -> mol frame A = R_A2M*A*R_A2M.'; else - R_A2M=eye(3); + R_A2M = eye(3); end % preparing the derivatives (specific for each electron) @@ -119,7 +119,7 @@ % Construct hyperfine Hamiltonian for c1 = 1:3 for c2 = 1:3 - tempProduct=sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + tempProduct = sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); Hhf = Hhf + A(c1,c2)*tempProduct; dHhfx = dHhfx + dAxM(c1,c2)*tempProduct; dHhfy = dHhfy + dAyM(c1,c2)*tempProduct; @@ -127,13 +127,20 @@ end end + dHhf{eSp,nSp} = {dHhfx,dHhfy,dHhfz}; end % for all specified nuclei - dHhf{eSp,nSp}={dHhfx,dHhfy,dHhfz}; % --> the indices in the derivative may not be needed? end % for all specified electrons Hhf = (Hhf+Hhf')/2; % hermitianise, e.g. guards against small imaginary remainders on the diagonal if ~useSparseMatrices Hhf = full(Hhf); % convert sparse to full + for eSp = elSpins + for nSp = nucSpins + for k = 1:3 + dHhf{eSp,nSp}{k} = full(dHhf{eSp,nSp}{k}); + end + end + end end end diff --git a/tests/ham_hf_deriv_size.m b/tests/ham_hf_deriv_size.m new file mode 100644 index 00000000..2dd70c7e --- /dev/null +++ b/tests/ham_hf_deriv_size.m @@ -0,0 +1,30 @@ +function ok = test() + +Sys(1).S = 1/2; +Sys(1).Nucs = '1H'; +Sys(1).A = [1 2 3]; + +Sys(2).S = 1/2; +Sys(2).Nucs = '1H,14N'; +Sys(2).A = [1 2 3; 4 5 6]; + +Sys(3).S = [1/2 1]; +Sys(3).Nucs = '1H'; +Sys(3).A = [1 2 3, 4 5 6]; +Sys(3).J = 100; + +Sys(4).S = [1/2 1]; +Sys(4).Nucs = '1H,13C'; +Sys(4).A = [1 2 3, 4 5 6; 7 8 9, 11 15 19]; +Sys(4).J = 100; + +for k = 1:numel(Sys) + [~,dHhf] = ham_hf(Sys(k)); + elSpins = numel(Sys(k).S); + I = nucspin(Sys(k).Nucs); + nucSpins = numel(I); + ok(k,1) = iscell(dHhf); + ok(k,2) = all(size(dHhf)==[elSpins nucSpins]); + nElements = cellfun(@(x)numel(x),dHhf); + ok(k,3) = all(nElements(:)==3); +end diff --git a/tests/ham_hf_deriv_value.m b/tests/ham_hf_deriv_value.m new file mode 100644 index 00000000..d60c99e9 --- /dev/null +++ b/tests/ham_hf_deriv_value.m @@ -0,0 +1,16 @@ +function ok = test() + +% Test whether hyperfine hamiltonian derivatives are correct, +% for a system with one electron and one nucleus + +Sys.S = 1/2; +Sys.Nucs = '14N'; +Sys.A = [1.435 2.765 3.9876]; + +[Hhf,dHhf] = ham_hf(Sys); + +% Calculate hamiltonian from derivatives +Hhf_2 = Sys.A(1)*dHhf{1}{1} + Sys.A(2)*dHhf{1}{2} + Sys.A(3)*dHhf{1}{3}; + +ok = areequal(Hhf_2,Hhf,1e-10,'abs'); + From 68e64b892c2f0f9401bb6dbca6f808a19893af06 Mon Sep 17 00:00:00 2001 From: fmentink Date: Mon, 7 Apr 2025 11:38:40 -0400 Subject: [PATCH 6/6] updated ham_ez, ham_hf, ham_ee, ham_zf as well as inclusion of tests functions for the derivatives. --- easyspin/ham_ee.m | 29 +++--- easyspin/ham_eeD.m | 146 ------------------------------ easyspin/ham_ez.m | 87 +++++++++--------- easyspin/ham_ezD.m | 173 ------------------------------------ easyspin/ham_hf.m | 40 +++++---- easyspin/ham_hfD.m | 139 ----------------------------- easyspin/ham_zf.m | 36 ++++---- easyspin/ham_zfD.m | 176 ------------------------------------- tests/ham_ee_deriv_size.m | 20 +++++ tests/ham_ee_deriv_value.m | 34 +++++++ tests/ham_ez_deriv_size.m | 28 ++++++ tests/ham_ez_deriv_value.m | 69 +++++++++++++++ tests/ham_hf_deriv_value.m | 22 ++++- tests/ham_zf_deriv_size.m | 20 +++++ tests/ham_zf_deriv_value.m | 32 +++++++ 15 files changed, 322 insertions(+), 729 deletions(-) delete mode 100644 easyspin/ham_eeD.m delete mode 100644 easyspin/ham_ezD.m delete mode 100644 easyspin/ham_hfD.m delete mode 100644 easyspin/ham_zfD.m create mode 100644 tests/ham_ee_deriv_size.m create mode 100644 tests/ham_ee_deriv_value.m create mode 100644 tests/ham_ez_deriv_size.m create mode 100644 tests/ham_ez_deriv_value.m create mode 100644 tests/ham_zf_deriv_size.m create mode 100644 tests/ham_zf_deriv_value.m diff --git a/easyspin/ham_ee.m b/easyspin/ham_ee.m index 37dc4faa..ea2f2f4f 100644 --- a/easyspin/ham_ee.m +++ b/easyspin/ham_ee.m @@ -1,11 +1,11 @@ % ham_ee Electron-electron spin interaction Hamiltonian % -% F = ham_ee(SpinSystem) -% F = ham_ee(SpinSystem,eSpins) -% F = ham_ee(SpinSystem,eSpins,'sparse') +% [F, dF] = ham_ee(SpinSystem) +% [F, dF] = ham_ee(SpinSystem,eSpins) +% [F, dF] = ham_ee(SpinSystem,eSpins,'sparse') % % Returns the electron-electron spin interaction (EEI) -% Hamiltonian, in MHz. +% Hamiltonian, in MHz and its derivatives. % % Input: % - SpinSystem: Spin system structure. EEI @@ -18,6 +18,8 @@ % Output: % - F: Hamiltonian matrix containing the EEI for % electron spins specified in eSpins. +% - dF: Derivative of the Hamiltonian matrix containing the EEI for +% electron spins specified in eSpins. function [F,dF] = ham_ee(System,Spins,opt) @@ -113,9 +115,17 @@ % Sum up Hamiltonian terms for c1 = 1:3 - so1 = sop(sys,[Pairs(iPair,1),c1],'sparse'); + if ~useSparseMatrices % sparse -> full + so1 = sop(sys,[Pairs(iPair,1),c1]); + else + so1 = sop(sys,[Pairs(iPair,1),c1],'sparse'); + end for c2 = 1:3 - so2 = sop(sys,[Pairs(iPair,2),c2],'sparse'); + if ~useSparseMatrices % sparse -> full + so2 = sop(sys,[Pairs(iPair,2),c2]); + else + so2 = sop(sys,[Pairs(iPair,2),c2],'sparse'); + end tempProduct=so1*so2; F = F + J(c1,c2)*tempProduct; dFdeex = dFdeex + deexM(c1,c2)*tempProduct; @@ -133,14 +143,13 @@ for c = 1:3 F2 = F2 + sop(sys,[Pairs(iPair,1),c;Pairs(iPair,2),c],'sparse'); end + if ~useSparseMatrices % sparse -> full + F2=full(F2); + end dFdJ=F2^2; F = F + System.ee2(iCoupling)*dFdJ; dF{iPair} = {dFdeex,dFdeey,dFdeez,dFdJ}; % derivatives end F = (F+F')/2; % Hermitianise -if ~useSparseMatrices - F = full(F); % sparse -> full -end - end diff --git a/easyspin/ham_eeD.m b/easyspin/ham_eeD.m deleted file mode 100644 index 37dc4faa..00000000 --- a/easyspin/ham_eeD.m +++ /dev/null @@ -1,146 +0,0 @@ -% ham_ee Electron-electron spin interaction Hamiltonian -% -% F = ham_ee(SpinSystem) -% F = ham_ee(SpinSystem,eSpins) -% F = ham_ee(SpinSystem,eSpins,'sparse') -% -% Returns the electron-electron spin interaction (EEI) -% Hamiltonian, in MHz. -% -% Input: -% - SpinSystem: Spin system structure. EEI -% parameters are in the ee, eeFrame, and ee2 fields. -% - eSpins: If given, specifies electron spins -% for which the EEI should be computed. If -% absent, all electrons are included. -% - 'sparse': If given, the matrix is returned in sparse format. - -% Output: -% - F: Hamiltonian matrix containing the EEI for -% electron spins specified in eSpins. - -function [F,dF] = ham_ee(System,Spins,opt) - -if nargin==0, help(mfilename); return; end - -if nargin<1 || nargin>3, error('Wrong number of input arguments!'); end -if nargout<0, error('Not enough output arguments.'); end -if nargout>2, error('Too many output arguments.'); end - -if nargin<3, opt = ''; end -if nargin<2, Spins = []; end -if ~ischar(opt) - error('Third input must be a string, ''sparse''.'); -end -useSparseMatrices = strcmp(opt,'sparse'); - -[System,err] = validatespinsys(System); -error(err); - -sys = spinvec(System); -n = prod(2*sys+1); - -% Special cases: only one spins, ee not given or all zero -F = sparse(n,n); -if (System.nElectrons==1), return; end -if ~any(System.ee(:)) && ~any(System.ee2(:)), return; end - -if isempty(Spins), Spins = 1:System.nElectrons; end - -% Some error checking on the second input argument -if numel(Spins)<2 - error('Spins (2nd argument) must contain at least 2 values!'); -end -if any(Spins<1) || any(Spins>System.nElectrons) - error('Spins (2nd argument) must contain values between 1 and %d!',System.nElectrons); -end -if numel(unique(Spins))~=numel(Spins) - error('Spins (2nd argument) contains double entries!'); -end - -F = sparse(n,n); - -% Compile list of wanted interactions -Spins = sort(Spins); -[idx1,idx2] = find(tril(ones(numel(Spins)),-1)); -idx = [idx2,idx1]; - -Pairs = Spins(idx); -nPairs = size(Pairs,1); -Coupl = Pairs(:,1) + (Pairs(:,2)-1)*System.nElectrons; - -% Compile list of all spin pairs -[e2,e1] = find(tril(ones(System.nElectrons),-1)); -allPairsIdx = e1 + (e2-1)*System.nElectrons; - -ee = System.ee; -if isfield(System,'eeFrame') - eeFrame = System.eeFrame; -else - eeFrame = []; -end - -%for conistency accross methods -nStates=System.nStates; -% Bilinear coupling term S1*ee*S2 -%---------------------------------------------------------------- -for iPair = 1:nPairs - iCoupling = find(Coupl(iPair)==allPairsIdx); - - % Construct matrix representing coupling tensor - if System.fullee - J = ee(3*(iCoupling-1)+(1:3),:); - else - J = diag(ee(iCoupling,:)); - end - if ~isempty(eeFrame) && any(eeFrame(iCoupling,:)) - R_M2ee = erot(eeFrame(iCoupling,:)); % mol frame -> ee frame - R_ee2M = R_M2ee.'; % ee frame -> mol frame - J = R_ee2M*J*R_ee2M.'; - else - R_ee2M = eye(3); - end - - % preparing the derivatives (specific for each electron) - dFdeex = sparse(nStates,nStates); - dFdeey = sparse(nStates,nStates); - dFdeez = sparse(nStates,nStates); - - % preparing the derivatives coefficient - deexM = R_ee2M(:,1)*R_ee2M(:,1).'; % rotate derivative wrt eex to molecular frame - deeyM = R_ee2M(:,2)*R_ee2M(:,2).'; % rotate derivative wrt eey to molecular frame - deezM = R_ee2M(:,3)*R_ee2M(:,3).'; % rotate derivative wrt eez to molecular frame - - % Sum up Hamiltonian terms - for c1 = 1:3 - so1 = sop(sys,[Pairs(iPair,1),c1],'sparse'); - for c2 = 1:3 - so2 = sop(sys,[Pairs(iPair,2),c2],'sparse'); - tempProduct=so1*so2; - F = F + J(c1,c2)*tempProduct; - dFdeex = dFdeex + deexM(c1,c2)*tempProduct; - dFdeey = dFdeey + deeyM(c1,c2)*tempProduct; - dFdeez = dFdeez + deezM(c1,c2)*tempProduct; - end - end - - dF{iPair} = {dFdeex,dFdeey,dFdeez}; % derivatives - - % Isotropic biquadratic exchange coupling term +ee2*(S1.S2)^2 - %----------------------------------------------------------------- - if System.ee2(iCoupling)==0, continue; end - F2 = 0; - for c = 1:3 - F2 = F2 + sop(sys,[Pairs(iPair,1),c;Pairs(iPair,2),c],'sparse'); - end - dFdJ=F2^2; - F = F + System.ee2(iCoupling)*dFdJ; - dF{iPair} = {dFdeex,dFdeey,dFdeez,dFdJ}; % derivatives -end - -F = (F+F')/2; % Hermitianise -if ~useSparseMatrices - F = full(F); % sparse -> full -end - -end diff --git a/easyspin/ham_ez.m b/easyspin/ham_ez.m index 4764deba..0fd0e2f9 100644 --- a/easyspin/ham_ez.m +++ b/easyspin/ham_ez.m @@ -6,9 +6,14 @@ % [mux, muy, muz] = ham_ez(SpinSystem) % [mux, muy, muz] = ham_ez(SpinSystem, eSpins) % [mux, muy, muz] = ham_ez(SpinSystem, eSpins, 'sparse') +% [mux, muy, muz, dmudgx, dmudgy, dmudgz] = ham_ez(SpinSystem) +% [mux, muy, muz, dmudgx, dmudgy, dmudgz] = ham_ez(SpinSystem, eSpins) +% [mux, muy, muz, dmudgx, dmudgy, dmudgz] = ham_ez(SpinSystem, eSpins, 'sparse') % % Returns the electron Zeeman interaction Hamiltonian matrix for -% the electron spins 'eSpins' of the spin system 'SpinSystem'. +% the electron spins 'eSpins' of the spin system 'SpinSystem'. It can +% return also the magnetic moment and its derivatives with respect to gx, +% gy, gz in the molecular frame. % % Input: % - SpinSystem: Spin system structure. @@ -25,13 +30,17 @@ % To get the full electron Zeeman Hamiltonian, use % H = -(mux*B(1)+muy*B(2)+muz*B(3)), where B is the magnetic field vector % in mT. +% - dmudgx, dmudgy, dmudgz, derivatives of the magnetic dipole moment +% with respect to dmudgi= dmu/dgi where i=x,y,z stored as a cell index as {spin,i} % - H: Electron Zeeman Hamiltonian matrix. function varargout = ham_ez(SpinSystem,varargin) if nargin==0, help(mfilename); return; end -if nargout~=1 && nargout~=6, error('Wrong number of output arguments!'); end +if nargout~=1 && nargout~=3 && nargout~=6 + error('Wrong number of output arguments!'); +end singleOutput = nargout==1; if singleOutput @@ -88,19 +97,20 @@ % Calculate prefactors pre = -bmagn/planck*Sys.g; % Hz/T pre = pre/1e9; % Hz/T -> GHz/T = MHz/mT +preDer = -bmagn/planck/1e9; % Loop over all electron spins -for i = eSpins +for eSp = eSpins % Get full g matrix if Sys.fullg - g = pre((i-1)*3+(1:3),:); + g = pre((eSp-1)*3+(1:3),:); else - g = diag(pre(i,:)); + g = diag(pre(eSp,:)); end % Transform g matrix to molecular frame - ang = Sys.gFrame(i,:); + ang = Sys.gFrame(eSp,:); if any(ang) R_M2g = erot(ang); % mol frame -> g frame R_g2M = R_M2g.'; % g frame -> mol frame @@ -109,64 +119,47 @@ R_g2M=eye(3); end - % preparing the derivatives (specific for each electron) - dMuxxM = sparse(nStates,nStates); - dMuxyM = sparse(nStates,nStates); - dMuxzM = sparse(nStates,nStates); - - dMuyxM = sparse(nStates,nStates); - dMuyyM = sparse(nStates,nStates); - dMuyzM = sparse(nStates,nStates); - - dMuzxM = sparse(nStates,nStates); - dMuzyM = sparse(nStates,nStates); - dMuzzM = sparse(nStates,nStates); + % preparing the derivatives of the magnetic moment + for index=1:3 + dmuMgx{index} = sparse(nStates,nStates); + dmuMgy{index} = sparse(nStates,nStates); + dmuMgz{index} = sparse(nStates,nStates); + end % preparing the derivatives coefficient - dgxM = R_g2M(:,1)*R_g2M(:,1).'; % rotate derivative wrt gx to molecular frame - dgyM = R_g2M(:,2)*R_g2M(:,2).'; % rotate derivative wrt gy to molecular frame - dgzM = R_g2M(:,3)*R_g2M(:,3).'; % rotate derivative wrt gz to molecular frame + dgxM = preDer*R_g2M(:,1)*R_g2M(:,1).'; % rotate derivative wrt gx to molecular frame (and scale with preDer) + dgyM = preDer*R_g2M(:,2)*R_g2M(:,2).'; % rotate derivative wrt gy to molecular frame (and scale with preDer) + dgzM = preDer*R_g2M(:,3)*R_g2M(:,3).'; % rotate derivative wrt gz to molecular frame (and scale with preDer) % Build magnetic dipole moment components in MHz/mT for k = 1:3 - Sk = sop(spins,[i,k],'sparse'); + if ~useSparseMatrices + Sk = sop(spins,[eSp,k]); + else + Sk = sop(spins,[eSp,k],'sparse'); + end muxM = muxM + g(1,k)*Sk; muyM = muyM + g(2,k)*Sk; muzM = muzM + g(3,k)*Sk; - dMuxxM = dMuxxM + dgxM(1,k)*Sk; - dMuyxM = dMuyxM + dgxM(2,k)*Sk; - dMuzxM = dMuzxM + dgxM(3,k)*Sk; - - dMuxyM = dMuxyM + dgyM(1,k)*Sk; - dMuyyM = dMuyyM + dgyM(2,k)*Sk; - dMuzyM = dMuzyM + dgyM(3,k)*Sk; - - dMuxzM = dMuxzM + dgzM(1,k)*Sk; - dMuyzM = dMuyzM + dgzM(2,k)*Sk; - dMuzzM = dMuzzM + dgzM(3,k)*Sk; + for index=1:3 + dmuMgx{index} = dmuMgx{index} + dgxM(index,k)*Sk; + dmuMgy{index} = dmuMgy{index} + dgyM(index,k)*Sk; + dmuMgz{index} = dmuMgz{index} + dgzM(index,k)*Sk; + end end + dmuMdgx{eSp} = {dmuMgx{1}, dmuMgx{2}, dmuMgx{3}}; + dmuMdgy{eSp} = {dmuMgy{1}, dmuMgy{2}, dmuMgy{3}}; + dmuMdgz{eSp} = {dmuMgz{1}, dmuMgz{2}, dmuMgz{3}}; end if isempty(B0) - % Return magnetic dipole moment components - if ~useSparseMatrices - muxM = full(muxM); - muyM = full(muyM); - muzM = full(muzM); - - dMuxM = {dMuxxM, dMuyxM, dMuzxM}; - dMuyM = {dMuxyM, dMuyyM, dMuzyM}; - dMuzM = {dMuxzM, dMuyzM, dMuzzM}; - end - varargout = {muxM, muyM, muzM, dMuxM, dMuyM, dMuzM}; + % Return magnetic dipole moment components and derivatives + varargout = {muxM, muyM, muzM, dmuMdgx, dmuMdgy, dmuMdgz}; else % Return Zeeman Hamiltonian H = -(muxM*B0(1) + muyM*B0(2) + muzM*B0(3)); - if ~useSparseMatrices - H = full(H); - end varargout = {H}; end diff --git a/easyspin/ham_ezD.m b/easyspin/ham_ezD.m deleted file mode 100644 index 4764deba..00000000 --- a/easyspin/ham_ezD.m +++ /dev/null @@ -1,173 +0,0 @@ -% ham_ez Electron Zeeman interaction Hamiltonian -% -% H = ham_ez(SpinSystem, B) -% H = ham_ez(SpinSystem, B, eSpins) -% H = ham_ez(SpinSystem, B, eSpins, 'sparse') -% [mux, muy, muz] = ham_ez(SpinSystem) -% [mux, muy, muz] = ham_ez(SpinSystem, eSpins) -% [mux, muy, muz] = ham_ez(SpinSystem, eSpins, 'sparse') -% -% Returns the electron Zeeman interaction Hamiltonian matrix for -% the electron spins 'eSpins' of the spin system 'SpinSystem'. -% -% Input: -% - SpinSystem: Spin system structure. -% - B: Magnetic field vector, in millitesla. -% - eSpins: Vector of indices for electron spins to include. If eSpins is -% omitted, all electron spins are included. -% - 'sparse': If given, results returned in sparse format. -% -% Output: -% - mux, muy, muz: Components of the magnetic dipole moment operator -% for the selected electron spins as defined by mui=-d(H)/d(B_i) -% i=x,y,z where B_i are the cartesian components of -% the external field in the molecular frame. Units are MHz/mT = GHz/T. -% To get the full electron Zeeman Hamiltonian, use -% H = -(mux*B(1)+muy*B(2)+muz*B(3)), where B is the magnetic field vector -% in mT. -% - H: Electron Zeeman Hamiltonian matrix. - -function varargout = ham_ez(SpinSystem,varargin) - -if nargin==0, help(mfilename); return; end - -if nargout~=1 && nargout~=6, error('Wrong number of output arguments!'); end -singleOutput = nargout==1; - -if singleOutput - if nargin<1 || nargin>4 - error('Wrong number of input arguments!'); - end - if nargin<2 - error('Field vector (second input, in mT) is missing.') - else - B0 = varargin{1}; - end - if nargin<3, eSpins = []; else, eSpins = varargin{2}; end - if nargin<4, opt = ''; else, opt = varargin{3}; end -else - B0 = []; - if nargin<2, eSpins = []; else, eSpins = varargin{1}; end - if nargin<3, opt = ''; else, opt = varargin{2}; end -end - -if ~ischar(opt) - error('Last input must be a string, ''sparse''.'); -end -useSparseMatrices = strcmp(opt,'sparse'); - -% Validate spin system -[Sys,err] = validatespinsys(SpinSystem); -error(err); - -% Vector of all spin quantum numbers -spins = Sys.Spins; -nElectrons = Sys.nElectrons; - -% No 'eSpins' specified -> use all -if isempty(eSpins), eSpins = 1:nElectrons; end - -% Validate magnetic field if given -if ~isempty(B0) - if numel(B0)~=3 - error('Magnetic field vector (2nd input) must be a 3-element array.'); - end -end - -% Validate third argument (eSpins) -if any(eSpins<1) || any(eSpins>nElectrons) - error('Electron spin indices (2nd input argument) invalid!'); -end - -% Initialize Zeeman interaction component matrices to zero -nStates = Sys.nStates; -muxM = sparse(nStates,nStates); -muyM = sparse(nStates,nStates); -muzM = sparse(nStates,nStates); - -% Calculate prefactors -pre = -bmagn/planck*Sys.g; % Hz/T -pre = pre/1e9; % Hz/T -> GHz/T = MHz/mT - -% Loop over all electron spins -for i = eSpins - - % Get full g matrix - if Sys.fullg - g = pre((i-1)*3+(1:3),:); - else - g = diag(pre(i,:)); - end - - % Transform g matrix to molecular frame - ang = Sys.gFrame(i,:); - if any(ang) - R_M2g = erot(ang); % mol frame -> g frame - R_g2M = R_M2g.'; % g frame -> mol frame - g = R_g2M*g*R_g2M.'; - else - R_g2M=eye(3); - end - - % preparing the derivatives (specific for each electron) - dMuxxM = sparse(nStates,nStates); - dMuxyM = sparse(nStates,nStates); - dMuxzM = sparse(nStates,nStates); - - dMuyxM = sparse(nStates,nStates); - dMuyyM = sparse(nStates,nStates); - dMuyzM = sparse(nStates,nStates); - - dMuzxM = sparse(nStates,nStates); - dMuzyM = sparse(nStates,nStates); - dMuzzM = sparse(nStates,nStates); - - % preparing the derivatives coefficient - dgxM = R_g2M(:,1)*R_g2M(:,1).'; % rotate derivative wrt gx to molecular frame - dgyM = R_g2M(:,2)*R_g2M(:,2).'; % rotate derivative wrt gy to molecular frame - dgzM = R_g2M(:,3)*R_g2M(:,3).'; % rotate derivative wrt gz to molecular frame - - % Build magnetic dipole moment components in MHz/mT - for k = 1:3 - Sk = sop(spins,[i,k],'sparse'); - muxM = muxM + g(1,k)*Sk; - muyM = muyM + g(2,k)*Sk; - muzM = muzM + g(3,k)*Sk; - - dMuxxM = dMuxxM + dgxM(1,k)*Sk; - dMuyxM = dMuyxM + dgxM(2,k)*Sk; - dMuzxM = dMuzxM + dgxM(3,k)*Sk; - - dMuxyM = dMuxyM + dgyM(1,k)*Sk; - dMuyyM = dMuyyM + dgyM(2,k)*Sk; - dMuzyM = dMuzyM + dgyM(3,k)*Sk; - - dMuxzM = dMuxzM + dgzM(1,k)*Sk; - dMuyzM = dMuyzM + dgzM(2,k)*Sk; - dMuzzM = dMuzzM + dgzM(3,k)*Sk; - end - -end - -if isempty(B0) - % Return magnetic dipole moment components - if ~useSparseMatrices - muxM = full(muxM); - muyM = full(muyM); - muzM = full(muzM); - - dMuxM = {dMuxxM, dMuyxM, dMuzxM}; - dMuyM = {dMuxyM, dMuyyM, dMuzyM}; - dMuzM = {dMuxzM, dMuyzM, dMuzzM}; - end - varargout = {muxM, muyM, muzM, dMuxM, dMuyM, dMuzM}; -else - % Return Zeeman Hamiltonian - H = -(muxM*B0(1) + muyM*B0(2) + muzM*B0(3)); - if ~useSparseMatrices - H = full(H); - end - varargout = {H}; -end - -end diff --git a/easyspin/ham_hf.m b/easyspin/ham_hf.m index 57f293cf..e877a6d8 100644 --- a/easyspin/ham_hf.m +++ b/easyspin/ham_hf.m @@ -1,11 +1,11 @@ -% ham_hf Hyperfine interaction Hamiltonian +% ham_hf Hyperfine interaction Hamiltonian and derivative % -% Hhf = ham_hf(System) -% Hhf = ham_hf(System,eSpins) -% Hhf = ham_hf(System,eSpins,nSpins) -% Hhf = ham_hf(System,eSpins,nSpins,'sparse') +% [Hhf,dHhf] = ham_hf(System) +% [Hhf,dHhf] = ham_hf(System,eSpins) +% [Hhf,dHhf] = ham_hf(System,eSpins,nSpins) +% [Hhf,dHhf] = ham_hf(System,eSpins,nSpins,'sparse') % -% Returns the hyperfine interaction Hamiltonian (in units of MHz) between +% Returns the hyperfine interaction Hamiltonian (in units of MHz) and its derivative between % electron spins 'eSpins' and nuclear spins 'nSpins' of the system % 'System'. eSpins=1 is the first electron spins, nSpins=1 is the first % nuclear spin. @@ -119,28 +119,32 @@ % Construct hyperfine Hamiltonian for c1 = 1:3 for c2 = 1:3 - tempProduct = sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + if ~useSparseMatrices + tempProduct = sop(SpinVec,[eSp c1; nElectrons+nSp c2]); + else + tempProduct = sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); + end Hhf = Hhf + A(c1,c2)*tempProduct; dHhfx = dHhfx + dAxM(c1,c2)*tempProduct; dHhfy = dHhfy + dAyM(c1,c2)*tempProduct; dHhfz = dHhfz + dAzM(c1,c2)*tempProduct; end end - + % Return the derivative of the hyperfine coupling for each electron-nuclear spin pair dHhf{eSp,nSp} = {dHhfx,dHhfy,dHhfz}; end % for all specified nuclei end % for all specified electrons Hhf = (Hhf+Hhf')/2; % hermitianise, e.g. guards against small imaginary remainders on the diagonal -if ~useSparseMatrices - Hhf = full(Hhf); % convert sparse to full - for eSp = elSpins - for nSp = nucSpins - for k = 1:3 - dHhf{eSp,nSp}{k} = full(dHhf{eSp,nSp}{k}); - end - end - end -end +% if ~useSparseMatrices +% Hhf = full(Hhf); % convert sparse to full +% for eSp = elSpins +% for nSp = nucSpins +% for k = 1:3 +% dHhf{eSp,nSp}{k} = full(dHhf{eSp,nSp}{k}); +% end +% end +% end +% end end diff --git a/easyspin/ham_hfD.m b/easyspin/ham_hfD.m deleted file mode 100644 index 45917045..00000000 --- a/easyspin/ham_hfD.m +++ /dev/null @@ -1,139 +0,0 @@ -% ham_hf Hyperfine interaction Hamiltonian -% -% Hhf = ham_hf(System) -% Hhf = ham_hf(System,eSpins) -% Hhf = ham_hf(System,eSpins,nSpins) -% Hhf = ham_hf(System,eSpins,nSpins,'sparse') -% -% Returns the hyperfine interaction Hamiltonian (in units of MHz) between -% electron spins 'eSpins' and nuclear spins 'nSpins' of the system -% 'System'. eSpins=1 is the first electron spins, nSpins=1 is the first -% nuclear spin. -% -% If 'sparse' is given, the matrix is returned in sparse format. - -function [Hhf,dHhf] = ham_hf(System,elSpins,nucSpins,opt) - -if nargin==0 - help(mfilename); - return -end - -switch nargin - case 1 - elSpins = []; - nucSpins = []; - opt = ''; - case 2 - nucSpins = []; - opt = ''; - case 3 - opt = ''; - case 4 - otherwise - error('Incorrect number of input arguments.') -end - -if ~ischar(opt) - error('Fourth input must be a string, ''sparse''.'); -end -useSparseMatrices = strcmp(opt,'sparse'); - -% Validate spin system. -[Sys,err] = validatespinsys(System); -error(err); - -% Get spin vector and space dimension -SpinVec = Sys.Spins; -nStates = Sys.nStates; -nElectrons = Sys.nElectrons; -nNuclei = Sys.nNuclei; - -Hhf = sparse(nStates,nStates); % sparse zero matrix - -% Special case: no nuclei present, so no hyperfine -if nNuclei==0 - if ~useSparseMatrices - Hhf = full(Hhf); - end - return -end - -% Get electron spin list -if isempty(elSpins) - elSpins = 1:nElectrons; -else - if any(elSpins<1) || any(elSpins>nElectrons) - error('Electron spins (2nd argument) contains out-of-range values!'); - end -end - -% Get electron spin list -if isempty(nucSpins) - nucSpins = 1:nNuclei; -else - if any(nucSpins<1) || any(nucSpins>nNuclei) - error('Nuclear spins (3rd argument) contains out-of-range values!'); - end -end - -fullAMatrix = size(System.A,1)>nNuclei; - -% Generate Hamiltonian for hyperfine interaction -for eSp = elSpins - eidx = (eSp-1)*3+(1:3); - for nSp = nucSpins - if Sys.I(nSp)==0, continue; end - - % Construct full hyperfine matrix - if fullAMatrix - A = Sys.A((nSp-1)*3+(1:3),eidx); - else - A = diag(Sys.A(nSp,eidx)); - end - - if ~any(A(:)) - continue - end - - % Transform matrix into molecular frame representation - ang = Sys.AFrame(nSp,eidx); - if any(ang) - R_M2A = erot(ang); % mol frame -> A frame - R_A2M = R_M2A.'; % A frame -> mol frame - A = R_A2M*A*R_A2M.'; - else - R_A2M=eye(3); - end - - % preparing the derivatives (specific for each electron) - dHhfx = sparse(nStates,nStates); - dHhfy = sparse(nStates,nStates); - dHhfz = sparse(nStates,nStates); - - % preparing the derivatives coefficient - dAxM = R_A2M(:,1)*R_A2M(:,1).'; % rotate derivative wrt Ax to molecular frame - dAyM = R_A2M(:,2)*R_A2M(:,2).'; % rotate derivative wrt Ay to molecular frame - dAzM = R_A2M(:,3)*R_A2M(:,3).'; % rotate derivative wrt Az to molecular frame - - % Construct hyperfine Hamiltonian - for c1 = 1:3 - for c2 = 1:3 - tempProduct=sop(SpinVec,[eSp c1; nElectrons+nSp c2],'sparse'); - Hhf = Hhf + A(c1,c2)*tempProduct; - dHhfx = dHhfx + dAxM(c1,c2)*tempProduct; - dHhfy = dHhfy + dAyM(c1,c2)*tempProduct; - dHhfz = dHhfz + dAzM(c1,c2)*tempProduct; - end - end - - end % for all specified nuclei - dHhf{eSp,nSp}={dHhfx,dHhfy,dHhfz}; % --> the indices in the derivative may not be needed? -end % for all specified electrons - -Hhf = (Hhf+Hhf')/2; % hermitianise, e.g. guards against small imaginary remainders on the diagonal -if ~useSparseMatrices - Hhf = full(Hhf); % convert sparse to full -end - -end diff --git a/easyspin/ham_zf.m b/easyspin/ham_zf.m index 91f3e5af..f95b67ce 100644 --- a/easyspin/ham_zf.m +++ b/easyspin/ham_zf.m @@ -15,7 +15,7 @@ % % If 'sparse' is given, the matrix is returned in sparse format. -function [H,dH] = ham_zf(SpinSystem,idxElectrons,opt) +function [Hzf,dHzf] = ham_zf(SpinSystem,idxElectrons,opt) if nargin==0, help(mfilename); return; end @@ -27,7 +27,7 @@ if ~ischar(opt) error('Third input must be a string, ''sparse''.'); end -sparseResult = strcmp(opt,'sparse'); +useSparseMatrices = strcmp(opt,'sparse'); if isempty(idxElectrons) idxElectrons = 1:Sys.nElectrons; @@ -38,7 +38,7 @@ end nStates=Sys.nStates; -H = sparse(nStates,nStates); +Hzf = sparse(nStates,nStates); Spins = Sys.Spins; @@ -68,9 +68,9 @@ end % preparing the derivatives (specific for each electron) - dHdDx = sparse(nStates,nStates); - dHdDy = sparse(nStates,nStates); - dHdDz = sparse(nStates,nStates); + dHzfdDx = sparse(nStates,nStates); + dHzfdDy = sparse(nStates,nStates); + dHzfdDz = sparse(nStates,nStates); % preparing the derivatives coefficient dDxM = R_D2M(:,1)*R_D2M(:,1).'; % rotate derivative wrt Dx to molecular frame @@ -79,21 +79,25 @@ % Construct spin operator matrices for c = 3:-1:1 - Sxyz{c} = sop(Spins,[iSpin,c],'sparse'); + if ~useSparseMatrices + Sxyz{c} = sop(Spins,[iSpin,c]); + else + Sxyz{c} = sop(Spins,[iSpin,c],'sparse'); + end end % Construct SDS term for c1 = 1:3 for c2 = 1:3 tempProduct=Sxyz{c1}*Sxyz{c2}; % storing the matrix product temporarily - H = H + D(c1,c2)*tempProduct; - dHdDx = dHdDx + dDxM(c1,c2)*tempProduct; - dHdDy = dHdDy + dDyM(c1,c2)*tempProduct; - dHdDz = dHdDz + dDzM(c1,c2)*tempProduct; + Hzf = Hzf + D(c1,c2)*tempProduct; + dHzfdDx = dHzfdDx + dDxM(c1,c2)*tempProduct; + dHzfdDy = dHzfdDy + dDyM(c1,c2)*tempProduct; + dHzfdDz = dHzfdDz + dDzM(c1,c2)*tempProduct; end end - dH{iSpin} = {dHdDx,dHdDy,dHdDz}; + dHzf{iSpin} = {dHzfdDx,dHzfdDy,dHzfdDz}; end % Fourth-order terms a and F @@ -151,18 +155,14 @@ % Build Hamiltonian q = k:-1:-k; for iq = find(Bk_M~=0) - H = H + Bk_M(iq)*stev(Spins,[k,q(iq),iSpin],'sparse'); + Hzf = Hzf + Bk_M(iq)*stev(Spins,[k,q(iq),iSpin],'sparse'); end end % for all electron spins specified end % for all tensor ranks -H = (H+H')/2; % Hermitianize -if ~sparseResult - H = full(H); -end - +Hzf = (Hzf+Hzf')/2; % Hermitianize end %=============================================================================== diff --git a/easyspin/ham_zfD.m b/easyspin/ham_zfD.m deleted file mode 100644 index 91f3e5af..00000000 --- a/easyspin/ham_zfD.m +++ /dev/null @@ -1,176 +0,0 @@ -% ham_zf Electronic zero field interaction Hamiltonian -% -% F = ham_zf(SpinSystem) -% F = ham_zf(SpinSystem,Electrons) -% F = ham_zf(SpinSystem,Electrons,'sparse') -% -% Returns the electronic zero-field interaction (ZFI) -% Hamiltonian of the system SpinSystem, in units of MHz. -% -% Can retrun the derivative of the ZFI if Sys.Dstrain has been given -% -% If the vector Electrons is given, the ZFI of only the -% specified electrons is returned (1 is the first, 2 the -% second, etc). Otherwise, all electrons are included. -% -% If 'sparse' is given, the matrix is returned in sparse format. - -function [H,dH] = ham_zf(SpinSystem,idxElectrons,opt) - -if nargin==0, help(mfilename); return; end - -[Sys,err] = validatespinsys(SpinSystem); -error(err); - -if nargin<2, idxElectrons = []; end -if nargin<3, opt = ''; end -if ~ischar(opt) - error('Third input must be a string, ''sparse''.'); -end -sparseResult = strcmp(opt,'sparse'); - -if isempty(idxElectrons) - idxElectrons = 1:Sys.nElectrons; -end - -if any(idxElectrons>Sys.nElectrons) || any(idxElectrons<1) - error('Electron spin index/indices (2nd argument) out of range!'); -end - -nStates=Sys.nStates; -H = sparse(nStates,nStates); - -Spins = Sys.Spins; - -% Quadratic term S*D*S -%------------------------------------------------------------------------------- -for iSpin = idxElectrons - - % Prepare full 3x3 D matrix - if Sys.fullD - D = Sys.D(3*(iSpin-1)+(1:3),:); - else - D = diag(Sys.D(iSpin,:)); - end - - if ~any(D(:)) - continue - end - - % Transform D tensor to molecular frame and obtain derivatives - ang = Sys.DFrame(iSpin,:); - if any(ang) - R_M2D = erot(ang); % mol frame -> D frame - R_D2M = R_M2D.'; % D frame -> mol frame - D = R_D2M*D*R_D2M.'; - else - R_D2M = eye(3); - end - - % preparing the derivatives (specific for each electron) - dHdDx = sparse(nStates,nStates); - dHdDy = sparse(nStates,nStates); - dHdDz = sparse(nStates,nStates); - - % preparing the derivatives coefficient - dDxM = R_D2M(:,1)*R_D2M(:,1).'; % rotate derivative wrt Dx to molecular frame - dDyM = R_D2M(:,2)*R_D2M(:,2).'; % rotate derivative wrt Dy to molecular frame - dDzM = R_D2M(:,3)*R_D2M(:,3).'; % rotate derivative wrt Dz to molecular frame - - % Construct spin operator matrices - for c = 3:-1:1 - Sxyz{c} = sop(Spins,[iSpin,c],'sparse'); - end - - % Construct SDS term - for c1 = 1:3 - for c2 = 1:3 - tempProduct=Sxyz{c1}*Sxyz{c2}; % storing the matrix product temporarily - H = H + D(c1,c2)*tempProduct; - dHdDx = dHdDx + dDxM(c1,c2)*tempProduct; - dHdDy = dHdDy + dDyM(c1,c2)*tempProduct; - dHdDz = dHdDz + dDzM(c1,c2)*tempProduct; - end - end - - dH{iSpin} = {dHdDx,dHdDy,dHdDz}; -end - -% Fourth-order terms a and F -%----------------------------------------------------------------------------- -% Spin system fields: a and F - -% Abragam/Bleaney, pp. 142, 437 -% Bleaney/Trenam, Proc.Roy.Soc.A, 223(1152), 1-14, (1954) -% Doetschman/McCool, Chem.Phys. 8, 1-16 (1975) -% Scullane, J.Magn.Reson. 47, 383-397 (1982) -% Jain/Lehmann, phys.stat.sol.(b) 159, 495-544 (1990) - -% These terms are no longer supported. -if isfield(Sys,'aF') - error('Sys.aF is no longer supported. Use Sys.B4 and Sys.B4Frame instead.'); -end - -% High-order terms in extended Stevens operator format -%----------------------------------------------------------------------------- -% Spin system fields: -% B1, B2, B3, ... (k = 1...12) -> processed to Sys.B -% B1Frame, B2Frame, B3Frame, ... -> processed to Sys.BFrame - -% Run over all ranks k -for k = 1:numel(Sys.B) - Bk = Sys.B{k}; - if isempty(Bk), continue; end - - % Run over all desired electron spins - for iSpin = idxElectrons - - % If D is used, skip rank-2 Stevens operator terms - D_present = isfield(Sys,'D') && ~isempty(D) && any(D(:)~=0); - if D_present && k==2, continue; end - - % Skip if all rank-k coefficients are zero - if all(Bk(iSpin,:)==0), continue; end - - % Apply transformation if non-zero tilt angles are given - % (Sys.BFrame is processed from Sys.B?Frame by validatespinsys) - tiltang = Sys.BFrame{k}(iSpin,:); - if any(tiltang) - % Calculate transformation matrix for ISTOs - Dk = wignerd(k,tiltang(1),tiltang(2),tiltang(3)); - Ck = isto2stev(k); - % Transform to Stevens operator basis - DB = Ck*Dk.'/Ck; - DB = removenumericalnoise(DB); - % Transform Bk - Bk_M = Bk(iSpin,:)*DB; % eigenframe of Bk -> molecular frame - else - Bk_M = Bk(iSpin,:); - end - - % Build Hamiltonian - q = k:-1:-k; - for iq = find(Bk_M~=0) - H = H + Bk_M(iq)*stev(Spins,[k,q(iq),iSpin],'sparse'); - end - - end % for all electron spins specified - -end % for all tensor ranks - -H = (H+H')/2; % Hermitianize -if ~sparseResult - H = full(H); -end - -end -%=============================================================================== - -function A_ = removenumericalnoise(A) -reA = real(A); -imA = imag(A); -thr = 1e-14; -reA(abs(reA) g frame + else + R_g2M=eye(3); + end + R_g2M = R_M2g.'; % g frame -> mol frame + gx{eSp} = Sys.g(eSp,1)*R_g2M(:,1)*R_g2M(:,1).'; + gy{eSp} = Sys.g(eSp,2)*R_g2M(:,2)*R_g2M(:,2).'; + gz{eSp} = Sys.g(eSp,3)*R_g2M(:,3)*R_g2M(:,3).'; +end + + +% Calculate magnetic moment from derivatives +muMx_2 = gx{1}(1,1)*dmuMx{1}{1} + gx{1}(1,2)*dmuMx{1}{2} + gx{1}(1,3)*dmuMx{1}{3} +... + gx{2}(1,1)*dmuMx{2}{1} + gx{2}(1,2)*dmuMx{2}{2} + gx{2}(1,3)*dmuMx{2}{3} +... + gy{1}(1,1)*dmuMy{1}{1} + gy{1}(1,2)*dmuMy{1}{2} + gy{1}(1,3)*dmuMy{1}{3} +... + gy{2}(1,1)*dmuMy{2}{1} + gy{2}(1,2)*dmuMy{2}{2} + gy{2}(1,3)*dmuMy{2}{3} +... + gz{1}(1,1)*dmuMz{1}{1} + gz{1}(1,2)*dmuMz{1}{2} + gz{1}(1,3)*dmuMz{1}{3} +... + gz{2}(1,1)*dmuMz{2}{1} + gz{2}(1,2)*dmuMz{2}{2} + gz{2}(1,3)*dmuMz{2}{3}; + +ok(2) = areequal(muMx_2,muMx,1e-10,'abs'); + +muMy_2 = gx{1}(2,1)*dmuMx{1}{1} + gx{1}(2,2)*dmuMx{1}{2} + gx{1}(2,3)*dmuMx{1}{3} +... + gx{2}(2,1)*dmuMx{2}{1} + gx{2}(2,2)*dmuMx{2}{2} + gx{2}(2,3)*dmuMx{2}{3} +... + gy{1}(2,1)*dmuMy{1}{1} + gy{1}(2,2)*dmuMy{1}{2} + gy{1}(2,3)*dmuMy{1}{3} +... + gy{2}(2,1)*dmuMy{2}{1} + gy{2}(2,2)*dmuMy{2}{2} + gy{2}(2,3)*dmuMy{2}{3} +... + gz{1}(2,1)*dmuMz{1}{1} + gz{1}(2,2)*dmuMz{1}{2} + gz{1}(2,3)*dmuMz{1}{3} +... + gz{2}(2,1)*dmuMz{2}{1} + gz{2}(2,2)*dmuMz{2}{2} + gz{2}(2,3)*dmuMz{2}{3}; + +ok(3) = areequal(muMy_2,muMy,1e-10,'abs'); + +muMz_2 = gx{1}(3,1)*dmuMx{1}{1} + gx{1}(3,2)*dmuMx{1}{2} + gx{1}(3,3)*dmuMx{1}{3} +... + gx{2}(3,1)*dmuMx{2}{1} + gx{2}(3,2)*dmuMx{2}{2} + gx{2}(3,3)*dmuMx{2}{3} +... + gy{1}(3,1)*dmuMy{1}{1} + gy{1}(3,2)*dmuMy{1}{2} + gy{1}(3,3)*dmuMy{1}{3} +... + gy{2}(3,1)*dmuMy{2}{1} + gy{2}(3,2)*dmuMy{2}{2} + gy{2}(3,3)*dmuMy{2}{3} +... + gz{1}(3,1)*dmuMz{1}{1} + gz{1}(3,2)*dmuMz{1}{2} + gz{1}(3,3)*dmuMz{1}{3} +... + gz{2}(3,1)*dmuMz{2}{1} + gz{2}(3,2)*dmuMz{2}{2} + gz{2}(3,3)*dmuMz{2}{3}; + +ok(4) = areequal(muMz_2,muMz,1e-10,'abs'); + + +ok=all(ok); \ No newline at end of file diff --git a/tests/ham_hf_deriv_value.m b/tests/ham_hf_deriv_value.m index d60c99e9..beb7cdde 100644 --- a/tests/ham_hf_deriv_value.m +++ b/tests/ham_hf_deriv_value.m @@ -1,4 +1,4 @@ -function ok = test() +function ok = ham_hf_deriv_value() % Test whether hyperfine hamiltonian derivatives are correct, % for a system with one electron and one nucleus @@ -12,5 +12,23 @@ % Calculate hamiltonian from derivatives Hhf_2 = Sys.A(1)*dHhf{1}{1} + Sys.A(2)*dHhf{1}{2} + Sys.A(3)*dHhf{1}{3}; -ok = areequal(Hhf_2,Hhf,1e-10,'abs'); +ok(1) = areequal(Hhf_2,Hhf,1e-10,'abs'); +% test for a larger spin system and different euler angles +Sys.S = [1/2 1]; +Sys.Nucs = '14N,1H'; +Sys.A = [1.435 2.765 3.9876, 1 2 3 ; 4 5 6, 7 8 9]; +Sys.AFrame = [35 26 0, 87 63 42; 81 23 60, 42 36 47]*pi/180; +Sys.J=1; + +[Hhf,dHhf] = ham_hf(Sys); + +% Calculate hamiltonian from derivatives +Hhf_2 = Sys.A(1,1)*dHhf{1,1}{1} + Sys.A(1,2)*dHhf{1,1}{2} + Sys.A(1,3)*dHhf{1,1}{3}+... + Sys.A(2,1)*dHhf{1,2}{1} + Sys.A(2,2)*dHhf{1,2}{2} + Sys.A(2,3)*dHhf{1,2}{3}+... + Sys.A(1,4)*dHhf{2,1}{1} + Sys.A(1,5)*dHhf{2,1}{2} + Sys.A(1,6)*dHhf{2,1}{3}+... + Sys.A(2,4)*dHhf{2,2}{1} + Sys.A(2,5)*dHhf{2,2}{2} + Sys.A(2,6)*dHhf{2,2}{3}; + +ok(2) = areequal(Hhf_2,Hhf,1e-10,'abs'); + +ok=all(ok); \ No newline at end of file diff --git a/tests/ham_zf_deriv_size.m b/tests/ham_zf_deriv_size.m new file mode 100644 index 00000000..9f2b0224 --- /dev/null +++ b/tests/ham_zf_deriv_size.m @@ -0,0 +1,20 @@ +function ok = ham_zf_deriv_size() + +Sys(1).S = [1]; +Sys(1).D = 100*[-1 -1 2]; + +Sys(2).S = [3/2 5/2]; +D = [1 2]'; +E = [0.5 0.3]'; +Sys(2).D = D*[-1,-1,2]/3 + E*[1,-1,0]; +Sys(2).DFrame = [23 96 30 ; 10 6 81]*pi/180; +Sys(2).J=0; + +for k = 1:numel(Sys) + [~,dHzf] = ham_zf(Sys(k)); + elSpins = numel(Sys(k).S); + ok(k,1) = iscell(dHzf); + ok(k,2) = all(size(dHzf)==[1 elSpins]); + nElements = cellfun(@(x)numel(x),dHzf); + ok(k,3) = all(sum(nElements(:)) == elSpins*3); +end diff --git a/tests/ham_zf_deriv_value.m b/tests/ham_zf_deriv_value.m new file mode 100644 index 00000000..68cbdacc --- /dev/null +++ b/tests/ham_zf_deriv_value.m @@ -0,0 +1,32 @@ +function ok = ham_zf_deriv_value() + +% Test whether ZFS hamiltonian derivatives are correct, +% for a system with 1 electron S = 1 + +Sys.S = 1; +Sys.D = 100*[-1 -1 2]; + +[Hzf,dHzf] = ham_zf(Sys); + +% Calculate hamiltonian from derivatives +Hzf_2 = Sys.D(1)*dHzf{1}{1} + Sys.D(2)*dHzf{1}{2} + Sys.D(3)*dHzf{1}{3}; + +ok(1) = areequal(Hzf_2,Hzf,1e-10,'abs'); + + +% test for a larger spin system and asymmetry and euler angles +clear Sys +Sys.S = [3/2 5/2]; +D = [1 2]'; +E = [0.5 0.3]'; +Sys.D = D*[-1,-1,2]/3 + E*[1,-1,0]; +Sys.DFrame = [23 96 30 ; 10 6 81]*pi/180; +Sys.J=0; + +[Hzf,dHzf] = ham_zf(Sys); + +% Calculate hamiltonian from derivatives +Hzf_2 = Sys.D(1,1)*dHzf{1}{1} + Sys.D(1,2)*dHzf{1}{2} + Sys.D(1,3)*dHzf{1}{3}+... + Sys.D(2,1)*dHzf{2}{1} + Sys.D(2,2)*dHzf{2}{2} + Sys.D(2,3)*dHzf{2}{3}; + +ok(2) = areequal(Hzf_2,Hzf,1e-10,'abs');