(* Li2 with external clouds, Kimball_Loebl orig, 13.0.2012 *)
Clear[u1,u2,u3,u4, c0, c1,c2,c3,T,Vnn,Vne,Vee]; Z = 3.0;
(* Kimball's parameters *)
c0={k1-> 1.0, k2-> 1.0,s1-> 0.4,s2-> 0.4};

Vne=-2*3*Z/Q-4*Z/(P+Q)-2*Z/(P+Q);
Vee=2*4/(P+Q)+2/(P+Q)+3*s2/P+2*3*s1/Q;
Vnn=Z*Z/(P+Q)/2;
T=2.25*k2/P^2+2*2.25*k1/Q^2;
func=T+Vne+Vee+Vnn /. c0
N[t=FindMinimum[func, {P,2.0},{Q,0.6}],10]

u1=t[[2]];
N[-(Vne+Vee+Vnn)/T /.u1 /.c0,10]
N[2*(P+Q)*0.529177 /.u1 /.c0,10]
Vne /.u1 /.c0
Vee /.u1 /.c0
Vnn /.u1 /.c0

plot0=Graphics[{Circle[{-(P+Q),0},Q],Circle[{P+Q,0},Q],Circle[{0,0},P],Disk[{-(P+Q),0},0.08],Disk[{(P+Q),0},0.08],{Thickness[0.01],Line[{{-(P+Q),0},{(P+Q),0}}]}}]/.u1
Show[plot0,{AspectRatio->Automatic,Axes->True,GridLines->Automatic,PlotRange->{{-5,5},{-4,4}},Frame->True}]