(* Cargado de integrales mono y bi electrónicas en la base atómica. Aquí 1==a y 2==b.*)

h=(({{-2.6442, -1.5113}, {-1.5113, -1.7201}}));
hh[a_,b_]:=h[[a,b]]

BI11=(({{1.0547, 0.4744}, {0.4744, 0.2469}}));
BI12=(({{0.4744, 0.5664}, {0.2469, 0.3504}}));
BI21=(({{0.4744, 0.2469}, {0.5664, 0.3504}}));
BI22=(({{0.2469, 0.3504}, {0.3504, 0.6250}}));
BIP={BI11,BI12, BI21, BI22};
alfa[x_,y_]:=2*(x-1)+y
BI[a_,b_,c_,d_]:=BIP[[alfa[a,b]]][[c,d]]

(*INICIO DEL CÁLCULO AUTOCONSISTENTE*)

(*Acá van los coeficientes de los orbitales moleculares en la base atómica, para la pasada inicial.*)

Coef1={1,0}
Coef2={0,1}
Coef1SQ=(Sqrt[((Coef1[[1]])^2)+2*0.5784*Coef1[[1]]*Coef1[[2]]+((Coef1[[2]])^2)]);
Coef2SQ=(Sqrt[((Coef2[[1]])^2)+2*0.5784*Coef2[[1]]*Coef2[[2]]+((Coef2[[2]])^2)]);
C1={Coef1[[1]]/Coef1SQ,Coef1[[2]]/Coef1SQ}
C2={Coef2[[1]]/Coef2SQ,Coef2[[2]]/Coef2SQ}

(**Con esa elección inicial, armo el operador de Fock y calculo sus elementos de matriz en la base atómica. Luego multiplico por la inversa de la matriz de overlap y calculo autovalores/autovectores del sistema. Los autovalores son las energías orbitales asociadas a ese paso iterativo; los autovectores son los orbitales moleculares (también para ese paso de la iteración, no son los finales hasta que no converja)*)

AA=Eigensystem[Inverse[(({{1, 0.5784}, {0.5784, 1}}))].(({{hh[1,1]+(∑_(i=1)^2)(∑_(j=1)^2)(2C1[[i]]*BI[1,i,1,j]*C1[[j]]-C1[[i]]*BI[1,i,j,1]*C1[[j]]), hh[1,2]+(∑_(i=1)^2)(∑_(j=1)^2)(2C1[[i]]*BI[1,i,2,j]*C1[[j]]-C1[[i]]*BI[1,i,j,2]*C1[[j]])}, {hh[2,1]+(∑_(i=1)^2)(∑_(j=1)^2)(2C1[[i]]*BI[2,i,1,j]*C1[[j]]-C1[[i]]*BI[2,i,j,1]*C1[[j]]), hh[2,2]+(∑_(i=1)^2)(∑_(j=1)^2)(2C1[[i]]*BI[2,i,2,j]*C1[[j]]-C1[[i]]*BI[2,i,j,2]*C1[[j]])}}))]
Coef1=AA[[2]][[1]];
Coef2=AA[[2]][[2]];
Coef1SQ=(Sqrt[((Coef1[[1]])^2)+2*0.5784*Coef1[[1]]*Coef1[[2]]+((Coef1[[2]])^2)]);
Coef2SQ=(Sqrt[((Coef2[[1]])^2)+2*0.5784*Coef2[[1]]*Coef2[[2]]+((Coef2[[2]])^2)]);
(*Estos son los coeficientes de los orbitales moleculares en la base atómica, ya normalizados adecuadamente*)
C1={Coef1[[1]]/Coef1SQ,Coef1[[2]]/Coef1SQ}
C2={Coef2[[1]]/Coef2SQ,Coef2[[2]]/Coef2SQ}

(*Si vuelven a calcular -apretando Shift+Enter en la última línea-, arma el operador de Fock colocando los electrones en los orbitales moleculares de menor energía y repite el procedimiento.*)