(* LiH on x-axis; with polarization, Li-cloud at X=0 *)
(* 1 Li (x from center of P), 2 H (Q from center of R), 3 Li-cloud (P),  4 H-cloud (R), Z=3 inserted *)
c2={k1-> 1.02480344, k2-> 1.232443643,s1-> 0.27877144,s2-> 0.28887416};
T = k2*2.25/R^2 +k1*2.25/P^2;
d12 =P+R+Q-x;
d13 = x;
d14 = P+R-x; 
d23 = P+R+Q;
d24 = Q;
d34 = P+R;

Vee = 3*s2/R+3*s1/P+4/d34;
Vnn = 3/d12;
Vne = -3*(3-(d13/P)^2)/P-6/d14-(3-(d24/R)^2)/R-2/d23;

func = T + Vee + Vnn + Vne /. c2
(* arbitrary center positions on x-axis, 4 variables *) 
N[t = FindMinimum[func, {R,{1.8,2.2}},{P,{0.5,0.6}},{Q,{0.4,0.5}},
{x,{0.01,0.02}},{MaxIterations -> 100}],20]  
u = t[[2]];
(* lengths *)
d12 /.u
d13 /. u
d14 /. u
d23 /. u
d24 /. u
d34 /.u
N[(Vnn+Vee+Vne)/T /. u /. c2, 10]  (* Virial *)
N[-6*x/P^3-3/d12^2+6/d14^2 /. u , 10]    (* Forces on Li3+ *)
N[2*Q/R^3+2/d23^2-3/d12^2 /. u, 10](* Forces on H+ *)

plot1 = Graphics[{Circle[{0,0},P],Disk[{x,0},0.03],Circle[{P + R,0},R],Disk[{P +R+ Q,0},0.03],{Dashing[{0.03,0.03}],Circle[{P+R,0},Q],
{Circle[{0,0},x]}}}]
/.u

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