(* HF with F[1s2 1sp, 3p LP],~no overlap, parametrized ES 05.01.2013
Variables:
       P  F-core radius
Q  F-H cloud radius
R  lone pair radius
S  distance proton to center of FH cloud
α angle lonePair-F-H
β angle LP-F-LP not a variable: cosβ = 1.0-1.5*Sin[α]^2 *)

Clear[wc,k1,k2,k3,s1,s2,s3,par];
Z=9.0; w3=Sqrt[3];

(*core*)
Tc=2.25*k1/P^2;
Vnec=-3.0*Z/P;
Veec=3.0*s1/P;

(*XH bond*)
Tb=2.25*k2/Q^2;
Vneb=-2.0*(1.5-0.5*(S/Q)^2)/Q;
Veeb=3.0*s2/Q;

(*lone pair*)
Tp=2.25*k3/R^2;
Veep=3.0*s3/R;

(*core-bond*)
Vnecb=-2.0*Z/(P+Q)-2/(P+Q+S);
Veecb=4.0/(P+Q);
Vnncb=Z/(P+Q+S);

(*core-lone pair*)
Vnecp=-2.0*Z/(P+R);
Veecp=4.0/(P+R);

(*bond-lone pair*)
(* d2 w(FH)-w(LP) *)
d2=Sqrt[(P+Q)^2+(P+R)^2-2*(P+Q)*(P+R)*Cos[α]];
Veebp = 4/d2;
(* d4 w(LP)-H *)
d4=Sqrt[(P+Q+S)^2+(P+R)^2-2*(P+Q+S)*(P+R)*Cos[α]];
Vnebp =-2.0/d4 ;

(*lone pair-lone pair*)
(* none for NH3, but for H2O and HF *)
(* HF has one angle as variable: H-F-LP  *)
(* d1 w(LP)-w(LP) *)
d1=(P+R)*Sin[α]*w3;
Veepp = 4.0/d1;
(*sum terms*)

T=Tc+           Tb+      3*Tp;
Vne=Vnec+Vneb+                     Vnecb+3*Vnecp+ 3*Vnebp;
Vee=Veec+Veeb+ 3*Veep+Veecb+3*Veecp+ 3*Veebp +3*Veepp;
Vnn=                                                 Vnncb;

(*Parameters: G3//6-311G(d,p) can be read for LiH ... HF
in table from www.kimball-model.org/Kimball/params.xls *)
par={k1->0.91058613,k2->1.41656890,k3->1.21497466,s1->0.51017912,s2->0.44166435,s3->0.36888374};

func=T+Vne+Vee+Vnn/.par;
t=FindMinimum[func,{P,0.16},{Q,1.0},{R, 0.83},{S,0.6}, {α,1.82},MaxIterations->500, Method -> "Newton"]

u=t[[2]];
Vne/.u
Vee/.u/.par
Vnn/.u
(Vne+Vee+Vnn)/T/.u/.par
(P+Q+S)*0.529177/.u
(* angle LP-FH *)
180/π*α /. u
(* angle LP-F-LP *)
180/π*ArcCos[1.0-1.5*Sin[α]^2] /. u
d1 /.u
d2 /.u
d4 /.u

wd = Sin[α] /. u;  wc=Cos[α] /. u;
plot3 = Graphics3D[{Sphere[{0,0,0}, P],Sphere[{0,0,0}, 0.1],{Opacity[0.5], Sphere[{0,0,P+Q}, Q]}, 
Sphere[{-wd*(P + R)/2, wd*w3*(P + R)/2,wc*(P+R)}, R], 
Sphere[{-wd*(P + R)/2, -wd*w3*(P + R)/2,wc*(P+Q)}, R], Sphere[{wd*(P + R), 0,wc*(P+R)}, R],{Darker[Red],
Sphere[{0,0,P+Q+S}, 0.06]}}] /. u;
Show[plot3,{AspectRatio -> Automatic, Boxed -> False}]

plot2=Graphics[{Circle[{0,0},P],Disk[{0,0},0.1],Circle[{0,P+Q},Q],Disk[{0,P+Q+S},0.1],Circle[{wd*(P+R),wc*(P+R)},R],{Dashing[{0.02}],Circle[{-wd*(P+R)/2,wc*(P+R)},R]},{Thickness[0.01],Line[{{0,0},{0,P+Q+S}}],{Dashing[{0.01}],Line[{{0,0},{-wd*(P+R)/2,wc*(P+R)}}]},{Dashing[{0.01}],Line[{{0,0},{P+R,wc*(P+R)}}]} }}]/.u;

Show[plot2,{AspectRatio->Automatic,PlotRange->{{-2.5,2},{-1.5,2.5}},Axes->True,GridLines->Automatic,Frame->True}]

plot1=Graphics[{Circle[{0,0},P],Disk[{0,0},0.1],Circle[{0,0},Q],Circle[{-wd*(P+R)/2,wd*w3*(P+R)/2},R],Circle[{-wd*(P+R)/2,-wd*w3*(P+R)/2},R],Circle[{wd*(P+R),0},R]}]/.u;
Show[plot1,{AspectRatio->Automatic,PlotRange->{{-1.5,2},{-2,2}},Axes->True,GridLines->Automatic,Frame->True}]