% This program calculates the gxx,gyy and gzz values of a S=1 iron(IV) system (d4)
% starting from crystal field parameters (Delta and V) and spin orbit 
% coupling constant (lambda). Based on the theory of W. T. Oosterhuis and 
% G. Lang, J. Chem. Phys. 58, 4757 (1973). All the equations and figures 
% present in the comments refer to the work of Bracci et al. Journal 
% Applied Magnetic Resonance 2020


% Variables: V and Delta = Crystal field rhombic and axial terms, respectively; lambda = spin-orbit coupling constant
% NOTE: in the case you are using the programs CpdI_simulations.m or loop_gs_FeIV.m or you have to disable (comment) V and Delta, so they can be properly loaded from the corresponding program.
    V=1000; 
    Delta=3000;
    lambda=400; %cm-1

% Building the crystal field matrix
    Crystal_1e=[0 0 0; 0 Delta-V/2 0; 0 0 Delta+V/2]; % Building the crystal field matrix for one electron on the basis set dxy, dxz, dyz
    Crystal_1h=-Crystal_1e; % Crystal field for one hole in the basis set dxy, dzx, dyz. changing from electron model to hole model.
    Crystal_2h=kron(Crystal_1h,diag([1 1 1]))+kron(diag([1 1 1]),Crystal_1h); %Building the crystal field matrix for two holes in the basis set dxydxy,dxydxz,dxydyz,dxzdyz,dxzdxz,dxzdyz,dyzdxy,dyzdxz,dyzdyz

% Defining the antisymmetric basis set in the (|dxy>,|dzx>,|dzy>)x(|dxy>,|dzx>,|dzy>) basis
    phix=(1i/sqrt(2))*[0 1 0 -1 0 0 0 0 0]'; 
    phiy=(1/sqrt(2))*[0 0 1 0 0 0 -1 0 0]';
    phiz=(1i/sqrt(2))*[0 0 0 0 0 -1 0 1 0]';

