% This program simulates the spectra of the proteins taken into account in 
% Bracci et al. (Journal Applied Magnetic Resonance 2020). This program calls 
% gs_from_crystal_field.m because it takes the g-values of Fe(IV) and the 
% zero field splitting parameters (D and E) from the program gs_from_crystal_field.m. 
% You have to introduce crystal field parameters Delta and V here as input values. 
% NOTE: to properly works you have to delete (or comment) the Delta and V variables inside 
% gs_from_crystal_field.m. In that way they can be properly loaded from
% here. 
% All the simulations have been performed using EasySpin

%% HRP using a J distribution
clear all

V=0.1; 
Delta=1830; %corresponding to D = 26 cm-1 

gs_from_crystal_field

%limit the parameters to the fourth decimal place
D = round(D,4); 
E = round(E,4);
gFex = round(gxx,4);
gFey = round(gyy,4);
gFez = round(gzz,4);

J = [4 -2 -2]*1e-6* boltzm/planck; % conversion from K to MHz % Same values of Schulz_1979.

HRP.S = [1 1/2];  % two electron spins
grad = gfree;
HRP.g = [ gFex gFey gFez; grad grad grad];
HRP.lw = 5;
HRP.D = [D*29979 E*29979; 0 0]; %conversion from cm-1 to MHz
ExpHRP= struct('mwFreq',9.2,'Range',[0 700],'nPoints',3000,'Temperature',4);
ExpHRP.Harmonic =0;

spectrum = 0;
x = linspace(-1,3,200); ya = gaussian(x,1,1.1); %x and ya determines the J distribution. Same standard deviation as Schulz 1979
int = zeros(length(x),ExpHRP.nPoints);
for i = 1:length(x)
    HRP.ee = J*x(i);
    [B,int(i,:)]=pepper(HRP,ExpHRP); 
    spectrum = int(i,:)*ya(i)+spectrum;
end

figure(1004) %figure 10d
hold on
plot(B,spectrum/max(spectrum)*1.8+0.5,'color',[0.25, 0.25, 0.25],'LineWidth', 2)
plot(B,int(60,:),B,int(70,:),B,int(80,:),B,int(120,:),'LineWidth', 2)
xlabel('[mT]')
xlim([240 600])
xticks([240 360 480 600 ])
set(gca,'TickDir', 'out','box', 'off','ytick',[], 'Linewidth', 1,'FontSize',15); 
legend('box', 'off', 'FontSize', 15,'Location','northeast');
legend('Simulation','red','yellow','purple','green')
hold off
%% CPO using a J distribution
clear all

V=120; %corresponding to E/D = -0.035
Delta=1280; %corresponding to D = 39 cm-1 
rarioJD = 1.02; %J/D ratio

gs_from_crystal_field

%limit the parameters to the fourth decimal place
D = round(D,4); 
E = round(E,4);
gFex = round(gxx,4);
gFey = round(gyy,4);
gFez = round(gzz,4);

J = rarioJD*D;

CPO.S = [1 1/2];  % two electron spins
grad = 2.08;
CPO.g = [gFex gFey gFez; grad grad grad];
CPO.lw = [0.9 0.7] ; 
CPO.D = [D*29979 E*29979; 0 0]; %from cm-1 to MHz
ExpCPO= struct('mwFreq',9.5,'Range',[250 600],'nPoints',3000,'Temperature',4); %9.5GHz Exp.mwfreq from Rutter1984
ExpCPO.Harmonic =0;

spectrum = 0;
x = linspace(0.9,1.1,500); ya = gaussian(x,1,0.052); %x and ya determines the J distribution. Same standard deviation as Rutter 1984
int = zeros(length(x),ExpCPO.nPoints);
for i = 1:length(x)
    CPO.J = J*x(i)*29979; %from cm-1 to MHz
    [B,int(i,:)]=pepper(CPO,ExpCPO); 
    spectrum = int(i,:)*ya(i)+spectrum;
end

