Overhaul of g/A/D/E strain broadening

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

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.