% This program simulates the interaction of the Fe(IV) S=1 and the free
% radical S=1/2. It makes the figures 6 (interaction of Fe(IV) with a free
% radical) and 7 (effective g values of Compound I as a funcition of J/D).
% Figura 601 = 6a, Figura 603 = 6c, Figura 701 = 7a, Figura 702 = 7b.  All 
% the equations and figures present in the comments refer to the work of 
% Bracci et al. Journal Applied Magnetic Resonance 2020

% Zero field splitting paramenters (D, E) and g values are taken from the
% program gs_from_crystal_field. Crystal field parameters and spin orbit
% copling constant should be introduced in that program. 
    close all
    clear all
    gs_from_crystal_field

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

    n=1000; %number of points used for the J/D ratio
    ratio=logspace(-2,1,n); %J/D ratio in a logaritmis scale. Change this values to change x scale of the figures 

% Spin operators for S(Fe)=1. generated by EasySpin (Pauli matrices)
    Sx=sop(1,'x');
    Sy=sop(1,'y');
    Sz=sop(1,'z');
% Spin operators for S(rad)=1/2. generated by EasySpin (Pauli matrices)
    sx=sop(1/2,'x');
    sy=sop(1/2,'y');
    sz=sop(1/2,'z');
% Combining the spin operators of Fe(IV) and the free radical
    Ssx=kron(Sx,sx);
    Ssy=kron(Sy,sy);
    Ssz=kron(Sz,sz);
    
    for j=1:n
        J=ratio(j)*D; %D is kept constant and J is varied according to J/D ratio
        Jx=J; Jy=J; Jz=J; %assuming an isotropic J tensor
        HSS=Jx.*Ssx+Jy.*Ssy+Jz.*Ssz; %Generating the matrix for the spin-spin interaction of the Fe-radical couppled system. 
        HZFS=kron([D/3 0 E;0 -2*D/3 0;E 0 D/3],diag([1 1])); %generating the ZFS matrix starting from the D and E parameters generated by gs_from_crystal_field
        H_s=HZFS+HSS; %generating the spin Hamiltonian matrix
        H_s = [H_s(:,1) H_s(:,4) H_s(:,5) H_s(:,6) H_s(:,3) H_s(:,2)]; H_s = [H_s(1,:); H_s(4,:); H_s(5,:); H_s(6,:); H_s(3,:); H_s(2,:)];% reordering the matrix to get eq. 23 in the basis set [phi1alpha,phi0beta,phi-1alpha,phi-1beta,phi0alpha,phi1beta] 
        [Vec,En]=eig(H_s,'vector','nobalance'); %Calculating the egienstates and the energies 
        % As eig matlab function can change the sign (phase) of the eigenstates, the sign should be restored according to eq. 24
        if Vec(2,1) <0 
            Vec = -Vec;
        end
        E1(1,j)=En(1); E2(1,j)=En(2); E3(1,j)=En(3); E4(1,j)=En(4); E5(1,j)=En(5); E6(1,j)=En(6) ;
        V1(j,:)=Vec(:,1)';  V2(j,:)=Vec(:,2)';
    end

% Terms proportional to magnetic field B in zeeman perturbation 
    tz=(kron(gzz.*Sz,diag([1 1]))+kron(diag([1 1 1]),gfree.*sz));
    tx=(kron(gxx.*Sx,diag([1 1]))+kron(diag([1 1 1]),gfree.*sx));
    ty=(kron(gyy.*Sy,diag([1 1]))+kron(diag([1 1 1]),gfree.*sy));
    
% reordering the matrices in the basis set [phi1alpha,phi0beta,phi-1alpha,phi-1beta,phi0alpha,phi1beta] 
    tz = [tz(:,1) tz(:,4) tz(:,5) tz(:,6) tz(:,3) tz(:,2)]; tz = [tz(1,:); tz(4,:); tz(5,:); tz(6,:); tz(3,:); tz(2,:)];
    tx = [tx(:,1) tx(:,4) tx(:,5) tx(:,6) tx(:,3) tx(:,2)]; tx = [tx(1,:); tx(4,:); tx(5,:); tx(6,:); tx(3,:); tx(2,:)];
    ty = [ty(:,1) ty(:,4) ty(:,5) ty(:,6) ty(:,3) ty(:,2)]; ty = [ty(1,:); ty(4,:); ty(5,:); ty(6,:); ty(3,:); ty(2,:)];
    
