Overhaul of g/A/D/E strain broadening

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

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
 

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