Clear["Global`*"]
Off[ClebschGordan::phy]
Off[ClebschGordan::tri]
Off[SixJSymbol::tri]

The Hamiltonian for two excitons on a dimer is ℋ = g μ B·((Subscript[S, A])+(Subscript[S, B])) + J (Subscript[S, A])·(Subscript[S, B]) +(Subscript[S, A])·(Subscript[D, A])·(Subscript[S, A])+(Subscript[S, B])·(Subscript[D, B])·(Subscript[S, B]) + (Subscript[S, A])·X·(Subscript[S, B])

Begin by writing the one exciton spin matrices in the Zeeman basis

sx = 1/Sqrt[2]{{0,1,0},{1,0,1},{0,1,0}};
%//MatrixForm
sy=-I/Sqrt[2]{{0,1,0},{-1,0,1},{0,-1,0}};
%//MatrixForm
sz={{1,0,0},{0,0,0},{0,0,-1}};
%//MatrixForm

These matrices are ordered using states  | ++ ⟩, | 0 ⟩, | – ⟩  with quantization along the z-axis of the applied magnetic field

The 9x9 triplet pair matrix forms of the matrices is formed using Kronecker products

sAx=KroneckerProduct[sx,IdentityMatrix[3]];
sAy=KroneckerProduct[sy,IdentityMatrix[3]];
sAz=KroneckerProduct[sz,IdentityMatrix[3]];
sBx=KroneckerProduct[IdentityMatrix[3],sx];
sBy=KroneckerProduct[IdentityMatrix[3],sy];
sBz=KroneckerProduct[IdentityMatrix[3],sz];

The resulting ordering of the vectors are | + + ⟩, | + 0 ⟩, | + – ⟩, | 0 + ⟩, | 0 0 ⟩, | 0 – ⟩, | – + ⟩, | – 0 ⟩, | – – ⟩

The first term of the Hamiltonian represents interaction of the triplets with an external magnetic field B, assuming the B field points along the lab frame z-axis we have

(* Here we assume that the magnetic field and quantization axis points along the lab-orietnted z-axis *)
Hzeeman = Bz(sAz+sBz);
%//MatrixForm

The next term in the Hamiltonian gives the isotropic exchange interaction of the two triplets