% calculating the g eff values
    gz=2*diag(V1*tz*V1');
    gx=2*diag(V2*tx*V1');
    gy=2*1i*diag(V2*ty*V1');

% Figure 601 (fig. 6a): J/D from negative to positive vs energy levels
    figure(601) %figure 6a.
    title('Figure 6a. enrgy vs J/D')
    hold on
    plot(ratio,E1,'k',ratio,E2,'k',ratio,E3,'k',ratio,E4,'k',ratio,E5,'k',ratio,E6,'k')
    hold off

% figure  701 (fig. 7a) (J<0), 702 (fig. 7b) (J>0) g_eff vs J/D
    figure(702) %figure 7b
    title('Figure 7b. J>0')
    hold on
    plot(ratio,abs(gz),'color',[0, 0.4470, 0.7410],'LineWidth', 2)
    plot(ratio,abs(gx),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2)
    plot(ratio,abs(gy),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2)
    ylabel('g-values')
    xlabel('|J|/D')
    xlim([0.01 10])
    legend('g_z','g_x','g_y')
    % legend('g_/_/','g_\perp'); %if E=0, gx = gy = gperp
    xticks([0.01 0.1 1 10])
    xticklabels([0.01 0.1 1 10])
    set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15,'XScale', 'log');
    legend('box', 'off', 'FontSize', 15,'Location','southwest');

% Figura 603 (fig. 6c): coefficents vs J/D. 
    figure(603) % figure 6c
    title('Figure 6c. coefficients vs J/D')
    hold on
    plot(ratio,V1(:,1:3),'k', 'LineWidth', 1)
    hold off

    clear Vec En H_s H HZFS HSS Jx Jy Jz J V1 V2 E1 E2 E3 E4 E5 E6 tz tx ty gz gx gy
for j=1:n
        J=-ratio(j)*D;
        Jx=J; Jy=J; Jz=J; %assuming an isotropic J tensor
        HSS=Jx.*Ssx+Jy.*Ssy+Jz.*Ssz; %Generating the matrix for the spin-spin interaction of the Fe-radical couppled system. 
        HZFS=kron([D/3 0 E;0 -2*D/3 0;E 0 D/3],diag([1 1])); %generating the ZFS matrix starting from the D and E parameters generated by gs_from_crystal_field
        H_s=HZFS+HSS; %generating the spin Hamiltonian matrix
        H_s = [H_s(:,1) H_s(:,4) H_s(:,5) H_s(:,6) H_s(:,3) H_s(:,2)]; H_s = [H_s(1,:); H_s(4,:); H_s(5,:); H_s(6,:); H_s(3,:); H_s(2,:)];% reordering the matrix to get eq. 23
        [Vec,En]=eig(H_s,'vector','nobalance'); %Calculating the egienstates and the energies 
        % As eig matlab function can change the sign (phase) of the eigenstates, the sign should be restored according to eq. 24
        if Vec(2,1) <0 
            Vec = -Vec;
        end
        E1(1,j)=En(1); E2(1,j)=En(2); E3(1,j)=En(3); E4(1,j)=En(4); E5(1,j)=En(5); E6(1,j)=En(6) ;
        V1(j,:)=Vec(:,1)';  V2(j,:)=Vec(:,2)';
end

% Terms proportional to magnetic field B in zeeman perturbation 
    tz=(kron(gzz.*Sz,diag([1 1]))+kron(diag([1 1 1]),gfree.*sz));
    tx=(kron(gxx.*Sx,diag([1 1]))+kron(diag([1 1 1]),gfree.*sx));
    ty=(kron(gyy.*Sy,diag([1 1]))+kron(diag([1 1 1]),gfree.*sy));
    
% reordering the matrices in the basis set [phi1alpha,phi0beta,phi-1alpha,phi-1beta,phi0alpha,phi1beta] 
    tz = [tz(:,1) tz(:,4) tz(:,5) tz(:,6) tz(:,3) tz(:,2)]; tz = [tz(1,:); tz(4,:); tz(5,:); tz(6,:); tz(3,:); tz(2,:)];
    tx = [tx(:,1) tx(:,4) tx(:,5) tx(:,6) tx(:,3) tx(:,2)]; tx = [tx(1,:); tx(4,:); tx(5,:); tx(6,:); tx(3,:); tx(2,:)];
    ty = [ty(:,1) ty(:,4) ty(:,5) ty(:,6) ty(:,3) ty(:,2)]; ty = [ty(1,:); ty(4,:); ty(5,:); ty(6,:); ty(3,:); ty(2,:)];
        
% calculating the g eff values
    gz=2*diag(V1*tz*V1');
    gx=2*diag(V2*tx*V1');
    gy=2*1i*diag(V2*ty*V1');
    
% Figure 601 (fig. 6a): J/D from negative to positive vs energy levels
    figure(601) %figure 6a
    hold on
    plot(-ratio,E1,'k',-ratio,E2,'k',-ratio,E3,'k',-ratio,E4,'k',-ratio,E5,'k',-ratio,E6,'k')
    hold off
    set(gca,'TickDir', 'out','box', 'off','ytick',[], 'Linewidth', 1);

% figure  701 (fig. 7a) (J<0), 702 (fig. 7b) (J>0) g_eff vs J/D
    figure(701) %figure 7a
    title('Figure 7a. J<0')
    hold on
    plot(ratio,abs(gz),'color',[0, 0.4470, 0.7410],'LineWidth', 2)
    plot(ratio,abs(gx),'color',[0.8500, 0.3250, 0.0980],'LineWidth', 2)
    plot(ratio,abs(gy),'color',[0.4660, 0.6740, 0.1880],'LineWidth', 2)
    ylabel('g-values')
    xlabel('|J|/D')
    xlim([0.01 10])
    legend('g_z','g_x','g_y')
    % legend('g_/_/','g_\perp'); %if E=0, gx = gy = gperp
    xticks([0.01 0.1 1 10])
    xticklabels([0.01 0.1 1 10])
    set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15,'XScale', 'log');
    legend('box', 'off', 'FontSize', 15,'Location','northwest');
    
% Figura 603 (fig. 6c): coefficents vs J/D. 
    figure(603) % figure 6c
    hold on
    plot(-ratio,V1(:,1:3),'k', 'LineWidth', 1)
    hold off
    yline(-sqrt(2/3),'--')
    yline(sqrt(2/3),'--')
    yline(sqrt(1/3),'--')
    yticks([-sqrt(2/3) sqrt(1/3) sqrt(2/3)])
    yticklabels({'-\surd2/3', '\surd1/3', '\surd2/3'})
    ylabel('Coefficients')
    xlabel('J/D')
    set(gca,'TickDir', 'out','box', 'off', 'Linewidth', 1,'FontSize',15);
