(* Li2 with external clouds, kim_loebl_param, with polarization 13.08.2012 *)
Clear[u1,u2,u3,u4, c0, c1,c2,c3,T,Vnn,Vne,Vee]; Z = 3.0;
(* parameters with 6-311G *)
c3={k1-> 1.09743496, k2-> 0.70863834,s1-> 0.20592185,s2-> 0.19895393};

Vne=-2*Z(3-(R/Q)^2)/Q-4*Z/(P+Q+R)-4*Z/(2*P+2*Q+R);
Vee=2*4/(P+Q)+2/(P+Q)+3*s2/P+2*3*s1/Q;
Vnn=Z*Z/(P+Q+R)/2;
T=2.25*k2/P^2+2*2.25*k1/Q^2;

func=T+Vne+Vee+Vnn /. c3
N[t=FindMinimum[func, {P,1.96822},{Q,0.58373},{R,-0.0272199}],10]
u4=t[[2]];
N[-(Vne+Vee+Vnn)/T /.u4 /.c3,10]
N[2*(P+Q+R)*0.529177 /.u4/.c3,10]
Vne /.u4 /.c3
Vee /.u4 /.c3
Vnn /.u4 /.c3
force = 2*R*Z/Q^3+2*Z/(P+Q+R)^2+2*Z/(2*P+2*Q+R)^2-Z^2/(2*(P+Q+R))^2 /.u4 /.c3

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