Simulate also hfi with three 1H

reports/oop_ciss_trEPR_stickSpectra.tex

@@ -24,26 +24,30 @@
 		\item g1 = [2.0030, 2.0026, 2.0023], g2 = [2.0062, 2.0051, 2.0022] for P700+ and A1- respectively;
 		\item Euler angles [-10, -128, -83] to transform from the reference frame of P700+ to the one of A1-;
 		\item dipolar interaction and exchange coupling dip = -0.177 mT and J = 0.001 mT (following the spin Hamiltonian convention used by Zech);
-		\item direction of the spin-spin interaction corresponding to the negative x-axis of the A1- frame of reference (hence Euler angles [0, 90, 0] to go from the dipolar frame of reference to the A1- frame of reference).
+		\item $z_D$, direction of the spin-spin interaction, corresponding to the negative x-axis of the A1- frame of reference (hence Euler angles [0, 90, 0] to go from the dipolar frame of reference to the A1- frame of reference);
+		\item hyperfine coupling of the A1- radical with three equivalent hydrogen nuclei, A = [9, 9, 12.8] MHz and [-60, 90, 0] Euler angles to transform from the frame of the hyperfine to the one of A1-.
 	\end{itemize}
 
 	\section{Transition frequencies and intensities}
-	To calculate the powder average, we keep our frame fixed with the frame of A1- and calculate the trEPR stick spectra for different orientations of the external magnetic field.\\
-	As explained in Zech's thesis \cite{zechPulsedTransientElectron1998} chapter 4, the stick spectrum at a given orientation of the magnetic field is going to be similar to what shown in Fig. \ref{fig:stickSpectrum_zechThesis}. The position of the resonances is given by the values of J, dip and:
+	To calculate the powder average, we keep our frame fixed with the frame of A1- and calculate the trEPR stick spectra for different orientations of the external magnetic field. As explained in Zech's thesis \cite{zechPulsedTransientElectron1998} chapter 4, the stick spectrum at a given orientation of the magnetic field is going to be similar to what shown in Fig. \ref{fig:stickSpectrum_zechThesis}.	As shown in Fig. \ref{fig:stickSpectrum_zechThesis}, the position of the resonances is given by the values of J, dip and:
 	\begin{equation}
-		\omega_0 = \frac{ \mu_B B_0}{h}\frac{g1 - g2}{2}
+		\omega_0 = \frac{ \mu_B B_0}{h}(g1 + g2)/2 + \frac{1}{2}\sum_{j} (A_{1, j} + A_{2, j})m_j,
+		\label{eq:formula_for_omega_0}
 	\end{equation}
 	and:
 	\begin{equation}
 		\Omega = \sqrt{\Delta\omega^2 + (J + d/2)^2},
+		\label{eq:formula_for_Omega}
 	\end{equation}
 	where:
 	\begin{equation}
-	\Delta\omega = \frac{ \mu_B B_0}{h}\frac{g1 - g2}{2}.
+	\Delta\omega = \frac{ \mu_B B_0}{h}(g1 - g2)/2 + \frac{1}{2}\sum_{j} (A_{1, j} - A_{2, j})m_j.
+	\label{eq:formula_for_Delta_omega}
 	\end{equation}
 	The absolute value of the intensity is the same for all the transitions at a given orientation. In particular:
 	\begin{equation}
 	I_{12} = -I_{34} = I_{13} = -I_{24} \propto \frac{\Delta\omega^2}{\Omega^2} = \sin^2(2\alpha),
+	\label{eq:formula_for_intensity_resonances}
 	\end{equation}
 	where $\alpha$ is the mixing angle used in Zech's thesis, defined as:
 	\begin{equation}
@@ -63,10 +67,11 @@
 
 	\section{Computational steps}
 	\begin{enumerate}
-		\item Calculate $d$ and effective g-values for each orientation in order to calculate $\omega_0$, $\Delta\omega$ and $\Omega$
-		\item Find the transition frequencies and intensities
-		\item Apply a line broadening to the stick spectrum
-		\item Average over all orientations
+		\item Calculate $d$, effective g-values and effective hyperfine coupling for each orientation
+		\item Calculate $\omega_0$, $\Delta\omega$ and $\Omega$
+		\item Find the transition frequencies
+		\item Calculate the lineshapes for each transition with a certain linewidth and the correct intensity
+		\item Sum all the lines and average over all orientations
 	\end{enumerate}
 
 	\subsection{Calculate orientation-dependent parameters}
@@ -79,7 +84,7 @@ nTheta = numel(thetas);
 phis = (0:3:360)*pi/180;
 nPhi = numel(phis);
 		\end{lstlisting}
-	In this case the array thetas will be thetas = [0, 3, 6, ..., 180]*pi/180 and the number of angles will be nTheta = 61.\\
+	In this case the array thetas will be thetas = [0, 3, 6, ..., 180]$^o$ and the number of angles will be nTheta = 61.\\
 	From here we can calculate the effective g-values at each orientations. First we calculate the versor $n$ indicating the direction of $B_0$:
 		\begin{lstlisting}[language=Octave]
 % Direction of B0
