RAB := 1.3
rA := {0, -RAB/2, 0}
rB := {0, RAB/2, 0}
r1 := {x1, y1, z1}
r2 := {x2, y2, z2}
R1A := Sqrt[x1^2+(-RAB/2-y1)^2 + z1^2];
R1B := Sqrt[x1^2+(RAB/2-y1)^2 + z1^2];
R2A := Sqrt[x2^2+(-RAB/2-y2)^2 + z2^2];
R2B := Sqrt[x2^2+(RAB/2-y2)^2 + z2^2];
R12 := Sqrt[(x1-x2)^2+ (y1-y2)^2+(z1-z2)^2]
R1Ae := Norm[rA-r1]
R1Be := Norm[rB-r1]
R2Ae := Norm[rA-r2]
R2Be := Norm[rB-r2]
R12e := Norm[r1 -r2]
ψA1 := 1/Sqrt[Pi]*Exp[-R1A]
ψA2 :=1/Sqrt[Pi]*Exp[-R2A]
ψB1 :=1/Sqrt[Pi]*Exp[-R1B]
ψB2 :=1/Sqrt[Pi]*Exp[-R2B]

int := 5

S := NIntegrate[ψA1*ψB1, {x1, -int, int},{y1, -int, int},{z1, -int, int},PrecisionGoal->2]

norm := 2*(1+S^2)

Cb := -NIntegrate[ψA1^2*1/R1B, {x1, -int, int},{y1, -int, int},{z1, -int, int}, Method->MultiDimensional, PrecisionGoal->2]

EWW := NIntegrate[ψA1^2*ψB2^2*1/R12, {x1, -int, int},{y1, -int, int},{z1, -int, int}, {x2, -int, int},{y2, -int, int},{z2, -int, int}, Method->MonteCarlo, MaxPoints->10^7,PrecisionGoal->2]

E0 := -0.5

E1 := 2*E0 + 2*Cb + EWW+1/RAB

AE := -NIntegrate[ψB1*ψA1*1/R1B, {x1, -int, int},{y1, -int, int},{z1, -int, int},PrecisionGoal->2]

EAW:= NIntegrate[ψB1*ψA2*ψA1*ψB2*1/R12, {x1, -int, int},{y1, -int, int},{z1, -int, int}, {x2, -int, int},{y2, -int, int},{z2, -int, int}, Method->MonteCarlo, MaxPoints->10^7, PrecisionGoal->2]

E2:=2*E0*S^2+2*AE*S+EAW+S^2*1/RAB

ETOT := 2*(E1+E2)/norm*27.2114

Table[{RAB, ETOT},{RAB, 4,6,1}]