% Crystal field contribution of the two-hole system
    Crystal2h_r=[phix'*Crystal_2h*phix phix'*Crystal_2h*phiy phix'*Crystal_2h*phiz;phiy'*Crystal_2h*phix phiy'*Crystal_2h*phiy phiy'*Crystal_2h*phiz;phiz'*Crystal_2h*phix phiz'*Crystal_2h*phiy phiz'*Crystal_2h*phiz]; % Restricted two holes ligand field matrix in the basis sets phix,phiy,phiz
    Crystal=Crystal2h_r-diag([Crystal2h_r(3,3) Crystal2h_r(3,3) Crystal2h_r(3,3)]); % Moving the energy of the ground state phiz to 0. 
    CF=kron(Crystal,diag([1 1 1])); % Building the crystal field matrix taking into account the orbital functions and spin parts, on the basis set of eq. 6.

% Defining spin operators for spin-orbit coupling (angular momentum part)
    I=diag([1 1 1]);
    lx=sop(2,'x');ly=sop(2,'y');lz=sop(2,'z'); %generated by EasySpin
    dxy=(-1i/sqrt(2)).*[1 0 0 0 -1]';dxz=(1/sqrt(2)).*[0 -1 0 1 0]';dyz=(1i/sqrt(2)).*[0 1 0 1 0]';dx2y2=(1/sqrt(2)).*[1 0 0 0 1]';dz2=[0 0 1 0 0]'; % Representation of the d orbitals in the basis set (|+2>,|+1>,|0>,|-1>,|-2>)
    R=[dxy, dxz, dyz, dx2y2, dz2]; % transformation matrix from (|+2>,|+1>,|0>,|-1>,|-2>) representation to (|dxy>,|dxz>,|dyz>,|dx2-y2>,|dz2>) basis set
    lx_1e=R'*lx*R; ly_1e=R'*ly*R; lz_1e=R'*lz*R; % Changing spin operator l from the (|+2>,|+1>,|0>,|-1>,|-2>) representation to (|dxy>,|dxz>,|dyz>,|dx2-y2>,|dz2>) basis set for one electron
    lx_1e=lx_1e(1:3,1:3); ly_1e=ly_1e(1:3,1:3); lz_1e=lz_1e(1:3,1:3); % One electron angular momentum matrices limited to (|dxy>,|dxz>,|dyz>) basis set (t2g subspace)
    L1x=kron(lx_1e,I); L2x=kron(I,lx_1e); L1y=kron(ly_1e,I); L2y=kron(I,ly_1e); L1z=kron(lz_1e,I); L2z=kron(I,lz_1e); % Generating the matrices for two electrons of angular momentum L=2 limited to t2g subspace

% Generating the angular momentum matrix of the individual electrons for the two-electron system in the basis of phix, phiy, phiz restricted to t2g 
    L1x_r=[phix'*L1x*phix phix'*L1x*phiy phix'*L1x*phiz; phiy'*L1x*phix phiy'*L1x*phiy phiy'*L1x*phiz; phiz'*L1x*phix phiz'*L1x*phiy phiz'*L1x*phiz];
    L2x_r=[phix'*L2x*phix phix'*L2x*phiy phix'*L2x*phiz; phiy'*L2x*phix phiy'*L2x*phiy phiy'*L2x*phiz; phiz'*L2x*phix phiz'*L2x*phiy phiz'*L2x*phiz];
    L1y_r=[phix'*L1y*phix phix'*L1y*phiy phix'*L1y*phiz; phiy'*L1y*phix phiy'*L1y*phiy phiy'*L1y*phiz; phiz'*L1y*phix phiz'*L1y*phiy phiz'*L1y*phiz];
    L2y_r=[phix'*L2y*phix phix'*L2y*phiy phix'*L2y*phiz; phiy'*L2y*phix phiy'*L2y*phiy phiy'*L2y*phiz; phiz'*L2y*phix phiz'*L2y*phiy phiz'*L2y*phiz];
    L1z_r=[phix'*L1z*phix phix'*L1z*phiy phix'*L1z*phiz; phiy'*L1z*phix phiy'*L1z*phiy phiy'*L1z*phiz; phiz'*L1z*phix phiz'*L1z*phiy phiz'*L1z*phiz];
    L2z_r=[phix'*L2z*phix phix'*L2z*phiy phix'*L2z*phiz; phiy'*L2z*phix phiy'*L2z*phiy phiy'*L2z*phiz; phiz'*L2z*phix phiz'*L2z*phiy phiz'*L2z*phiz];
    
% Defining spin operators for spin-orbit coupling (Spin part)
    I=diag([1 1]);
    sx=sop(1/2,'x');sy=sop(1/2,'y');sz=sop(1/2,'z'); %generated by EasySpin (Pauli matrices)
    S1x=kron(sx,I); S2x=kron(I,sx); S1y=kron(sy,I); S2y=kron(I,sy); S1z=kron(sz,I); S2z=kron(I,sz);

% Basis for the spin part, the symmetric states
    x1=[1 0 0 0]'; x0=(1/sqrt(2)).*[0 1 1 0]'; x_1=[0 0 0 1]';
    S1x_r=[x1'*S1x*x1 x1'*S1x*x0 x1'*S1x*x_1; x0'*S1x*x1 x0'*S1x*x0 x0'*S1x*x_1; x_1'*S1x*x1 x_1'*S1x*x0 x_1'*S1x*x_1];
    S2x_r=[x1'*S2x*x1 x1'*S2x*x0 x1'*S2x*x_1; x0'*S2x*x1 x0'*S2x*x0 x0'*S2x*x_1; x_1'*S2x*x1 x_1'*S2x*x0 x_1'*S2x*x_1];
    S1y_r=[x1'*S1y*x1 x1'*S1y*x0 x1'*S1y*x_1; x0'*S1y*x1 x0'*S1y*x0 x0'*S1y*x_1; x_1'*S1y*x1 x_1'*S1y*x0 x_1'*S1y*x_1];
    S2y_r=[x1'*S2y*x1 x1'*S2y*x0 x1'*S2y*x_1; x0'*S2y*x1 x0'*S2y*x0 x0'*S2y*x_1; x_1'*S2y*x1 x_1'*S2y*x0 x_1'*S2y*x_1];
    S1z_r=[x1'*S1z*x1 x1'*S1z*x0 x1'*S1z*x_1; x0'*S1z*x1 x0'*S1z*x0 x0'*S1z*x_1; x_1'*S1z*x1 x_1'*S1z*x0 x_1'*S1z*x_1];
    S2z_r=[x1'*S2z*x1 x1'*S2z*x0 x1'*S2z*x_1; x0'*S2z*x1 x0'*S2z*x0 x0'*S2z*x_1; x_1'*S2z*x1 x_1'*S2z*x0 x_1'*S2z*x_1];

% Generating the spin orbit coupling matrix
    SO=-lambda*(kron(L1x_r,S1x_r)+kron(L2x_r,S2x_r)+kron(L1y_r,S1y_r)+kron(L2y_r,S2y_r)+kron(L1z_r,S1z_r)+kron(L2z_r,S2z_r));   % The negative lambda is necessary to change from two-electron system to the two-hole system. 

% Diagonalyze the Hamiltonian to get D and E
    H=CF+SO;        % Hamiltonian in eq. 6
    [Vec,En]=eig(H,'nobalance'); % Calaculating the eigenstates and the energies
    Energy=diag(En); %Generating the energy levels as in figure 3.
    E=-(Energy(3)-Energy(2))/2; %Zero field splitting parameters
    D=(Energy(3)+Energy(2))/2-Energy(1);

% Zeeman perturbation
    psi0=Vec(:,1); psi_1=Vec(:,2); psi1=Vec(:,3); % As defined in eq. 11, psi are the three lowest energy levels. 

% As eig matlab function can change the sign (phase) of the eigenstates, the sign should be changed according to eq. 11
    if psi0(8)<0 
        psi0 = -psi0;
    end
    if psi_1(9)<0
        psi_1 = -psi_1;
    end
    if psi1(9)>0
        psi1 = -psi1;
    end
    
% Writing the crystal field and spin orbit terms already calculated in a spin Hamiltonian zero field splitting term for Seff = 1.
    a1=1.*[0 1 0]'; a3=1/sqrt(2).*[1 0 -1]'; a2=1/sqrt(2).*[1 0 1]'; A=[a1 a2 a3]; % Eigenstates of the ZFS matrix as defined in eq. 9, in the basis set of Seff = 1 chi1, chi0, chi-1. For E<0.

% Defining li and si operators in the basis sets phi x chi (eq. 6) for two electrons
    l1x_r=kron(L1x_r,diag([1 1 1]));   s1x_r=kron(diag([1 1 1]),S1x_r);
    l1y_r=kron(L1y_r,diag([1 1 1]));   s1y_r=kron(diag([1 1 1]),S1y_r);
    l1z_r=kron(L1z_r,diag([1 1 1]));   s1z_r=kron(diag([1 1 1]),S1z_r);
    l2x_r=kron(L2x_r,diag([1 1 1]));   s2x_r=kron(diag([1 1 1]),S2x_r);
    l2y_r=kron(L2y_r,diag([1 1 1]));   s2y_r=kron(diag([1 1 1]),S2y_r);
    l2z_r=kron(L2z_r,diag([1 1 1]));   s2z_r=kron(diag([1 1 1]),S2z_r);

% Calculating the zeeman perturbation in the psi states for two electrons
    t1z=[psi0'*(l1z_r+gfree.*s1z_r)*psi0 psi0'*(l1z_r+gfree.*s1z_r)*psi_1 psi0'*(l1z_r+gfree.*s1z_r)*psi1; psi_1'*(l1z_r+gfree.*s1z_r)*psi0 psi_1'*(l1z_r+gfree.*s1z_r)*psi_1 psi_1'*(l1z_r+gfree.*s1z_r)*psi1; psi1'*(l1z_r+gfree.*s1z_r)*psi0 psi1'*(l1z_r+gfree.*s1z_r)*psi_1 psi1'*(l1z_r+gfree.*s1z_r)*psi1];
    t1y=[psi0'*(l1y_r+gfree.*s1y_r)*psi0 psi0'*(l1y_r+gfree.*s1y_r)*psi_1 psi0'*(l1y_r+gfree.*s1y_r)*psi1; psi_1'*(l1y_r+gfree.*s1y_r)*psi0 psi_1'*(l1y_r+gfree.*s1y_r)*psi_1 psi_1'*(l1y_r+gfree.*s1y_r)*psi1; psi1'*(l1y_r+gfree.*s1y_r)*psi0 psi1'*(l1y_r+gfree.*s1y_r)*psi_1 psi1'*(l1y_r+gfree.*s1y_r)*psi1];
    t1x=[psi0'*(l1x_r+gfree.*s1x_r)*psi0 psi0'*(l1x_r+gfree.*s1x_r)*psi_1 psi0'*(l1x_r+gfree.*s1x_r)*psi1; psi_1'*(l1x_r+gfree.*s1x_r)*psi0 psi_1'*(l1x_r+gfree.*s1x_r)*psi_1 psi_1'*(l1x_r+gfree.*s1x_r)*psi1; psi1'*(l1x_r+gfree.*s1x_r)*psi0 psi1'*(l1x_r+gfree.*s1x_r)*psi_1 psi1'*(l1x_r+gfree.*s1x_r)*psi1];
    t2z=[psi0'*(l2z_r+gfree.*s2z_r)*psi0 psi0'*(l2z_r+gfree.*s2z_r)*psi_1 psi0'*(l2z_r+gfree.*s2z_r)*psi1; psi_1'*(l2z_r+gfree.*s2z_r)*psi0 psi_1'*(l2z_r+gfree.*s2z_r)*psi_1 psi_1'*(l2z_r+gfree.*s2z_r)*psi1; psi1'*(l2z_r+gfree.*s2z_r)*psi0 psi1'*(l2z_r+gfree.*s2z_r)*psi_1 psi1'*(l2z_r+gfree.*s2z_r)*psi1];
    t2y=[psi0'*(l2y_r+gfree.*s2y_r)*psi0 psi0'*(l2y_r+gfree.*s2y_r)*psi_1 psi0'*(l2y_r+gfree.*s2y_r)*psi1; psi_1'*(l2y_r+gfree.*s2y_r)*psi0 psi_1'*(l2y_r+gfree.*s2y_r)*psi_1 psi_1'*(l2y_r+gfree.*s2y_r)*psi1; psi1'*(l2y_r+gfree.*s2y_r)*psi0 psi1'*(l2y_r+gfree.*s2y_r)*psi_1 psi1'*(l2y_r+gfree.*s2y_r)*psi1];
    t2x=[psi0'*(l2x_r+gfree.*s2x_r)*psi0 psi0'*(l2x_r+gfree.*s2x_r)*psi_1 psi0'*(l2x_r+gfree.*s2x_r)*psi1; psi_1'*(l2x_r+gfree.*s2x_r)*psi0 psi_1'*(l2x_r+gfree.*s2x_r)*psi_1 psi_1'*(l2x_r+gfree.*s2x_r)*psi1; psi1'*(l2x_r+gfree.*s2x_r)*psi0 psi1'*(l2x_r+gfree.*s2x_r)*psi_1 psi1'*(l2x_r+gfree.*s2x_r)*psi1];

% Summing up the contribution of the first and second electron 
    tz = t1z + t2z;
    ty = t1y + t2y;
    tx = t1x + t2x;
    
% Terms proportional to magnetic field B in zeeman perturbation in chi1, chi0, chi-1
    gz=A*tz*A';
    gy=A*ty*A';
    gx=A*tx*A';

    Sx=sop(1,'x'); Sy=sop(1,'y'); Sz=sop(1,'z'); %generated by EasySpin
    
% Finding the gxx gyy and gzz
    gxx = gx(2)./Sx(2);
    gyy = gy(2)./Sy(2);
    gzz = gz(1)./Sz(1);
 
    clear A	CF	Crystal	Crystal2h_r	Crystal_1e	Crystal_1h	Crystal_2h	En	En_r	Energy	H	I	L1x	L1x_r	L1y	L1y_r	L1z	L1z_r	L2x	L2x_r	L2y	L2y_r	L2z	L2z_r	R	S1x	S1x_r	S1y	S1y_r	S1z	S1z_r	S2x	S2x_r	S2y	S2y_r	S2z	S2z_r	SO	Sx	Sy	Sz	Vec	ZFS	a1	a2	a3	dx2y2	dxy	dxz	dyz	dz2	gx	gy	gz	l1x_r	l1y_r	l1z_r	l2x_r	l2y_r	l2z_r	lx	lx_1e	ly	ly_1e	lz	lz_1e	phix	phiy	phiz	psi0	psi1	psi_1	s1x_r	s1y_r	s1z_r	s2x_r	s2y_r	s2z_r	sx	sy	sz	t1x	t1y	t1z	t2x	t2y	t2z	tx	ty	tz	x0	x1	x_1