figure(1104) %figure 11d
plot(B,spectrum/max(spectrum)*1.8+0.5,'color',[0.25, 0.25, 0.25],'LineWidth', 2)
xlabel('[mT]')
xlim([320 450])
xticks([320 339.4 392 413.9 450]) %B=339.4 392 413.9 correspond to a g = 2 1.73 1.64
set(gca,'TickDir', 'out','box', 'off','ytick',[], 'Linewidth', 1,'FontSize',15);
legend('box', 'off', 'FontSize', 15,'Location','northeast');
legend('Simulation')
%% P450 using a J distribution
clear all

V=0.1; 
Delta=1370; %corresponding to D = 36 cm-1 
rarioJD = 1.25; %J/D ratio

gs_from_crystal_field

%limit the parameters to the fourth decimal place
D = round(D,4); 
E = round(E,4);
gFex = round(gxx,4);
gFey = round(gyy,4);
gFez = round(gzz,4);

J = rarioJD*D;

P450.S = [1 1/2];  % two electron spins
grad = 2.08;
P450.g = [ gFex gFey gFez; grad grad grad];
P450.lw = [0 1.8];
P450.D = [D*29979 E*29979; 0 0]; %from cm-1 to MHz
ExpP450= struct('mwFreq',9.25,'Range',[250 500],'nPoints',3000,'Temperature',4);
ExpP450.Harmonic =1;

spectrum = 0;
x = linspace(0.8,1.2,600); ya = gaussian(x,1,0.13);%x and ya determines the J distribution.
int = zeros(length(x),ExpP450.nPoints);
for i = 1:length(x)
    P450.ee = J*x(i)*29979; %from cm-1 to MHz
    [B,int(i,:)]=pepper(P450,ExpP450); 
    spectrum = int(i,:)*ya(i)+spectrum;
end

figure(1203) %figure 12c
plot(B,spectrum/max(spectrum)*1.8+3,'color',[0.25, 0.25, 0.25],'LineWidth', 2)
xlabel('[mT]')
xlim([300 400])
xticks([300 350 400])
set(gca,'TickDir', 'out','box', 'off','ytick',[], 'Linewidth', 1,'FontSize',15);
legend('box', 'off', 'FontSize', 15,'Location','northeast');
legend('Simulation')
%% Catalase using a J distribution
clear all

V=0.1; 
Delta=2680; %corresponding to D = 17 cm-1 
rarioJD = -0.5; %J/D ratio

gs_from_crystal_field

%limit the parameters to the fourth decimal place
D = round(D,4); 
E = round(E,4);
gFex = round(gxx,4);
gFey = round(gyy,4);
gFez = round(gzz,4);

J = rarioJD*D; 

Cat.S = [1 1/2];  % two electron spins
grad = gfree;
Cat.g = [ gFex gFey gFez; grad grad grad];
Cat.lw = 10;
Cat.D = [D*29979 E*29979; 0 0];
Cat.J = J*29979;
CatExp= struct('mwFreq',9.525,'Range',[0 600],'nPoints',1500,'Temperature',2); %9.525 GHz Exp.mwfreq from Benecky 1993
CatExp.Harmonic =0;

spectrum = 0;
x = linspace(0,2,600); ya = gaussian(x,1,0.8);  %x and ya determines the J distribution.
int = zeros(length(x),CatExp.nPoints);
for i = 1:length(x)
    Cat.J = J*x(i)*29979;
    [B,int(i,:)]=pepper(Cat,CatExp);
    spectrum = int(i,:)*ya(i)+spectrum;
end

figure(1304) %figure 13d
plot(B,spectrum/max(spectrum)*1.8+0.5,'color',[0.25, 0.25, 0.25],'LineWidth', 2)
xlabel('[mT]')
xlim([140 420])
xticks([140 204.982 220 300 380]) % B=204.982 corresponds to g = 3.32
set(gca,'TickDir', 'out','box', 'off','ytick',[], 'Linewidth', 1,'FontSize',15);
legend('box', 'off', 'FontSize', 15,'Location','northeast');
legend('Simulation')