@@ -92,6 +97,7 @@ nVers(3, :, :) = cos(thetas')*ones(1, nPhi);
 	Next step: calculate g2, that is the effective g-value of A1-. Since the g2 tensor is diagonal in this reference frame, for each orientation it will be:
 	\begin{equation}
 		g\textsubscript{eff} = \sqrt{(g_{xx}\cdot n_x)^2 + (g_{yy}\cdot n_y)^2 + (g_{zz}\cdot n_z)^2}.
+		\label{eq:formula_for_g_effective}
 	\end{equation}
 	We basically want to multiply element-wise g2 by the value of nVers for each orientation and then take the square root of the sum of the squares of the projections. The code looks like:
 			\begin{lstlisting}[language=Octave]
@@ -101,9 +107,173 @@ g2 = squeeze(sqrt( sum( (Sys.g(2, :)'.*nVers).^2)));
 	First, the element-wise multiplication is carried out: Sys.g(2, :)' is the column vector 3 x 1 of [$g_{xx}, g_{yy}, g_{zz}$]. The operator .* does the element-wise multiplication at each fixed $\theta$ and $\phi$ value. As a consequence, the result of Sys.g(2, :)'.*nVers is still a 3 x nTheta x nPhi matrix. Afterwards these values are squared (also element wise), then they are summed along the first dimension and then the square root is calculated. At this point the result is a 1 x nTheta x nPhi matrix. The function squeeze() makes it a nTheta x nPhi matrix corresponding to the effective g-values for each orientation.\\
 	For the P700+ radical the process to calculate the effective g-values is the same, but we first need to transform the g-tensor to the frame of A1-. The code is:
 			\begin{lstlisting}[language=Octave]	
-g1TensorInFrame2 = erot(eAngles)'*diag(Sys.g(1, :))*erot(eAngles);
+euMatrixg1 = erot(Sys.gFrame(1, :));
+g1TensorInFrame2 = euMatrixg1'*diag(Sys.g(1, :))*euMatrixg1;
 g1 = squeeze(sqrt( sum( (pagemtimes(g1TensorInFrame2, nVers)).^2)));
 			\end{lstlisting}
-	Here erot(eAngles) is the rotation matrix with Euler angles [-10, -128, -83]. It is applied transposed on the left (' sign) and not transposed on the right side of the diagonal matrix (diag(Sys.g(1, :))). Afterwards we calculate the effective g-value multiplying the g-tensor by the nVers for each orientation. This is done by the function pagemtimes(), that enables multiplications of 3 x 3 matrices by 3 x 1 vectors for each value of $\theta$ and $\phi$. The rest is analogous to the A1- case.\\
+	Here erot(Sys.gFrame(1, :)) is the rotation matrix with Euler angles Sys.gFrame(1, :) = [-10, -128, -83]. It is applied transposed on the left (' sign) and not transposed on the right side of the diagonal matrix (diag(Sys.g(1, :))). Afterwards we calculate the effective g-value multiplying the g-tensor by the nVers for each orientation. This is done by the function pagemtimes(), that enables multiplications of 3 x 3 matrices by 3 x 1 vectors for each value of $\theta$ and $\phi$. The rest is analogous to the A1- case.\\
+	For the dipolar interaction:
+			\begin{lstlisting}[language=Octave]	
+% Dipolar interaction
+zD = erot(Sys.eeFrame)*[0, 0, 1]';
+dd = squeeze(dipFunc(dip, nVers, zD));
+			\end{lstlisting}
+	We first transform the $z_D$ vector to the A1- frame using the Euler angles Sys.eeFrame = [0, 90, 0], therefore obtaining zD = [-1, 0, 0]. Afterwards we the function dipFunc to calculate the value of the dipolar coupling for every orientation. The function is defined as:
+			\begin{lstlisting}[language=Octave]	
+% Dipolar interaction
+% Makes use of the fact that cos(thetaD) = dotProduct(B0, zD)
+% Expected input:
+%   dip:    1 x 1
+%   nVers:  3 x nTheta x nPhi
+%   zD:     3 x 1
+dipFunc = @(dip, nVers, zD) dip*((sum(nVers.*zD)).^2 - 1/3);	
+			\end{lstlisting}
+	Finally we calculate the effective hyperfine interaction for each orientation.
+			\begin{lstlisting}[language=Octave]	
+%
+% Hyperfine interaction with a number nNuc of equal nuclei
+%
+% For spin 2
+euMatrixHfi2 = erot(Sys.AFrame(1, 4:6));
+Ahfi2TensorInFrame2 = euMatrixHfi2'*diag(Sys.A(1, 4:6))*euMatrixHfi2;
+Ahfi2 = squeeze(sqrt( sum( (pagemtimes(Ahfi2TensorInFrame2', nVers)).^2)));
+			\end{lstlisting}
+	We calculate it only for spin 2, being the radical on A1-, because that is the electron that should interact with the nNuc = 3 identical hydrogen atoms. erot(Sys.AFrame(1, 4:6)) is the rotation matrix using Euler angles [-60, 90, 0] and it is used to transform the values of Sys.A(1, 4:6) = [9, 9, 12.8] MHz to the reference frame of A1-. Afterwards, the effective hyperfine Ahfi2 is calculated, analogously to Eq. \ref{eq:formula_for_g_effective}.
+	
+	\subsection{Calculate $\omega_0$, $\Delta\omega$, $\Omega$ and the transition frequencies}
+	We want to calculate $\omega_0$, $\Delta\omega$ and $\Omega$ using Eq. \ref{eq:formula_for_omega_0}, Eq. \ref{eq:formula_for_Delta_omega} and Eq. \ref{eq:formula_for_Omega}. First we calculate $\omega_0$ and $\Delta\omega$ neglecting the hyperfine couplings:
+			\begin{lstlisting}[language=Octave]	
+w0 = g2wFunc(B0, (g1 + g2)/2);  % GHz
+deltaw = g2wFunc(B0, (g1 - g2)/2);  % GHz	
+			\end{lstlisting}	
+	where:
+			\begin{lstlisting}[language=Octave]		
+% Expected input: B0 in mT
+% Output: frequency nu in GHz
+g2wFunc = @(B0, g) bmagn*B0/planck*g*1e-12;	
+			\end{lstlisting}
+	Afterwards we define $A_+ = A_1 + A_2$ (in our case $A_1 = 0$ because there is no interaction between electron 1, that is the radical on P700+, and the nuclei), we adjust the shapes of the vectors and then we add to $\omega_0$ the following values: $[-\text{nNuc}/2]A_+$, $ \cdot [-(\text{nNuc} - 1)/2]A_+$, ... , $[(\text{nNuc} - 1)/2]A_+$, $[\text{nNuc}/2]A_+$. We obtain therefore the values of Eq. \ref{eq:formula_for_omega_0} for each orientation.\\
+	We follow analogous steps to obtain $\Delta\omega$ as in Eq. \ref{eq:formula_for_Delta_omega}.\\
+	The code look like:
+			\begin{lstlisting}[language=Octave]			
+% Shift resonances due to hyperfine interaction
+Aplus = (Ahfi1 + Ahfi2)/2*1e-3;
+% Make w0 have size [nHfiLine, nTheta, nPhi]
+w0 = insertDimensionPos1(w0, nHfiLine);
+% Make Aplus have size [1, nTheta, nPhi]
+Aplus = reshape(Aplus, [1, size(Aplus)]);
+% Multiply by all possible values of sum of quantum number of the nuclei
+% Final size = [nHfiLine, nTheta, nPhi]
+Aplus = pagemtimes((-nNuc/2:1:nNuc/2)', Aplus);
+w0 = w0 + Aplus;
+% Same for deltaw and Aminus
+Aminus = (Ahfi1 - Ahfi2)/2*1e-3;
+deltaw = insertDimensionPos1(deltaw, nHfiLine); 
+Aminus = reshape(Aminus, [1, size(Aminus)]);
+Aminus = pagemtimes((-nNuc/2:1:nNuc/2)', Aminus);
+deltaw = deltaw + Aminus;	
+			\end{lstlisting}
+	In our case for example, we have nNuc = 3, which means that we expect a number of lines equal to $\text{nHfiLine} = 2\cdot\text{nNuc}\cdot 1/2 + 1 = \text{nNuc} + 1 = 4$. The distance of the lines from $\omega_0$ will be: $-3/2\cdot A_+$, $-1/2\cdot A_+$, $1/2\cdot A_+$, $3/2\cdot A_+$. That is what the code is doing.\\
+	We must take into consideration that the intensity of the lines will follow the Pascal triangle. Therefore we calculate the Pascal triangle:
+			\begin{lstlisting}[language=Octave]				
+% Pascal factors because some values of the sum of the quantum number come
+% up more than others
+pascalMatrix = pascal(nHfiLine);
+pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);  % Antidiag
+			\end{lstlisting}
+	With the Matlab built-in function pascal(n) we obtain a square matrix, where each super anti-diagonal corresponds to a row of the Pascal triangle. Since we are interested in the row number nHfiLine, we generate the square pascal(nHfiLine) and then we calculate the antidiagonal and store it in the vector pascalFactor. This will be useful later to properly calculate the intensity of the lines.
+	Finally we calculate $\Omega$:
+				\begin{lstlisting}[language=Octave]				
+Omega = OmegaFunc(deltaw, JJ, insertDimensionPos1(dd, nHfiLine));  % GHz
+				\end{lstlisting}
+	where:
+				\begin{lstlisting}[language=Octave]				
+% Expected input: deltaw in GHz, J in MHz, d in MHz
+% Output: Omega in GHz
+OmegaFunc = @(deltaw, J, d) sqrt(deltaw.^2 + (J*1e-3 + d/2*1e-3).^2);		
+				\end{lstlisting}
+	The position of each resonance depends on the difference between the energies of the eigenstates of the Hamiltonian. Following Zech's thesis \cite{zechPulsedTransientElectron1998}:
+				\begin{lstlisting}[language=Octave]					
+dipInteraction = (JJ - insertDimensionPos1(dd, nHfiLine))*1e-3;  % GHz
+wReson(1, :, :, :) = w0 - dipInteraction - (Omega);  % w12
+wReson(2, :, :, :) = w0 + dipInteraction - (Omega);  % w34
+wReson(3, :, :, :) = w0 - dipInteraction + (Omega);  % w13
+wReson(4, :, :, :) = w0 + dipInteraction + (Omega) ;  % w24
+				\end{lstlisting}
+	The function insertDimensionPos1 is used to have proper dimensions and can be ignored here. The matrix wReson has dimensions [4, nHfiLine, nTheta, nPhi].\\
+
+			
+	\subsection{Calculate the lineshapes for each transition with a certain linewidth and the correct intensity}
+		At this point we calculate the lines for each of these transitions:
+	\begin{lstlisting}[language=Octave]							
+% Lineshapes for each transition
+trlw = trlw*1e-3;  % GHz
+trSignal = ...
+	gaussiantransitions(xxSim', wReson, trlw, "fwhm");	
+	\end{lstlisting}
+	where guassiantransitions returns the gaussian lineshapes in the range xxSim centered at wReson with FWHM equal to trlw (which in our case would be around 12 MHz).\\
+	The intensity of each of these lines is determined by Eq. \ref{eq:formula_for_intensity_resonances}:
+	\begin{lstlisting}[language=Octave]						
+intensityReson = 1/8*(deltaw.^2)./(Omega.^2);
+% Sign of intensityReson alternates for the transitions w12, w34, w13, w24
+intensityReson = [1; -1; 1; -1].*insertDimensionPos1(intensityReson, 4); 	
+	\end{lstlisting}	
+	The normalization factor such that the integral of each of these lines is equal to the intensities of Eq. \ref{eq:formula_for_intensity_resonances} is equal to:
+	\begin{lstlisting}[language=Octave]						
+% Normalization to account for possible different linewidths
+intensityNorm = insertDimensionPos1(intensityReson, 1) ./ sum(trSignal);
+	\end{lstlisting}	
+	Finally we multiply the lineshapes by the normalization factor and the Pascal triangle factors:
+		\begin{lstlisting}[language=Octave]	
+trSignal = trSignal.*intensityNorm.*insertDimensionPos1(pascalFactor, 1);
+		\end{lstlisting}
+	
+	
+	\subsection{Sum all the lines and average over all orientations}
+	The final part of the code is:
+		\begin{lstlisting}[language=Octave]			
+trSignalSumHfi = squeeze(sum(trSignal, 3));  % Sum over hfi lines
+trSignalSum = squeeze(sum(trSignalSumHfi, 2));  % Sum over transitions
+
+solidAngleWeight = sin(thetas)/sum(sin(thetas')*ones(1, nPhi), 'all');
+trSignalPowder = sum(squeeze(sum(trSignalSum.*solidAngleWeight, 3)), 2);	
+		\end{lstlisting}
+	The matrix trSignalSumHfi is a [nAx, 4, nTheta, nPhi] matrix, where the lineshapes due to the hyperfine have been summed. Here nAx is the number of points of the frequency sweep that we are calculating (usually more than 100, I use 1024 or 2048). Each line of trSignalSumHfi is long nAx and is equivalent to the spectrum of the separate transitions $1\rightarrow2$ or $3\rightarrow4$ or $1\rightarrow3$ or $2\rightarrow4$ for each orientation. Note that the Pascal triangle factors were already taken into consideration.\\
+	The matrix trSignalSum is a [nAx, nTheta, nPhi] matrix, where the lineshapes due to the four transitions have been summed. This matrix shows the crystal spectrum for each orientation.\\
+	Finally, trSignalPowder is an array long nAx, where all the contribution have been averaged over all $\theta$ and $\phi$ values with weight equal to $\sin(\theta)$ (normalized by the sum of $\sin(\theta)$ over the all sphere, neglecting global proportionality factors that play no role in determining the lineshape of the spectrum). 
+		
+			
+	\section{Comparison with Easyspin}
+	We calculate the spectrum using our code with parameters:
+	\begin{itemize}
+		\item dip = -0.177 mT = -4.9 MHz and JJ = 0.001 mT = 0.028 MHz;
+		\item linewidth trlw = 0.42 mT = 12 MHz;
+		\item all the other parameters as specified at the beginning of this document.
+	\end{itemize}
+	We compare it with a simulation obtained using Easyspinusing the same parameters (changing the sign of the dipolar coupling and dividing it by 1.5 due to a different definition of the spin Hamiltonian). The comparison is in Fig. \ref{fig:zech_psiOopEseem_trEPRstickSpectra_EasyspinComparison}.
+	\begin{figure}
+		\centering
+		\includegraphics[width=.9\textwidth]{../images/zech_psiOopEseem_trEPRstickSpectra_EasyspinComparison}
+		\caption{Comparison between Easyspin and our spectrum.}
+		\label{fig:zech_psiOopEseem_trEPRstickSpectra_EasyspinComparison}
+	\end{figure}
+	
+	\section{Comparison with experimental data}
+	The calculation from our code is compared to experimental data in Fig. \ref{fig:zech_psiOopEseem_trEPRstickSpectra_ExperimentComparison} (the spectrum is flipped and the axis is transformed to magnetic field axis). The high field part of the spectrum is not well reproduced. In \cite{kamlowskiRadicalPairState} Fig. 5, the spectrum is fitting better the high field part. There they used different linewidths: the one for P700+ equal to 0.7 mT and the one for A1- equal to 0.335 mT. In our code we used the same value of 0.42 mT for both.
+	\begin{figure}
+		\centering
+		\includegraphics[width=.9\textwidth]{../images/zech_psiOopEseem_trEPRstickSpectra_ExperimentComparison}
+		\caption{Comparison between experimental data and our spectrum.}
+		\label{fig:zech_psiOopEseem_trEPRstickSpectra_ExperimentComparison}
+	\end{figure}
+	
+	
+	
+	
+	
+	
 	\printbibliography
 \end{document}
+
+
+

reports/reportOopCissCalc.bib

@@ -45,6 +45,13 @@
   file = {/home/gianlum33/Zotero/storage/R26T3BKN/Zech - 1998 - Pulsed and transient electron paramagnetic resonan.pdf}
 }
 
+@article{kamlowskiRadicalPairState,
+  title = {The {{Radical Pair State P7}}•+{{00A1}}•- in {{Photosystem I Single Crystals}}: {{Orientation Dependence}} of the {{Transient Spin-Polarized EPR Spectra}}},
+  author = {Kamlowski, Andreas and Zech, Stephan G and Fromme, Petra and Bittl, Robert and Lubitz, Wolfgang and Witt, Horst T and Stehlik, Dietmar},
+  langid = {english},
+  file = {/home/gianlum33/Zotero/storage/SPJSNCLM/Kamlowski et al. - The Radical Pair State P7•+00A1•- in Photosystem I.pdf}
+}
+
 @article{timmelSpincorrelatedRadicalPairs1998,
   title = {Spin-Correlated Radical Pairs: Microwave Pulse Effects on Lifetimes, Electron Spin Echo Envelope Modulations, and Optimum Conditions for Detection by Electron Spin Echo Spectroscopy},
   shorttitle = {Spin-Correlated Radical Pairs},

zech_psiOopEseem_powderAverage_inhBroadening.m

@@ -2,7 +2,7 @@
 clearvars
 
 % Import
-% cd('/net/storage/gianlum33/projects/oop_ciss_calculations/')
+cd('/net/storage/gianlum33/projects/oop_ciss_calculations/')
 expName = 'data_extracted/zech_p46_oopEseem_expData.csv';
 fitName = 'data_extracted/zech_p46_oopEseem_fit.csv';
 

zech_psiOopEseem_trEPRstickSpectra.m

@@ -3,6 +3,7 @@ clearvars
 
 %% Easyspin simulation
 
+tic
 % Taken from Easyspin example
 clear('Sys', 'Exp', 'Vary', 'VaryExp')
 
@@ -13,14 +14,21 @@ Sys.gFrame = [-10 -128 -83; ...
                 0    0   0]*pi/180;
 Sys.eeFrame = [0 90 0]*pi/180;  % zD direction is along -x of A1-
 
-Sys.J = unitconvert(-1e-3,'mT->MHz'); % MHz
-Sys.dip = unitconvert(+0.177,'mT->MHz'); % MHz
-Sys.lw = 10;  % MHz
+% Factor 1/1.5 coming from different definition of spin Hamiltonian?
+Sys.J = unitconvert(-1e-3,'mT->MHz')/1.5; % MHz
+Sys.dip = unitconvert(+0.177,'mT->MHz')/1.5; % MHz
+Sys.lw = 12;  % MHz
+% Sys.lw = mhz2mt(11);  % MHz
 Sys.initState = 'singlet';
-% Vary.dip = 1;
-% Vary.J = 0.1;
-% Vary.lw = 5;
-% VaryExp.mwCenterSweep = [0.1, 0];
+
+nNuc = 3;
+Sys.Nucs = '1H';
+Sys.Nucs = repmat(append(Sys.Nucs, ', '), [1, nNuc]);
+Sys.Nucs = Sys.Nucs(1:end - 2);
+Sys.A = [zeros(1, 3), 9, 9, 12.8];
+Sys.A = repmat(Sys.A, [nNuc, 1]);
+Sys.AFrame = [0 0 0 -60 90 0]*pi/180;
+Sys.AFrame = repmat(Sys.AFrame, [nNuc, 1]);
 
 B0 = 345.9;  % mT
 mwFreq = 9.7;
@@ -39,30 +47,49 @@ Exp.nPoints = numel(xxSim1);
 Exp.Harmonic = 0;
 
 % Opt.GridSize = 10;
+% Opt.separate = 'transitions';
 FitOpt.Verbosity = 0;
+FitOpt.maxTime = 0.5;
 
+% Exp.CrystalSymmetry = 1;       % space group Pmmm = #47
+% Exp.MolFrame = [0 0*pi/32 0];      % molecular frame aligned with crystal/sample frame
+% Exp.SampleFrame = [0 0 0]*pi/180; 
+
+% ySimNoHfi = pepper(Sys, Exp);
 ySim = pepper(Sys, Exp);
 
 hold on
 clf
+% plot(xxSim1, ySim/max(ySim), xxSim1 - 0.21, flip(trSignalPowder/max(trSignalPowder)))
 plot(xxSim, ySim/max(ySim), xxSim, trSignalPowder/max(trSignalPowder))
+% plot(xxSim, ySimNoHfi/max(ySimNoHfi), xxSim, ySim/max(ySim))
 legend('Easyspin', 'Gianluca')
 
 Ham = ham(Sys, [0, 0, B0]);
 [Ham0, mux, muy, muz] = ham(Sys);
+toc
 
 %%
 
-esfit(trSignalPowder, @pepper, {Sys, Exp}, {Vary}, FitOpt)
+% Fit = esfit(trSignalPowder/max(trSignalPowder), @pepper, {Sys, Exp}, {Vary}, FitOpt);
+% clf
+% plot(xxSim, ySim/max(ySim), xxSim, trSignalPowder/max(trSignalPowder)) %, xxSim, Fit.fit)
 
-%% My powder average
+%% My powder average (only for equivalent nuclei)
 
-% For some reason I need to multiply by these random factors in order to
-% have a good match with the spectrum simulated using Easyspin
-dip = - Sys.dip*1.4;  % MHz
-JJ = - Sys.J/2.4;  % MHz
-% Ahfi = [9, 9, 12.8];  % MHz, only for A1- radical
-eAngles = Sys.gFrame(1, :);
+tic
+dip = unitconvert(-0.177,'mT->MHz');  % MHz, -4.9 MHz
+JJ = - Sys.J/2;  % MHz
+if isfield(Sys, 'Nucs')
+    nNuc = numel(strsplit(Sys.Nucs, ','));  % Number of equivalent nuclei
+else
+    nNuc = 0;
+end
+nHfiLine = nNuc + 1;
+trlw = Sys.lw;  % ONLY gaussian lw implemented for now
+if isfield(Exp, 'mwFreq')
+    trlw = mt2mhz(trlw);
+end
 
 % Expected input: B0 in mT
 % Output: frequency nu in GHz
@@ -77,6 +104,11 @@ OmegaFunc = @(deltaw, J, d) sqrt(deltaw.^2 + (J*1e-3 + d/2*1e-3).^2);
 %   nVers:  3 x nTheta x nPhi
 %   zD:     3 x 1
 dipFunc = @(dip, nVers, zD) dip*((sum(nVers.*zD)).^2 - 1/3);
+% Input vec of dimensions a1 x a2 x a3 x ... and return vec of dimensions 
+% n x a1 x a2 x a3 x ... with vec(i, :, :, ...) = vec(j, :, :, ...) for
+% every i and j
+insertDimensionPos1 = @(vec, n) ...
+    repmat(reshape(vec, [1, size(vec)]), [n, ones(1, numel(size(vec)))]);
 
 % Theta, phi grid
 thetas = (0:3:180)*pi/180;
@@ -93,51 +125,229 @@ nVers(3, :, :) = cos(thetas')*ones(1, nPhi);
 g2 = squeeze(sqrt( sum( (Sys.g(2, :)'.*nVers).^2)));
 % nVersInFrame1 = pagemtimes(erot(eAngles), nVers);
 % g1 = squeeze(sqrt( sum( (Sys.g(1, :)'.*nVersInFrame1).^2)));
-g1TensorInFrame2 = erot(eAngles)'*diag(Sys.g(1, :))*erot(eAngles);
+euMatrixg1 = erot(Sys.gFrame(1, :));
+g1TensorInFrame2 = euMatrixg1'*diag(Sys.g(1, :))*euMatrixg1;
 g1 = squeeze(sqrt( sum( (pagemtimes(g1TensorInFrame2, nVers)).^2)));
 % Dipolar interaction
 zD = erot(Sys.eeFrame)*[0, 0, 1]';
 dd = squeeze(dipFunc(dip, nVers, zD));
-% Atensor = squeeze(sqrt( sum( (Ahfi'.*nVers).^2)));
 
+%
+% Hyperfine interaction with a number nNuc of equal nuclei
+%
+% For spin 1
+euMatrixHfi1 = erot(Sys.AFrame(1, 1:3));
+Ahfi1TensorInFrame2 = euMatrixHfi1'*diag(Sys.A(1, 1:3))*euMatrixHfi1;
+Ahfi1 = squeeze(sqrt( sum( (pagemtimes(Ahfi1TensorInFrame2', nVers)).^2)));
+% For spin 2
+euMatrixHfi2 = erot(Sys.AFrame(1, 4:6));
+Ahfi2TensorInFrame2 = euMatrixHfi2'*diag(Sys.A(1, 4:6))*euMatrixHfi2;
+Ahfi2 = squeeze(sqrt( sum( (pagemtimes(Ahfi2TensorInFrame2', nVers)).^2)));
+
+%
+% Calculate w0, deltaw, Omega
+%
 w0 = g2wFunc(B0, (g1 + g2)/2);  % GHz
 deltaw = g2wFunc(B0, (g1 - g2)/2);  % GHz
-Omega = OmegaFunc(deltaw, JJ, dd);  % GHz
 
+% Shift resonances due to hyperfine interaction
+Aplus = (Ahfi1 + Ahfi2)/2*1e-3;
+% Make w0 have size [nHfiLine, nTheta, nPhi]
+w0 = insertDimensionPos1(w0, nHfiLine);
+% Make Aplus have size [1, nTheta, nPhi]
+Aplus = reshape(Aplus, [1, size(Aplus)]);
+% Multiply by all possible values of sum of quantum number of the nuclei
+% Final size = [nHfiLine, nTheta, nPhi]
+Aplus = pagemtimes((-nNuc/2:1:nNuc/2)', Aplus);
+w0 = w0 + Aplus;
+% Same for deltaw and Aminus
+Aminus = (Ahfi1 - Ahfi2)/2*1e-3;
+deltaw = insertDimensionPos1(deltaw, nHfiLine); 
+Aminus = reshape(Aminus, [1, size(Aminus)]);
+Aminus = pagemtimes((-nNuc/2:1:nNuc/2)', Aminus);
+deltaw = deltaw + Aminus;
+
+% Pascal factors because some values of the sum of the quantum number come
+% up more than others
+pascalMatrix = pascal(nHfiLine);
+pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);  % Antidiag
+
+Omega = OmegaFunc(deltaw, JJ, insertDimensionPos1(dd, nHfiLine));  % GHz
+
+%
+% Calculate the resonance positions in GHz
+%
 clear('wReson', 'intensityOfLine', 'intensityReson4D')
 
-wReson(1, :, :) = w0 - (JJ - dd)*1e-3 - (Omega);  % w12
-wReson(2, :, :) = w0 + (JJ - dd)*1e-3 - (Omega);  % w34
-wReson(3, :, :) = w0 - (JJ - dd)*1e-3 + (Omega);  % w13
-wReson(4, :, :) = w0 + (JJ - dd)*1e-3 + (Omega);  % w24
+dipInteraction = (JJ - insertDimensionPos1(dd, nHfiLine))*1e-3;  % GHz
+wReson(1, :, :, :) = w0 - dipInteraction - (Omega);  % w12
+wReson(2, :, :, :) = w0 + dipInteraction - (Omega);  % w34
+wReson(3, :, :, :) = w0 - dipInteraction + (Omega);  % w13
+wReson(4, :, :, :) = w0 + dipInteraction + (Omega) ;  % w24
 % BReson = ...
 %     unitconvert(wReson*1e-6, 'MHz->mT', ...
 %     permute(repmat((g1 + g2)/2, [1, 1, 4]), [3, 1, 2]));
 
+% Lineshapes for each transition
+trlw = trlw*1e-3;  % GHz
+trSignal = ...
+    gaussiantransitions(xxSim', wReson, trlw, "fwhm");
 
-intensityOfLine = 1/8*(deltaw.^2)./(Omega.^2);
-intensityReson = repmat([1; -1; 1; -1], [1, nTheta, nPhi]); 
-% This is correct only if the lw is isotropic (not dependent on angles)
-intensityReson = ...
-    intensityReson.*reshape(intensityOfLine, [1, nTheta, nPhi]);  
+% Intensities
+intensityReson = 1/8*(deltaw.^2)./(Omega.^2);
+% Sign of intensityReson alternates for the transitions w12, w34, w13, w24
+intensityReson = [1; -1; 1; -1].*insertDimensionPos1(intensityReson, 4);  
 
-lw1 = repmat(Sys.lw*1e-3, [4, nTheta, nPhi]);
-%lwHfi = repmat(Atensor*1e-3, [4, 1, 1]);
-trSignal = gaussiantransitions(xxSim', wReson, lw1, "fwhm");
-
-trSignal = trSignal.*reshape(intensityReson, [1, size(intensityReson)]);
-trSignalSum = squeeze(sum(trSignal, 2));  % Sum over transitions
+% Normalization to account for possible different linewidths
+intensityNorm = insertDimensionPos1(intensityReson, 1) ./ sum(trSignal);
+trSignal = trSignal.*intensityNorm.*insertDimensionPos1(pascalFactor, 1);
+trSignalSumHfi = squeeze(sum(trSignal, 3));  % Sum over hfi lines
+trSignalSum = squeeze(sum(trSignalSumHfi, 2));  % Sum over transitions
 
 solidAngleWeight = sin(thetas)/sum(sin(thetas')*ones(1, nPhi), 'all');
 trSignalPowder = sum(squeeze(sum(trSignalSum.*solidAngleWeight, 3)), 2);
 
-hold on
+% hold on
+clf
+% plot(xxSim1, ySim/max(ySim), xxSim1 - 0.21, flip(trSignalPowder/max(trSignalPowder)))
+plot(xxSim, ySim/max(ySim), xxSim, trSignalPowder/max(trSignalPowder))
+toc
+
+%% Plot
+
+figure()
+clf
+plot(xxSim, ySim/max(ySim), xxSim, trSignalPowder/max(trSignalPowder))
+ylim(setaxlim(trSignalPowder/max(trSignalPowder), 1.05))
+yticks(0)
+labelaxesfig(gca, 'Microwave frequency / GHz', 'trEPR signal')
+legend('Easyspin', 'Gianluca')
+
+% Using Sys.dip = unitconvert(+0.177,'mT->MHz')/1.5; % MHz
+figPath = "/net/storage/gianlum33/projects/oop_ciss_calculations/images/";
+figPath = figPath + "zech_psiOopEseem_trEPRstickSpectra_EasyspinComparison";
+% exportgraphics(gcf, append(figPath, '.pdf'))
+
+%% Reproduce experimental data
+
+load('/net/storage/gianlum33/projects/zech_psi/data/processed/ZePSI-E-007015.mat');
 % clf
-plot(xxSim, trSignalPowder/max(trSignalPowder))
+% h = ScrollableAxes();
+% plot(h, x{2}, x{1}, y');
+xdata = x{2}/10;  % mT
+ydata = y(500, :);
+ydata = ydata/max(ydata);
+
+%%
+
+figure()
+clf
+plot(xdata, ydata, xxSim1-5.7, flip(trSignalPowder/max(trSignalPowder)))
+xlim([min([xdata; xxSim1'-5.7]), max([xdata; xxSim1'-5.7])])
+ylim(setaxlim(trSignalPowder/max(trSignalPowder), 1.05))
+xticks(339:342)
+yticks(0)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR signal')
+legend('Exp. data', 'Gianluca')
+
+% Using Sys.dip = unitconvert(+0.177,'mT->MHz')/1.5; % MHz
+figPath = "/net/storage/gianlum33/projects/oop_ciss_calculations/images/";
+figPath = figPath + "zech_psiOopEseem_trEPRstickSpectra_ExperimentComparison";
+% exportgraphics(gcf, append(figPath, '.pdf'))
 
 %%
 
 %{
+tic
+clear('Sys', 'Exp', 'Vary', 'VaryExp')
+
+Sys.S = [1/2 1/2];
+Sys.g = [2.0030 2.0026 2.0023; ... % P700+
+         2.0062 2.0051 2.0022];    % A1-
+Sys.gFrame = [-10 -128 -83; ...
+                0    0   0]*pi/180;
+Sys.eeFrame = [0 90 0]*pi/180;  % zD direction is along -x of A1-
+
+Sys.J = unitconvert(-1e-3,'mT->MHz'); % MHz
+Sys.dip = unitconvert(+0.177,'mT->MHz'); % MHz
+% Sys.lw = 9;  % MHz
+Sys.lw = mhz2mt(10);  % MHz
+Sys.initState = 'singlet';
+
+nNuc = 3;
+Sys.Nucs = '1H';
+Sys.Nucs = repmat(append(Sys.Nucs, ', '), [1, nNuc]);
+Sys.Nucs = Sys.Nucs(1:end - 2);
+Sys.A = [zeros(1, 3), 9, 9, 12.8];
+Sys.A = repmat(Sys.A, [nNuc, 1]);
+Sys.AFrame = [0 0 0 -60 90 0]*pi/180;
+Sys.AFrame = repmat(Sys.AFrame, [nNuc, 1]);
+
+% Vary.dip = 1;
+% Vary.J = 0.1;
+Vary.lw = 5;
+% VaryExp.mwCenterSweep = [0.1, 0];
+
+B0 = 345.9;  % mT
+mwFreq = 9.7;
+xxSim = linspace(9.65, 9.75, 2048);  % GHz
+xxSim1 = mhz2mt(xxSim*1e3);  % mT
+
+% Magnetic field sweep
+Exp.CenterSweep = [mean(xdata) + 1.95, max(xdata) - min(xdata)];
+Exp.mwFreq = 9.6;
+Exp.nPoints = numel(xdata);
+Exp.Harmonic = 0;
+
+VaryExp.CenterSweep = [0.5, 0];
+
+FitOpt.Verbosity = 0;
+FitOpt.maxTime = 5;
+
+% 
+% % Opt.GridSize = 10;
+% FitOpt.Verbosity = 1;
+% FitOpt.maxTime = 10;
+% % FitOpt.delta = 1;
+
+yFit0 = pepper(Sys, Exp);
+
+plot(xdata, ydata, xdata, yFit0/max(yFit0))
+toc
+
+%%
+
+esfit(ydata, @pepper, {Sys, Exp}, {Vary, VaryExp}, FitOpt)
+
+%%
+plot(xdata, ydata, xdata, Fit.fit)
+
+%%
+clf
+plot(xxSim, ySim/max(ySim))
+hold on
+plot(xxSim, trSignalSum(:, 1, 1)/max(trSignalSum(:, 1, 1)))
+% plot(xxSim, trSignal(:, 1, 1, 1)/max(trSignal(:, 1, 1, 1)))
+% plot(xxSim, yLow(:, 1, 1, 1)/max(yLow(:, 1, 1, 1)))
+% plot(xxSim, yHigh(:, 1, 1, 1)/max(yHigh(:, 1, 1, 1)))
+
+%%
+ith = 1;
+iph = 1;
+% clf
+hold on
+plot(xxSim, trSignalSum(:, ith, iph)/max(trSignalSum(:, ith, iph)))
+
+%%
+for itrans = 1:4
+    plot(xxSim, trSignal(:, itrans, ith, iph))
+    hold on
+end
+
+
+
+%%
+
 yA = gaussiantransitions(xxSim', wReson, ...
     repmat(Sys.lw*1e-3, [4, nTheta, nPhi]), "fwhm");
 yB = gaussiantransitions(xxSim', wReson, ...
@@ -268,41 +478,18 @@ colorbar
 xlabel('X')
 title('g1')
 
-%%
-
-trxx1 = planck*trxx/bmagn/gfree;
-ith = 10;
-iph = 300;
-ith = 62;
-iph = 185;
-ith = 1;
-iph = 1;
-clf
-for ith = 1:nTheta
-    for iph = 1:nPhi
-for itrans = 1:4
-    plot(xxSim, trSignal(:, itrans, ith, iph))
-    hold on
-end
-    end
-end
-plot(xxSim, trSignal(:, 1, 30, 60))
-xlim([-1, 1]*(Omega(ith, iph) + JJ - dd(ith, iph) + 2*trlw))
-xline(0)
-title(string(dd(ith, iph)))
-
 %}
 
 %%
 
 function y1 = gaussiantransitions(xx, x0, c, mode)
     % Output:
-    % y1:   nAx x 4 x nTheta x nPhi
+    % y1:
     arguments
-        xx  (:, 1) double       % nA x 1 equally spaced values
-        x0  (4, :, :) double    % 4 x nTheta x nPhi
-        c   (4, :, :) double    % 4 x nTheta x nPhi
-        mode string = "var"
+        xx      % [nA, 1] equally spaced values
+        x0      % [4, nHfiLine, nTheta, nPhi]
+        c   (1, 1) double        
+        mode string = "fwhm"
     end
     
     if strcmp(mode, "fwhm")
@@ -310,15 +497,27 @@ function y1 = gaussiantransitions(xx, x0, c, mode)
         c = c/(2*sqrt(2*log(2)));
     end
     
-    nAx = numel(xx);
-    xx = repmat(xx, [1, size(x0)]);
-    x0 = repmat( ...
-            reshape(x0, [1, size(x0)]), ...
-            [nAx, 1, 1, 1]);
-    c = repmat( ...
-            reshape(c, [1, size(c)]), ...
-            [nAx, 1, 1, 1]);
-    y1 = exp(-1/2 * (xx - x0).^2 ./ c.^2);
+    % Adjust dimensions of xx
+    [nAx, dim2] = size(xx);
+    if nAx < dim2
+        xx = xx';
+        nAx = dim2;
+    end
+    if min(size(xx)) ~= 1
+        error('gaussiantransitions: xx should have size [n, 1] or [1, n].')
+    end
+
+    xx = repmat(xx, [1, size(x0, [1, 3, 4])]);  % Size [nAx, 4, nTheta, nPhi]
+    y1 = zeros([nAx, size(x0)]);  % Size [nAx, 4, nHfiLine, nTheta, nPhi]
+    for iNuc = 1:size(x0, 2)  % 1:(nHfiLine)
+        x0_ = squeeze(x0(:, iNuc, :, :));
+        x0_ = repmat( ...
+            reshape(x0_, [1, size(x0_)]), ...
+            [nAx, ones(1, numel(size(x0_)))]);
+        xx_ = xx - x0_;
+        y1(:, :, iNuc, :, :) = exp(-1/2 * (xx_).^2 ./ c^2);
+    end
+    % y1 = exp(-1/2 * (xx - x0).^2 ./ c.^2);
     % Normalization over a range equal to 10*fwhm such that the area is
     % equal to 1
 %     dxx = xx(2) - xx(1);
@@ -331,29 +530,6 @@ function y1 = gaussiantransitions(xx, x0, c, mode)
 
 end
 
-function [y1, yc, yhfi] = gaussiantransitionshfi(xx, x0, c, hfi, mode)
-    % Output:
-    % y1:   nAx x 4 x nTheta x nPhi
-    arguments
-        xx  (:, 1) double       % nA x 1 equally spaced values
-        x0  (4, :, :) double    % 4 x nTheta x nPhi
-        c   (4, :, :) double    % 4 x nTheta x nPhi
-        hfi   (4, :, :) double    % 4 x nTheta x nPhi
-        mode string = "var"
-    end
-    
-    if strcmp(mode, "fwhm")
-        % If the input is fwhm, multiply the input parameters to have var
-        c = c/(2*sqrt(2*log(2)));
-        hfi = hfi/(2*sqrt(2*log(2)));
-    end
-    
-    % Convolution of N_sigma1 and N_sigma2 is N_(sqrt(sigma1^2 + sigma2^2))
-    y1 = gaussiantransitions(xx, x0, sqrt(c.^2 + hfi.^2), "var");
-    yc = gaussiantransitions(xx, x0, c, "var");
-    yhfi = gaussiantransitions(xx, x0, hfi, "var");
-end
-