diff --git a/easyspin/ham_ee.m b/easyspin/ham_ee.m index 89b04276..4427a8e5 100644 --- a/easyspin/ham_ee.m +++ b/easyspin/ham_ee.m @@ -1,11 +1,11 @@ -% ham_ee Electron-electron spin interaction Hamiltonian +% ham_ee Electron-electron spin interaction Hamiltonian % -% Hee = ham_ee(SpinSystem) -% Hee = ham_ee(SpinSystem,eSpins) -% Hee = 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,14 +18,16 @@ % Output: % - Hee: 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 Hee = 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 @@ -50,7 +52,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); @@ -77,11 +79,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),:); @@ -92,30 +96,57 @@ 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'); + 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'); - Hee = Hee + so1*J(c1,c2)*so2; + 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; + 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 Hee2 = 0; for c = 1:3 Hee2 = Hee2 + sop(sys,[Pairs(iPair,1),c;Pairs(iPair,2),c],'sparse'); end - Hee = Hee + System.ee2(iCoupling)*Hee2^2; -end - -Hee = (Hee+Hee')/2; % Hermitianise -if ~useSparseMatrices - Hee = full(Hee); % sparse -> full + 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 end diff --git a/easyspin/ham_ez.m b/easyspin/ham_ez.m index 20c02be7..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~=3, 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,49 +97,69 @@ % 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 g = R_g2M*g*R_g2M.'; + else + R_g2M=eye(3); + end + + % 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 = 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; + + 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); - end - varargout = {muxM, muyM, muzM}; + % 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_hf.m b/easyspin/ham_hf.m index 2f641eed..e877a6d8 100644 --- a/easyspin/ham_hf.m +++ b/easyspin/ham_hf.m @@ -1,18 +1,18 @@ -% 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. % % 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,21 +102,49 @@ 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'); + 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 -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_zf.m b/easyspin/ham_zf.m index f26f31e4..640d07cd 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 % % Hzf = ham_zf(SpinSystem) % Hzf = 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 [Hzf,dHzf] = ham_zf(SpinSystem,idxElectrons,opt) if nargin==0, help(mfilename); return; end @@ -25,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; @@ -35,13 +37,15 @@ error('Electron spin index/indices (2nd argument) out of range!'); end -H = sparse(Sys.nStates,Sys.nStates); +nStates=Sys.nStates; +Hzf = 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,26 +57,47 @@ 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) + 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 + 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'); + 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 - H = H + D(c1,c2)*(Sxyz{c1}*Sxyz{c2}); + tempProduct=Sxyz{c1}*Sxyz{c2}; % storing the matrix product temporarily + 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 + dHzf{iSpin} = {dHzfdDx,dHzfdDy,dHzfdDz}; end % Fourth-order terms a and F @@ -110,7 +135,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,22 +151,18 @@ 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'); + 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/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..8e4aa357 --- /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, 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 D frame + if any(Sys.DFrame(iEl,:)) + 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; + 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_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 new file mode 100644 index 00000000..7765ad19 --- /dev/null +++ b/easyspin/private/strains_ga.m @@ -0,0 +1,85 @@ +% 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 + %%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.'); + 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/resfields.m b/easyspin/resfields.m index a1c16477..845607cf 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 @@ -718,90 +715,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 @@ -924,7 +861,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 @@ -1282,31 +1219,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 @@ -1320,18 +1256,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 3c3b8592..e5cfc296 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(:)) @@ -604,12 +608,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 0ceca3d8..977540ee 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. @@ -532,7 +532,8 @@ % D strain %----------------------------------------------- - [useDStrain,dHdD,dHdE] = getdstrainops(CoreSys); + sigmadHdDE = strains_de(CoreSys); + useDStrain = ~isempty(sigmadHdDE); % g-A strain %------------------------------------------------- @@ -545,7 +546,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 @@ -558,7 +560,7 @@ end else for iEl = CoreSys.nElectrons:-1:1 - gStrainMatrix{iEl} = 0; + gStrainMatrix{iEl} = zeros(3); end end @@ -642,9 +644,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); @@ -657,18 +659,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 @@ -775,32 +777,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 @@ -829,7 +833,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/ham_ee_deriv_size.m b/tests/ham_ee_deriv_size.m new file mode 100644 index 00000000..44488414 --- /dev/null +++ b/tests/ham_ee_deriv_size.m @@ -0,0 +1,20 @@ +function ok = ham_ee_deriv_size() + +% simple spin system +Sys(1).S = [1/2 1/2]; +Sys(1).ee = [1 2 3]; +Sys(1).ee2 = 100; + +% more complex one +Sys(2).S = [1/2 1 1]; +Sys(2).ee = [1 2 3; 4 5 6 ; 7 8 9]; +Sys(2).ee2 = [100 200 300]; + +for k = 1:numel(Sys) + [~,dHee] = ham_ee(Sys(k)); + elSpins = numel(Sys(k).S); + ok(k,1) = iscell(dHee); %are these cell? + ok(k,2) = all(size(dHee)==[1 elSpins*(elSpins-1)/2]); % what about their size + nElements = cellfun(@(x)numel(x),dHee); + ok(k,3) = all(sum(nElements(:)) == elSpins*(elSpins-1)/2*4); %ensure the size are correct +end diff --git a/tests/ham_ee_deriv_value.m b/tests/ham_ee_deriv_value.m new file mode 100644 index 00000000..8ccc729a --- /dev/null +++ b/tests/ham_ee_deriv_value.m @@ -0,0 +1,34 @@ +function ok = ham_ee_deriv_value() + +% Test whether electron-electron hamiltonian derivatives are correct, +% for a system with 2 electrons + +Sys.S = [1/2 1/2]; +J=10; +Sys.ee = [1 1 -2]+ J; + +[Hee,dHee] = ham_ee(Sys); + +% Calculate hamiltonian from derivatives +Hee_2 = Sys.ee(1)*dHee{1}{1} + Sys.ee(2)*dHee{1}{2} + Sys.ee(3)*dHee{1}{3}; + +ok(1) = areequal(Hee_2,Hee,1e-10,'abs'); + +clear Sys +% test for a larger spin system and bilinear term +Sys.S = [1/2 1 1/2]; +dip=[[1 1 -2];2*[1 1 -2];3*[1 1 -2]]; +J=[1 2 3]; +Sys.ee=dip+J'; +Sys.eeFrame=[10 23 65; 78 42 136; -52 86 276]*pi/180; +Sys.ee2=[100 200 300]; +[Hee,dHee] = ham_ee(Sys); + +% Calculate hamiltonian from derivatives +Hee_2 = Sys.ee(1,1)*dHee{1}{1} + Sys.ee(1,2)*dHee{1}{2} + Sys.ee(1,3)*dHee{1}{3}+... + Sys.ee(2,1)*dHee{2}{1} + Sys.ee(2,2)*dHee{2}{2} + Sys.ee(2,3)*dHee{2}{3}+... + Sys.ee(3,1)*dHee{3}{1} + Sys.ee(3,2)*dHee{3}{2} + Sys.ee(3,3)*dHee{3}{3}; + +Hee_2 = Hee_2 + Sys.ee2(1)*dHee{1}{4} + Sys.ee2(2)*dHee{2}{4} + Sys.ee2(3)*dHee{3}{4}; + +ok(2) = areequal(Hee_2,Hee,1e-10,'abs'); diff --git a/tests/ham_ez_deriv_size.m b/tests/ham_ez_deriv_size.m new file mode 100644 index 00000000..602a3cd3 --- /dev/null +++ b/tests/ham_ez_deriv_size.m @@ -0,0 +1,28 @@ +function ok = ham_ez_deriv_size() + +% check the output with a simple case +Sys(1).S = 1/2; +Sys(1).g = 2*[1 1 1]; + +% check the output with a two electrons and a relative orientation between g tensors +Sys(2).S = [1/2 1/2]; +Sys(2).g = [2*[1 1.9 1.1]; 3*[1 1 2]]; +Sys(2).gFrame = [23 96 30 ; 10 6 81]*pi/180; +Sys(2).J=0; + +for k = 1:numel(Sys) + [muMx,muMy,muMz,dmuMx,dmuMy,dmuMz] = ham_ez(Sys(k)); + elSpins = numel(Sys(k).S); + ok(k,1) = iscell(dmuMx); + ok(k,2) = iscell(dmuMy); + ok(k,3) = iscell(dmuMz); + ok(k,4) = all(size(dmuMx)==[1 elSpins]); + ok(k,5) = all(size(dmuMy)==[1 elSpins]); + ok(k,6) = all(size(dmuMz)==[1 elSpins]); + nElements = cellfun(@(x)numel(x),dmuMx); + ok(k,7) = all(sum(nElements(:)) == elSpins*3); + nElements = cellfun(@(x)numel(x),dmuMy); + ok(k,8) = all(sum(nElements(:)) == elSpins*3); + nElements = cellfun(@(x)numel(x),dmuMz); + ok(k,9) = all(sum(nElements(:)) == elSpins*3); +end diff --git a/tests/ham_ez_deriv_value.m b/tests/ham_ez_deriv_value.m new file mode 100644 index 00000000..461cd90a --- /dev/null +++ b/tests/ham_ez_deriv_value.m @@ -0,0 +1,69 @@ +function ok = ham_ez_deriv_value() + +% Test whether magnetic moment and derivatives are correct, +% for a system with one electron + +Sys.S = 1/2; +Sys.g = 2*[1 1 1]; + +[muMx,muMy,muMz,dmuMx,dmuMy,dmuMz] = ham_ez(Sys); + + +% Calculate hamiltonian from derivatives +mu = muMx + muMy + muMz; +mu_2 = Sys.g(1)*dmuMx{1}{1} + Sys.g(2)*dmuMy{1}{2} + Sys.g(3)*dmuMz{1}{3}; + +ok(1) = areequal(mu_2,mu,1e-10,'abs'); + +% test for a larger spin system and different euler angles +Sys.S = [1/2 1/2]; +Sys.g = [2*[1 1.9 1.1]; 3*[1 1 2]]; +Sys.gFrame = [35 26 0; 81 23 60]*pi/180; +Sys.J=1; + +[muMx,muMy,muMz,dmuMx,dmuMy,dmuMz] = ham_ez(Sys); + +for eSp=1:numel(Sys.S) + ang = Sys.gFrame(eSp,:); + if any(ang) + R_M2g = erot(ang); % mol frame -> 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_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..beb7cdde --- /dev/null +++ b/tests/ham_hf_deriv_value.m @@ -0,0 +1,34 @@ +function ok = ham_hf_deriv_value() + +% 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(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'); 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 -