Dear Stefan,
Taking you up on your advice, I looked at the examples.
The closest one I could find was the 2pEcho simulation using the density matrix. I tried looping through various field values of a field-swept experiment, but I guess this is not what you suggest. Since the relevant part of the echo simulation does not consider the MW frequency, resonance conditions are not calculated even if the system and experiment are otherwise defined (the script at the end of the experiment).
Since saffron exports results in time (ESEEM) or frequency (ENDOR) domains, how would one go about calculating echo intensities of a simple 2pEcho experiment? "Trick" saffron during an ESEEM simulation by npoints = 1 , dt = 0 and tau = the field-sweep experiment tau? Would that be legitimate?
Thanks in advance!
Thanasis
Code: Select all
% Simulation of a two-pulse echo using the density matrix
%==========================================================================
% Here we show how one can simulate a 2-pulse echo transient
% using the density matrix formalism. This is an example of
% advanced usage.
clear all; close all;
cm=100*clight/1e6; % Conversion constant from cm-1 to MHz
%% Define parameters of pulse experiment
n = 300; % signal length
dt = 0.005; % step time [�s]
FWHM = 20; % line width [MHz]
tau = 0.7; % interpulse distance [�s]
% We compute a set of 2-pulse-echo signals, differing in the length of the pulses used.
% The static Hamiltonian is not neglected during the pulse.
offset = 1.2*FWHM*linspace(-1,1,200);
weights = gaussian(offset,0,FWHM);
% Variable tp or fixed?
tp = 0.020;
% tp = 0:0.020:0.180; % pulse length [�s]
%% Define system
gx = 2.002;
gy = 2.002;
gz = 2.002;
gx1 = 2.002;
gy1 = 2.002;
gz1 = 2.002;
gx2 = 2.006;
gy2 = 2.006;
gz2 = 2.006;
r12 = 10; % In Angstrom
R12 = [0; 0; 1]; R12 = R12/norm(R12); % Directional vector between spins Mi-Mj is defined along z
J = -0.01; % In cm-1, +JSiSj formalism
% A1 = [5 5 20]; % Hyperfine S1-A1
% A2 = [10 10 80]; % Hyperfine S2-A2
A1 = 0*[5 5 5]; % Hyperfine S1-A1
A2 = 0*[40 40 40]; % Hyperfine S2-A2
A12 = [0 0 0]; % Hyperfine S1-A2
A21 = [0 0 0]; % Hyperfine S2-A1
% Spin system
Sys1.S= [1/2 1/2];
% Sys1.g=[gx gy gz; gx gy gz];
Sys1.g=[gx1 gy1 gz1; gx2 gy2 gz2];
Sys1.lwpp=[1 0];
Sys1.gFrame = [0 90 -90; 0 90 -90]*pi/180;
if sum(A1) + sum(A2) ~= 0
% A-frames
Sys1.A = [A1 A12; A21 A2];
Sys1.Nucs = '14N,14N';
% Sys1.AFrame = [Aframe1 0 0 0; 0 0 0 Aframe2;]*pi/180;
end
% -------------Describe local g-tensors in molecular frame-------------
% Euler angles for T->M (reverse-opposite of gFrame)
euler1 = -fliplr(Sys1.gFrame(1,:));
euler2 = -fliplr(Sys1.gFrame(2,:));
% Rotation matrices for T->M
R_T2M_1 = erot(euler1);
R_T2M_2 = erot(euler2);
% Express g-vectors as tensors in the T-frame
g1 = diag(Sys1.g(1,:));
g2 = diag(Sys1.g(2,:));
% Transform g-tensors from T to M-frame
g1M = R_T2M_1 * g1 * R_T2M_1.';
g2M = R_T2M_2 * g2 * R_T2M_2.';
% Dipolar exchange with Sys.ee and WITH Easyspin frames in the molecular frame
D12 = transpose(g1M)*g2M - 3 * (transpose(g1M)*R12) * (transpose(R12)*g2M);
dip = -12993 * r12^-3 * D12; % MHz in the +JSiSj formalism
jiso = J*eye(3)*cm;
Sys1.ee = dip + jiso;
%% Define FS experimental conditions
% Define CW experiment and plot a CW spectrum
Exp.mwFreq=9.7; Exp.Range=[340 355]; Exp.nPoints=50; Exp.Harmonic = 1;
Opt.Threshold=0; Opt.Transitions = 'all';
[B,spc_new_eeframes] = pepper(Sys1,Exp,Opt);
figure(1)
plot(B,spc_new_eeframes,'Color','k','DisplayName',"eeD g-Molecular ee frame")
axis tight
xlim([330 360])
% Calculate pulse at each field
for jj = 1:length(B)
% For the spin expectation values of nth state define ... (nth = 1 for the ground state)
nth = 1;
Heig = B(jj);
% Calculate probability amplitudes
H = sham(Sys1,[0,0,Heig]); % Generate the Hamiltonian matrix of the system
[Vt,Ener]=eig(H); % Generate the right-hand vectors (Vt = |>) and the energy eigenvalues (Ener) of the Hamiltonian (in MHz)
[S1x,S1y,S1z,S2x,S2y,S2z] = sop(Sys1,'x1','y1','z1','x2','y2','z2');
[S1p,S1m,S2p,S2m] = sop(Sys1,'+1','-1','+2','-2');
[S1pS2m,S1mS2p,S1zS2z] = sop(Sys1,'+1-2','-1+2','z1z2');
S1 = S1x + S1y + S1z;
S2 = S2x + S2y + S2z;
S1S2 = 1/2 * (S1pS2m + S1mS2p) + S1zS2z; % S1*S2
S_2 = S1*S1 + S2*S2 + 2*S1S2; %ST^2
Sx = S1x+S2x;
Sy = S1y+S2y;
Sz = S1z+S2z;
flipAngle = pi/2;
signal = zeros(n,length(tp));
for p = 1:length(tp)
for k = 1:length(offset)
Ham0 = offset(k)*Sz;
Pulse = expm(-1i*(flipAngle*Sx+2*pi*tp(p)*Ham0));
TauEvol = expm(-2i*pi*tau*Ham0);
U = Pulse^2*TauEvol*Pulse;
signal(:,p) = signal(:,p) + ...
weights(k)*real(evolve(U*Sz*U',Sy,Ham0,n,dt));
end
end
% Result: The longer the pulses, the broader the
% echo.
figure(2)
plot((0:n-1)*dt,signal);
hold on
title('The shape of the primary echo depending on the pulse length');
xlabel('time after the pi pulse [μs]');
ylabel('echo signal');
end