(* LiH molecule after Kimball, no overlap,
parametrized with RHF/6 - 311  g *)

Clear[k1];  Clear[k2];  Clear[s1];  Clear[s2];

T = 2.25*k1/P^2 + 2.25*k2/Q^2;

Vne = -9.0/P - (3.0 - R^2/Q^2)/Q - 2.0/(P + Q + R) - 
          2.0*3.0/(P + Q);

Vee = 3.0*s1/P + 3.0*s2/Q + 4.0/(P + Q);

Vnn = 3.0/(P + Q + R);

func = T + Vne + Vee + Vnn  /. c

(* Parameters in list c have been determined by "least \
squares fit" to Etot,  Vne,  Vee,  and Re *)

c = {k1  ->  1.02805521,  k2 ->  1.23258566,  
          s1  ->  0.27438341,  s2  ->  0.28887266};

N[t = FindMinimum[func, {P, 0.56}, {Q, 2.022}, {R, 0.449}], 
      10]

u = t[[2]];

N[Vne  /. u, 10]

N[Vee  /. u  /. c, 9]

N[Vnn  /.  u, 8]

(P + Q + R)*0.529177  /. u

N[-(Vne + Vee + Vnn)/T  /. u /. c, 9]

\(plot1 = 
      Graphics[{Circle[{0, 0}, P], 
            Disk[{0, 0}, 0.1], 
			Circle[{P + Q, 0}, Q], 
            Disk[{P + Q + R, 0}, 0.1], 
			{Thickness[0.012], 
              Line[{{0, 0}, {P + Q + R, 0}}]}}]
\ \ \ \ \ \ \  /. \ 
        u\)

Show[
        plot1,  {AspectRatio -> Automatic,  
          PlotRange -> {{-1, 6}, {-3.5, 3.5}}, Axes -> True, 
          GridLines -> Automatic, Frame -> True}]