(* NH3 with N[1s2 3sp, 1 sp LP], ~no overlap, parametrized ES 12.11.2012
    Variables : 
        P N - core radius
        Q N - H cloud radius
        R lone pair radius
        S distance proton to center of NH cloud
        \[Alpha] angle lonePair-N-H
        \[Beta] angle H-N-H, 
        not a variable : cos\[Beta] = 1.0-1.5*Sin[\[Alpha]]^2 *)
Clear[wc];
Z = 7.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)*Sin[\[Alpha]]*w3;
(* d3 H - H *)
d3 = (P + Q + S)*Sin[\[Alpha]]*w3;
(* d5 w(NH) - H, \[Beta] is angle HNH *)
cos\[Beta] = 1.0 - 1.5*Sin[\[Alpha]]^2;
d5 = Sqrt[(P + Q + S)^2 + (P + Q)^2 - 2*(P + Q + S)*(P + Q)*cos\[Beta]];
Veebb = 4.0/d1;
Vnebb = -4.0/d5;
Vnnbb = 1.0/d3;

(*bond - lone pair*)
(* d2 w(NH) - w(LP) *)
d2 = Sqrt[(P + Q)^2 + (P + R)^2 - 2*(P + Q)*(P + R)*Cos[\[Alpha]]];
Veebp = 4/d2;
(* d4 w(LP) - H *)
d4 = Sqrt[(P + Q + S)^2 + (P + R)^2 - 2*(P + Q + S)*(P + R)*Cos[\[Alpha]]];
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, 
  or LP-O-H *)
(* HF is identical to NH3 with LP/bond reversed, one angle *)

(*sum terms*)
T   = Tc   + 3*Tb   + Tp;
Vne = Vnec + 3*Vneb +        3*Vnecb + Vnecp + 3*Vnebb + 3*Vnebp;
Vee = Veec + 3*Veeb + Veep + 3*Veecb + Veecp + 3*Veebb + 3*Veebp;
Vnn =                        3*Vnncb +         3*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.95894892, k2 -> 1.34927913, k3 -> 1.23272312, s1 ->
       0.42320792, s2 -> 0.29390809, s3 -> 0.39204779};

func = T + Vne + Vee + Vnn /. par;
t = FindMinimum[func, {P, 0.2125}, {Q, 1.12}, {R, 1.355}, {S, 
    0.58}, {\[Alpha], 1.97}, MaxIterations -> 500]

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/\[Pi]*ArcCos[cos\[Beta]] /. u
(* angle LP - N - H *)
180/\[Pi]*\[Alpha] /. u
d1 /. u
d3 /. u
d5 /. u
d2 /. u
d4 /. u

wd = Sin[\[Alpha]] /. u;  wc = Cos[\[Alpha]] /. 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}, Q], 
        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}]

(* 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.05}], 
            Line[{{0, 0}, {0, P + R}}]} }}] /. u
    
Show[plot2, {AspectRatio -> Automatic, PlotRange -> {{-3, 3}, {-3, 3}}, 
    Axes -> True, GridLines -> Automatic, Frame -> True}]

(******************** RESULTS ***********************)
{-56.5069, {P -> 0.212559, Q -> 1.12155, R -> 1.3546, 
    S -> 0.579361, \[Alpha] -> 1.97065}}
-156.312
31.3404
11.9575
-2.
1.01256
105.824
112.91
2.12847
3.05279
2.61408
2.4215
2.90725