(* H2O with O[1s2 2sp, 2sp LP],~no overlap, parametrized ES 04.01.2013
Variables:
       P  O-core radius
Q  O-H cloud radius
R  lone pair radius
S  distance proton to center of OH cloud
α  angle lonePair-O-lonePair
β  angle H-O-H *)

Clear[wc,k1,k2,k3,s1,s2,s3,par];
Z=8.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-bond*)
(* d1 w(NH)-w(NH) *)
d1=(P+Q)*Sqrt[2*(1-Cos[β])];
Veebb=4.0/d1;
(* d3 H-H *)
d3=(P+Q+S)*Sqrt[2*(1-Cos[β])];
Vnnbb=1.0/d3;
(* d5 w(OH)-H, β is angle HOH *)
d5=Sqrt[(P+Q+S)^2+(P+Q)^2-2*(P+Q+S)*(P+Q)*Cos[β]];
Vnebb=-2.0/d5;

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

(*lone pair-lone pair*)
(* none for NH3, but for H2O and HF *)
(* H2O has two independent angles as variables: H-O-H and LP-O-LP, *)
d6 = (P+R)*Sqrt[2*(1-Cos[α])];
Veepp=4/d6;

(*sum terms*)

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

(*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.92438558,k2->1.42313501,k3->1.32139847,s1->0.43732345,s2->0.31779336,s3->0.44741336};

(*par = {k1->0.92438558,k2->1.42313501,k3->1.32139847,s1->0.43732345,s2->0.31779336,s3->0.44741336};
par*)

func=T+Vne+Vee+Vnn/.par;
t=FindMinimum[func,{P,0.18},{Q,1.0},{R, 1.1},{S,0.6}, {α,1.97},{β,1.97},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 HNH *)
180/π*β/.u
(* angle LP-N-LP *)
180/π*α /. u
d1 /.u
d3 /.u
d5 /.u
d2 /.u
d4 /.u

wd = Sin[α] /. u;  wc=Cos[α] /. u;
(* Projection onto HHH-plane *)

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

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

(* Projection onto LP-HH-plane *)

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

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