(* Assuming that the isotropic exchange is negative (as given by Hund's rule) and J is the absolute value *)
Hexchange = -J(sAx.sBx + sAy.sBy+sAz.sBz);
%//MatrixForm

For now we focus our attention on these first two terms which dominate the Hamiltonian   (Subscript[ℋ, 0])=g μ B·((Subscript[S, A])+(Subscript[S, B])) + J (Subscript[S, A])·(Subscript[S, B])

(* Define (Subscript[H, 0]) and convert to the dimensionless quantity ϵ *)
H0 = FullSimplify[(Hzeeman+Hexchange)/J /.Bz->δ J];
%//MatrixForm

This part of the Hamiltonian is diagonal in the total spin basis, using Clebsch-Gordan coefficients we can write a unitary matrix to transform to this basis

UtotSpin = {{0,0,1/Sqrt[3],0,-1/Sqrt[3],0,1/Sqrt[3],0,0}, (* |0 0⟩ *)
				{0,1/Sqrt[2],0,-1/Sqrt[2],0,0,0,0,0}, (* |1 +1⟩ *)
				{0,0,1/Sqrt[2],0,0,0,-1/Sqrt[2],0,0}, (* |1 0⟩ *)
				{0,0,0,0,0,-1/Sqrt[2],0,1/Sqrt[2],0}, (* |1 -1⟩ *)
				{1,0,0,0,0,0,0,0,0}, (* |2 +2⟩ *)
				{0,1/Sqrt[2],0,1/Sqrt[2],0,0,0,0,0}, (* |2 +1⟩ *)
				{0,0,1/Sqrt[6],0,2/Sqrt[6],0,1/Sqrt[6],0,0}, (* |2 0⟩ *)
				{0,0,0,0,0,1/Sqrt[2],0,1/Sqrt[2],0}, (* |2 -1⟩ *)
				{0,0,0,0,0,0,0,0,1} (* |2 -2⟩ *)};
Transpose[%]//MatrixForm
H0TotSpin=FullSimplify[UtotSpin.H0.Transpose[UtotSpin]];
%//MatrixForm

Plotting the energies, we have

plotEnergies = Simplify[Diagonal[H0TotSpin]];
legend[s_,m_]=Style[StringForm["S=``, M=``",s,m],20];
plotLegend := {legend[0,0],legend[1,+1],legend[1,0],legend[1,-1],legend[2,+2],legend[2,+1],legend[2,0],legend[2,-1],legend[2,-2]};
plotColor = {{Gray,Thick},{Darker[Green],Thick},{Blend[{Gray,Green},0.2],Thick},{Darker[Green],Dashed,Thick},{Cyan,Thick},{Blue,Thick},{Blend[{Gray,Blue},0.2],Thick},{Blue,Dashed,Thick},{Cyan,Dashed,Thick}};
pointSize=0.02;
transitionPoints = {{Red, PointSize@pointSize,Point[{1,2}]},{Red, PointSize@pointSize,Point[{1,1}]},{Red, PointSize@pointSize,Point[{1,0}]},{Red, PointSize@pointSize,Point[{2/3,1-2/3}]},{Red, PointSize@pointSize,Point[{2,1}]},{Red, PointSize@pointSize,Point[{2,-1}]}};
selectionRules = {{Red,Dashed,Line[{{1,-1},{1,0}}]}(*20->21/11*),{Red,Dashed,Line[{{2/3,-1+2/3},{2/3,1-2/3}}]}(*21->22/11*)}(*for parallel choromophores*);

export=False;
Plot[plotEnergies,{δ,0,10},PlotLabel->Style["(Subscript[H, 0]) Energy Levels",Black,28],PlotRange->{{0,2.5},{-3,3}},PlotLegends->Placed[plotLegend,{1.1,0.6}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,28],None},{Style["(Subscript[B, z])/J",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->Join[transitionPoints,selectionRules]*)]
If[export,
Export["C:\\Users\\Sina\\Google Drive\\Research\\Singlet Fission Project\\H0_spectrum_image.png",%,ImageResolution->500],Print["Image not exported"]];

Adding in now the dipole/zfs interaction...

(*In the principal frame of the molecule *)
(* Due to Mathematica reserved functions D -> d, E -> e, can also consider the lowercase quantities to be D/J, E/J respectively*) 
Ddipole[d_,e_] = {{-d/3 + e, 0, 0},{0,-d/3 - e,0},{0,0,2d/3}};

(*Defining active/passive rotation operators*)
Ra[angles_,actsOn_] := EulerMatrix[angles].actsOn.Transpose[EulerMatrix[angles]];
Rp[angles_,actsOn_]:=Ra[{-angles[[3]],-angles[[2]],-angles[[1]]},actsOn]
?EulerMatrix

Visualization of the Euler angles {α,β,γ} -> {ϕ,θ,ψ}

PlotEulerAngles[phi_,theta_,psi_,label_,convention_]:=Module[{xi,eta,zeta,xiprime,etaprime,zetaprime,xprime,yprime,zprime,etheta,epsi,phiarc,thetaarc,psiarc},
(* Start with z rotation *)
xi = RotationMatrix[phi,{0,0,1}].{1,0,0};
eta = RotationMatrix[phi,{0,0,1}].{0,1,0};
zeta = RotationMatrix[phi,{0,0,1}].{0,0,1};

If[convention=="x",
xiprime =RotationMatrix[theta,xi].xi;
etaprime =RotationMatrix[theta,xi].eta;
zetaprime =RotationMatrix[theta,xi].zeta;
etheta = RotationMatrix[theta/2,xi].zeta;
epsi = RotationMatrix[psi/2,zetaprime].xiprime,


(* else the convention == "y" *)
xiprime =RotationMatrix[theta,eta].xi;
etaprime =RotationMatrix[theta,eta].eta;
zetaprime =RotationMatrix[theta,eta].zeta;
etheta = RotationMatrix[theta/2,eta].zeta;
epsi = RotationMatrix[psi/2,zetaprime].etaprime;
];

xprime = RotationMatrix[psi,zetaprime].xiprime;
yprime = RotationMatrix[psi,zetaprime].etaprime;
zprime = RotationMatrix[psi,zetaprime].zetaprime;


phiarc = Table[0.7RotationMatrix[tx,{0,0,1}].{1,0,0},{tx,0,Max[phi,0.01],0.05}];

If[convention=="x",
thetaarc = Table[0.7RotationMatrix[tx,xi].zeta,{tx,0,Max[theta,0.01],0.05}];psiarc = Table[0.7RotationMatrix[tx,zetaprime].xiprime,{tx,0,Max[psi,0.01],0.05}],

(* else the convention == "y" *)
thetaarc = Table[0.7RotationMatrix[tx,eta].zeta,{tx,0,Max[theta,0.01],0.05}];
psiarc = Table[0.7RotationMatrix[tx,zetaprime].etaprime,{tx,0,Max[psi,0.01],0.05}]];

Show[
(* x,y,z axes *)
Graphics3D[{Red,Arrow[Tube[{{0,0,0},{1,0,0}},0.02]]}],
Graphics3D[{Green,Arrow[Tube[{{0,0,0},{0,1,0}},0.02]]}],
Graphics3D[{Blue,Arrow[Tube[{{0,0,0},{0,0,1}},0.02]]}],
(* e1,e2,e3 axes *)
Graphics3D[{Red,Tube[{{0,0,0},xprime},0.02],Sphere[xprime,0.1]}],
Graphics3D[{Green,Tube[{{0,0,0},yprime},0.02],Sphere[yprime,0.1]}],
Graphics3D[{Blue,Tube[{{0,0,0},zprime},0.02],Sphere[zprime,0.1]}],
(* original equator *)
Graphics3D[{Gray,Cylinder[{{0,0,-0.005},{0,0,0.005}},0.7]}],
(* new equator *)
Graphics3D[{White,Cylinder[{-0.005zprime,0.005zprime},0.7]}],
(* intersection of the two *)
Graphics3D[{Black,Tube[{{0,0,0},0.7eta}]}],
(* e1a *)
Graphics3D[{Black,Tube[{{0,0,0},0.7xi}]}],
(* ϕ angle *)
Graphics3D[{Red,Tube[phiarc,0.03]}],
(* θ angle *)
Graphics3D[{Blue,Tube[thetaarc,0.03]}],
(* ψ angle *)
Graphics3D[{Green,Tube[psiarc,0.03]}],
If[label,
Graphics3D[{
Text[Style["x",Italic,FontSize->30],{1.2,0,0}],
Text[Style["y",Italic,FontSize->30],{0,1.2,0}],
Text[Style["z",Italic,FontSize->30],{0,0,1.2}],
Text[Style["x'",Italic,FontSize->30],1.3xprime],
Text[Style["y'",Italic,FontSize->30],1.3yprime],
Text[Style["z'",Italic,FontSize->30],1.3zprime]}],{}],
Graphics3D[{
Text[Style["θ",FontSize->30],0.85 etheta],
Text[Style["ϕ",FontSize->30],{0.85Cos[phi/2],0.85Sin[phi/2],0}],
Text[Style["ψ",FontSize->30],0.85epsi]}],
PlotRange->{{-1.4,1.4},{-1.4,1.4},{-1.4,1.4}},ViewPoint->{2,2,1},ImageSize->{450,480}, SphericalRegion-> True]];

Manipulate[
PlotEulerAngles[phi,theta,psi,label,convention],
Row[{"phi (", Style["ϕ"], ")"}],
{{phi, .5, ""},0,π,.01, Appearance->"Labeled" , ImageSize-> Tiny},
Row[{"theta (", Style["θ"], ")"}],
{{theta, .8, ""},0,π,.01, Appearance->"Labeled", ImageSize-> Tiny },
Row[{"psi (", Style["ψ"], ")"}],
{{psi, 1.2, ""},0,π,.01, Appearance->"Labeled", ImageSize-> Tiny },
Delimiter,
"axes labels",
{{label,False,""},{True,False}},
"convention",
{{convention, "x", ""},{"x"-> Style["x", Italic],"y"-> Style["y", Italic]}},
SaveDefinitions->True, ControlPlacement-> Left]

For linked dimers we first need to rotate molecule B and molecule A into a shared ‘dimer’ frame, this corresponds to a rotation {180, -(180-β)/2,0} for molecule B and {0,-(180-β)/2,0} for molecule A, where β is the bridging angle

(* Switch between 'chasing' orientation of x-axes and orientation where x-axes point towards joined edge of molecules*)
chasing=False;
(*Set of angles that describes angle between B field and shared molecular axis*)
BFieldAngles={ϕ,θ,0};
(*Set of angles for molecule B that describes relation between principle frame of molecules*) 
BAngles = {If[chasing,0,Pi],If[chasing,1,-1](Pi-β)/2,0};
(* Set of angles for molecule A that describes relation between principle frame of molecules *)
AAngles ={0,-(Pi-β)/2,0};

(* Rotate from principle frame (Ddipole[d,e]) to lab frame *)
DBLab =  Rp[BFieldAngles,Ra[BAngles,Ddipole[dB,eB]]];
DBLab = FullSimplify[%];
%//MatrixForm
DALab =  Rp[BFieldAngles,Ra[AAngles,Ddipole[dA,eA]]];
DALab = FullSimplify[%];
%//MatrixForm

Put together the parts to form the dipole term (Subscript[H, dipole])=(Subscript[S, A])·(Subscript[D, A])·(Subscript[S, A])+(Subscript[S, B])·(Subscript[D, B])·(Subscript[S, B])=(Subscript[Σ, αβ])(Subscript[(Subscript[S, A,α])(((Subscript[D, A])^((L)))), αβ])(Subscript[S, A, β])+(Subscript[S, B,α])(Subscript[(((Subscript[D, B])^((L)))), αβ])(Subscript[S, B, β])

S = {sx,sy,sz};
assume = {β∈Reals,β>0,β≤Pi,ϕ∈Reals,ϕ≥ 0,ϕ≤2Pi,θ∈Reals,θ≥ 0,θ<Pi,d>e,d>0,e≥0,d∈Reals,e∈Reals};
Hdipole[ϕ_,θ_,β_] =Simplify[KroneckerProduct[Sum[Sum[DALab[[i]][[j]]*S[[i]].S[[j]],{i,1,3}],{j,1,3}],IdentityMatrix[3]] +KroneckerProduct[IdentityMatrix[3],Sum[Sum[DBLab[[i]][[j]]*S[[i]].S[[j]],{i,1,3}],{j,1,3}]]/.{dA->d,dB->d,eA->e,eB->e},Assumptions->{θ∈ Reals,ϕ∈Reals,β∈Reals,d>0,e<0}];
%//MatrixForm

Alternately we can use the result of the derivation in my notes

Be aware the Mathematica uses the convention WignerD[{j,m1,m2}, ψ, θ ϕ] = (e^(i m1 ψ))(e^(i m2 ϕ))WignerD[{j,m1,m2},0,θ,0] which is opposite the phase of the active rotation matrix element defined in Mueller. Mathematica ALSO uses WignerD[{j,m1,m2},θ] = (Subscript[((d^j)), -m1,-m2])(θ). To fix this I will define my own WignerD matrix elements

?Conjugate

myWignerD[{j_,m1_,m2_},ψ_,θ_,ϕ_]:=WignerD[{j,-m1,-m2},ψ,θ,ϕ]

Dtensor={e,0,Sqrt[2/3]d,0,e};
HZFS[Spinp_,Mp_,Spin_,M_]:=Sum[(-1)^w Sum[Sum[(-1)^(-Mp)Sqrt[5(2Spin+1)(2Spinp+1)]ThreeJSymbol[{2,-w},{Spin,M},{Spinp,-Mp}]SixJSymbol[{1,1,2},{Spinp,Spin,1}]Dtensor[[j+3]]myWignerD[{2,w,q},-BFieldAngles[[3]],-BFieldAngles[[2]],-BFieldAngles[[1]]]((-1)^(Spin+Spinp)myWignerD[{2,q,j},AAngles[[1]],AAngles[[2]],AAngles[[3]]]+myWignerD[{2,q,j},BAngles[[1]],BAngles[[2]],BAngles[[3]]]),{j,-2,2}],{q,-2,2}],{w,-2,2}];
HZFSKoriModified[Spinp_,Mp_,Spin_,M_]:=Sum[(-1)^q Sum[Sum[(-1)^(-Mp)Sqrt[5(2 Spin+1)(2 Spinp+1)]ThreeJSymbol[{2,q},{Spin,M},{Spinp,-Mp}]SixJSymbol[{1,1,2},{Spinp,Spin,1}]Dtensor[[qpp+3]]
myWignerD[{2,-q,qp},-BFieldAngles[[3]],-BFieldAngles[[2]],-BFieldAngles[[1]]]
((-1)^(Spin+Spinp)myWignerD[{2,qp,qpp},AAngles[[1]],AAngles[[2]],AAngles[[3]]]+myWignerD[{2,qp,qpp},BAngles[[1]],BAngles[[2]],BAngles[[3]]]),{qpp,-2,2}],{qp,-2,2}],{q,-2,2}]
HZFSMatrix=Flatten[Table[Reverse[Table[Flatten[Table[Reverse[Table[HZFS[Sp,Mp,S,M],{M,-S,S}]],{S,0,2}]],{Mp,-Sp,Sp}]],{Sp,0,2}],1];
HZFSKoriMatrix=Flatten[Table[Reverse[Table[Flatten[Table[Reverse[Table[HZFSKoriModified[Sp,Mp,S,M],{M,-S,S}]],{S,0,2}]],{Mp,-Sp,Sp}]],{Sp,0,2}],1];

These should give the same eigenvalues ( but not necessarily the same matrix )

HdipoleTotSpin=FullSimplify[UtotSpin.Hdipole[ϕ,θ,β].Transpose[UtotSpin]/.{dA->d,dB->d,eA->e,eB->e},assume];
HZFSKoriMatrix = FullSimplify[HZFSKoriMatrix,Assumptions->assume];
HZFSMatrix=FullSimplify[HZFSMatrix,Assumptions->assume];

(*HZFSMatrix//MatrixForm
HdipoleTotSpin//MatrixForm
FullSimplify[HdipoleTotSpin-HZFSMatrix]//MatrixForm
Simplify[Table[%[[i]][[j]]*Conjugate[%[[i]][[j]]],{i,1,9},{j,1,9}],Assumptions->{e∈Reals,d∈Reals,ϕ∈Reals,θ∈Reals,β∈Reals}]//MatrixForm*)

The full Hamiltonian is then

(*Assume dA, dB and eA, eB are the same (same molecule) *)
H[ϕ_,θ_,β_,δ_]=H0TotSpin+HdipoleTotSpin;
HTest[ϕ_,θ_,β_,δ_]=H0TotSpin+HZFSMatrix;
HKoriTest[ϕ_,θ_,β_,δ_]=H0TotSpin+HZFSKoriMatrix;
FullSimplify[H[0,0,β,δ]]//MatrixForm(*magnetic field along shared molecule/dimer z axis*)

(*Isolate the blocks to be diagonalized *)
U = {{1,0,0,0,0,0,0,0,0},
{0,0,0,1,0,0,0,0,0},
{0,0,0,0,0,1,0,0,0},
{0,0,0,0,0,0,0,1,0},
{0,1,0,0,0,0,0,0,0},
{0,0,1,0,0,0,0,0,0},
{0,0,0,0,1,0,0,0,0},
{0,0,0,0,0,0,1,0,0},
{0,0,0,0,0,0,0,0,1}};
%//MatrixForm
Simplify[U.HTest[ϕ,θ,β,δ].Transpose[U]]//MatrixForm
Hsimp = Simplify[(U.HTest[ϕ,θ,β,δ].Transpose[U])[[1;;4,1;;4]],Assumptions->{ϕ∈Reals,θ∈Reals,β∈Reals,δ∈Reals,δ>0,d>0,e<0}];
%//MatrixForm
eigs=FullSimplify[Eigenvalues[Hsimp/.{β->Pi/2}],Assumptions->{ϕ∈Reals,θ∈Reals,β∈Reals,δ∈Reals,δ>0,d>0,e<0}]

?Chop

FullSimplify[Chop[eigs[[2]]-eigs[[3]]]/.{d->0.1,e->-0.1/4,β->Pi,θ->0},Assumptions->{δ>0,δ<3,ϕ>0,ϕ<2Pi,ϕ∈Reals}]

U.KroneckerProduct[{"|0,0⟩","|1,1⟩","|1,0⟩","|1,-1⟩","|2,2⟩","|2,1⟩","|2,0⟩","|2,-1⟩","|2,-2⟩"},{"⟨0,0|","⟨1,1|","⟨1,0|","⟨1,-1|","⟨2,2|","⟨2,1|","⟨2,0|","⟨2,-1|","⟨2,-2|"}].Transpose[U]//MatrixForm

(* Make it easy to access different blocks of the Hamiltonian *)
blockHamiltonian ={};
vecList={"|0,0⟩","|1,1⟩","|1,0⟩","|1,-1⟩","|2,2⟩","|2,1⟩","|2,0⟩","|2,-1⟩","|2,-2⟩"};
states={};
For[i=1,i<9,i=i+1,
For[j=i+1,j<10,j=j+1,
blockHamiltonian=Append[blockHamiltonian,{{HTest[ϕ,θ,β,δ][[i]][[i]],HTest[ϕ,θ,β,δ][[i]][[j]]},{HTest[ϕ,θ,β,δ][[j]][[i]],HTest[ϕ,θ,β,δ][[j]][[j]]}}] 
]]
For[i=1;k=1,i<9,i=i+1,
For[j=i+1,j<10,j=j+1,
states=Append[states,{vecList[[i]],vecList[[j]],k}];k=k+1
]]
Transpose[states]//MatrixForm

blockHamiltonian[[23]]//MatrixForm
Simplify[blockHamiltonian[[23]][[1]][[1]]-blockHamiltonian[[23]][[2]][[2]]]

FullSimplify[blockHamiltonian[[6]]/.{β->Pi/2,ϕ->Pi/2,θ->Pi/2}]//MatrixForm

Making a plot

export=False;
diagEnergies[d_,e_,β_,δ_] = Diagonal[H[0,0,β,δ]];
legend[s_,m_]=Style[StringForm["S=``, M=``",s,m],20];
plotLegend = {legend[0,0],legend[1,+1],legend[1,0],legend[1,-1],legend[2,+2],legend[2,+1],legend[2,0],legend[2,-1],legend[2,-2]};
plotColor = {{Gray,Thick},{Darker[Green],Thick},{Blend[{Gray,Green},0.2],Thick},{Darker[Green],Dashed,Thick},{Cyan,Thick},{Blue,Thick},{Blend[{Gray,Blue},0.2],Thick},{Blue,Dashed,Thick},{Cyan,Dashed,Thick}};
pointSize=0.02;
transitionPoints = {{Red, PointSize@pointSize,Point[{1.46667,2}]},{Pink, PointSize@pointSize,Point[{0.1,1.06667}]},{Pink, PointSize@pointSize,Point[{3.0333,2}]},{Yellow, PointSize@pointSize,Point[{1.9,2.86667}]},{Yellow, PointSize@pointSize,Point[{2.0333,-1.06663}]},{Yellow, PointSize@pointSize,Point[{1,1.06667}]},{Yellow, PointSize@pointSize,Point[{2.1,1.06667}]},{Yellow, PointSize@pointSize,Point[{1,-0.033333}]},{Black,PointSize@pointSize,Point[{1.03333,2}]},{Black,PointSize@pointSize,Point[{0.633333,0.333337}]}};
selectionRules = {{Red,Dashed,Line[{{1,-1},{1,0}}]}(*20->21/11*),{Red,Dashed,Line[{{2/3,-1+2/3},{2/3,1-2/3}}]}(*21->22/11*)}(*for parallel choromophores*);
Plot[Evaluate[{diagEnergies[0.1,0.1/3,0,δ]}],{δ,0,10},PlotLabel->Style["",Black,28],PlotRange->{{0,3.5},{-3,3}},PlotLegends->Placed[plotLegend,{1.1,0.6}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,28],None},{Style["Bz/J",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor,Epilog->transitionPoints]
If[export,Export["C:\\Users\\Sina\\Google Drive\\Research\\Singlet Fission Project\\H_avoidedCrossingsIdentified_spectrum_image.png",%,ImageResolution->500];,Print["Image not exported"]]

?Grid

d=0.1;
e=-d/4;
y1=Simplify[Eigenvalues[H[0,0,Pi,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y2=Simplify[Eigenvalues[H[0,0,Pi/2,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y3=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y4= Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/2,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
Grid[{{Plot[Evaluate[y1],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.5,-0.5}},PlotLabel->Style["Quintet B||z, β=π",Black],ImageSize->Medium],
Plot[Evaluate[y1],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.5,1.5}},PlotLabel->Style["Triplet B||z, β=π",Black],ImageSize->Medium]},{
Plot[Evaluate[y2],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.5,-0.5}},PlotLabel->Style["Quintet B||z, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y2],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.5,1.5}},PlotLabel->Style["Triplet B||z, β=π/2",Black],ImageSize->Medium]}},ItemSize->Full]
Grid[{{Plot[Evaluate[y1],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.2,-0.8}},PlotLabel->Style["Quintet B||z, β=π",Black],ImageSize->Medium],
Plot[Evaluate[y1],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.8,1.2}},PlotLabel->Style["Triplet B||z, β=π",Black],ImageSize->Medium]},{
Plot[Evaluate[y3],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.2,-0.8}},PlotLabel->Style["Quintet B||y, β=π",Black],ImageSize->Medium],
Plot[Evaluate[y3],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.8,1.2}},PlotLabel->Style["Triplet B||y, β=π",Black],ImageSize->Medium]}},ItemSize->Full]
Grid[{{Plot[Evaluate[y2],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.2,-0.8}},PlotLabel->Style["Quintet B||z, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y2],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.8,1.2}},PlotLabel->Style["Triplet B||z, β=π/2",Black],ImageSize->Medium]},{
Plot[Evaluate[y4],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{-1.2,-0.8}},PlotLabel->Style["Quintet B||y, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y4],{δ,0,0.2},PlotLegends->Automatic,PlotRange->{{0,0.2},{0.8,1.2}},PlotLabel->Style["Triplet B||y, β=π/2",Black],ImageSize->Medium]}},ItemSize->Full](*
FullSimplify[y2[[7]]-y2[[2]]]/.δ->0
FullSimplify[y3[[3]]-y3[[4]]]/.δ->0
FullSimplify[y4[[4]]-y4[[3]]]/.δ->0*)
d=.
e=.

d=0.1;
e=-d/4;
y5=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/3,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y6=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/4,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
d=.
e=.

Grid[{{Plot[Evaluate[y1],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||z, β=π",Black],ImageSize->Medium],
Plot[Evaluate[y2],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||z, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y3],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π",Black],ImageSize->Medium]},{
Plot[Evaluate[y4],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y5],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/3",Black],ImageSize->Medium],
Plot[Evaluate[y6],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/4",Black],ImageSize->Medium]}},ItemSize->Full]

Solve[y1[[5]]==y1[[2]],δ]
FindRoot[Chop[y2[[6]]]-Chop[y2[[2]]],{δ,1}]
Solve[y3[[3]]==y3[[4]],δ]
{FindMaximum[y4[[3]],{δ,1}],FindMinimum[y4[[4]],{δ,1}]}
{FindMaximum[y5[[2]],{δ,1}],FindMinimum[y5[[3]],{δ,1}]}
{FindMaximum[y6[[2]],{δ,1}],FindMinimum[y6[[3]],{δ,1}]}

d=0.01;
e=-d/4;
y7=Simplify[Eigenvalues[H[0,0,Pi,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y8=Simplify[Eigenvalues[H[0,0,Pi/2,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y9=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y10=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/2,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y11=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/3,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
y12=Simplify[Eigenvalues[H[Pi/2,Pi/2,Pi/4,δ]],Assumptions->{d∈Reals, e∈Reals,d>0,e<0}];
d=.
e=.

Grid[{{Plot[Evaluate[y7],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||z, β=π",Black],ImageSize->Medium],
Plot[Evaluate[y8],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||z, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y9],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π",Black],ImageSize->Medium]},{
Plot[Evaluate[y10],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/2",Black],ImageSize->Medium],
Plot[Evaluate[y11],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/3",Black],ImageSize->Medium],
Plot[Evaluate[y12],{δ,0.8,1.2},PlotLegends->Automatic,PlotRange->{{0.8,1.2},{-0.2,0.2}},PlotLabel->Style["Crossing B||y, β=π/4",Black],ImageSize->Medium]}},ItemSize->Full]

Solve[y7[[3]]==y7[[4]],δ]
FindRoot[Chop[y8[[6]]]-Chop[y8[[2]]],{δ,1}]
Solve[y9[[3]]==y9[[4]],δ]
{FindMaximum[y10[[3]],{δ,1}],FindMinimum[y10[[4]],{δ,1}]}
{FindMaximum[y11[[2]],{δ,1}],FindMinimum[y11[[3]],{δ,1}]}
{FindMaximum[y12[[2]],{δ,1}],FindMinimum[y12[[3]],{δ,1}]}

eigs=Eigenvalues[H[0,0,Pi,δ]]/.{d->0.1,e->-0.1/3};
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
bruteforce=Plot[Mapping,{δ,0,10},PlotLabel->Style["Brute Force Method",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,0.05},{-2,2.5}},PlotLegends->Placed[plotLegend,{1.1,0.6}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];
eigs=Eigenvalues[HTest[0,0,Pi,δ]]/.{d->0.1,e->-0.1/3};
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
SinaTest=Plot[Mapping,{δ,0,10},PlotLabel->Style["Sina's Expression",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,2},{-2,2.5}}, Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];
eigs=Eigenvalues[HKoriTest[0,0,Pi,δ]/.{d->0.1,e->-0.1/3}];
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
KoriTest=Plot[Mapping,{δ,0,10},PlotLabel->Style["Kori's Expression",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,2},{-2,2.5}}, Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];

Show[GraphicsGrid[{{bruteforce},{SinaTest,KoriTest}},ImageSize->Full],PlotLabel->Style["Expression Comparison for β = π, B || z",Black,26,FontFamily->"Helvetica"]]

eigs=Eigenvalues[H[Pi/2,Pi/2,Pi/2,δ]]/.{d->0.1,e->-0.1/3};
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
bruteforce=Plot[Mapping,{δ,0,10},PlotLabel->Style["Brute Force Method",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,2},{-2,2.5}},PlotLegends->Placed[plotLegend,{1.1,0.6}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];
eigs=Eigenvalues[HTest[Pi/2,Pi/2,Pi/2,δ]]/.{d->0.1,e->-0.1/3};
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
SinaTest=Plot[Mapping,{δ,0,10},PlotLabel->Style["Sina's Expression",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,2},{-2,2.5}}, Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];
eigs=Eigenvalues[HKoriTest[Pi/2,Pi/2,Pi/2,δ]/.{d->0.1,e->-0.1/3}];
Mapping={eigs[[9]],eigs[[5]],eigs[[1]],eigs[[3]],eigs[[8]],eigs[[4]],eigs[[7]],eigs[[2]],eigs[[6]]};
KoriTest=Plot[Mapping,{δ,0,10},PlotLabel->Style["Kori's Expression",Black,22,FontFamily->"Helvetica"],PlotRange->{{0,2},{-2,2.5}}, Frame->{{True,False},{True,False}},FrameLabel->{{Style["Energy/J",Black,18,FontFamily->"Helvetica"],None},{Style["Bz/J",Black,18,FontFamily->"Helvetica"], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->plotColor(*,Epilog->transitionPoints*)];

Show[GraphicsGrid[{{bruteforce},{SinaTest,KoriTest}},ImageSize->Full],PlotLabel->Style["Expression Comparison for β = π/2, B || y",Black,26,FontFamily->"Helvetica"]]

expression=Evaluate[Eigenvalues[HTest[0,0,0,δ]]-HTest[0,0,0,δ][[7]][[7]]]/.{d->0.1,e->-0.1/4};
parallelList=Table[expression,{δ,0,3,0.01}];
ListPlot[Transpose[%]]
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/fullSpectrum_Eigenvalues_beta180.txt",parallelList,"Table"];

expression=Evaluate[Eigenvalues[HTest[Pi/2,Pi/2,Pi/2,δ]]-HTest[Pi/2,Pi/2,Pi/2,δ][[7]][[7]]]/.{d->0.1,e->-0.1/4};
bentList=Table[expression,{δ,0,3,0.01}];
ListPlot[Transpose[%]]
Export["/Users/sina/Google Drive/Research/Singlet Fission Project/fullSpectrum_Eigenvalues_beta90.txt",bentList,"Table"];

Plot[Evaluate[(Eigenvalues[HTest[Pi/2,Pi/2,Pi/2,δ]]-HTest[Pi/2,Pi/2,Pi/2,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}],{δ,0,2}]
Plot[Evaluate[(Eigenvalues[HTest[Pi,Pi/4,Pi/2,δ]]-HTest[Pi,Pi/4,Pi/2,δ][[7]][[7]])/.{d->0.1,e->-0.1/4}],{δ,0,2}]
(*ListPlot[Transpose[%]]*)

Chop[FullSimplify[HTest[ϕ,θ,β,δ]/.{d->0.1,e->-0.1/4,δ->1}]]//MatrixForm