(* BH3 molecule with B[1s2 3sp], overlap, parametrized ES 31.12.2012 *)
(* Variables: P,Q,R,S *)
(* with core exchange *)

Clear[k1,k2,s1,s2, Z, P,Q,R,S,Vne,Vee,Vnn,Veec,Vnec,Veecb,Vnecb,Vnncb,Veebb,Vnebb,Vnebb,d1,d2,d3];

Z=5.0; w3=Sqrt[3];
R=1.11; S=1.14;
(* B core *)
Tc=2.25*k1/P^2;
Vnec=-3.0*Z/P;
Veec=3.0*s1/P;

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

(* core-bond *)
Vnecb=-(2.0*Z/Q)*(1.5-0.5*(R/Q)^2)-2.0/(R+S);
Veecb=(4/Q)*(1.5-0.5*(R/Q)^2-0.3*(P/Q)^2);
Vex  = -1.2*(P/Q)^3/Q;
Veecb = Veecb + Vex;
Vnncb=Z/(R+S);

(* bond-bond *)
d1=(R+S)*w3;
d2=Sqrt[R^2+(R+S)^2+R*(R+S)];
d3=R*w3;
Vnebb = -4/d2;
Veebb = (4/Q)*(1.2-0.5*(d3/Q)^2+(3/16)*(d3/Q)^3-(1/160)*(d3/Q)^5);
Vnnbb =  1/d1;

(* sum terms *)
T = Tc+3*Tb;
Vne = Vnec+3*Vneb+3*Vnecb+3*Vnebb;
Vee = Veec+3*Veeb+3*Veecb+3*Veebb;
Vnn = 3*Vnncb+3*Vnnbb;

(* Parameters in list c *)

(* c={k1 -> 1.0, k2-> 1.0, s1 -> 0.3, s2 -> 0.3}; *)  
(*c={k1 -> 1.01613950, k2 -> 0.701, s1 -> 0.25742249,s2 -> 0.22958899};*)
(*c={k1 -> 1.04011439, k2 -> 0.75422468, s1 -> 0.21404979,s2 -> 0.07673395};*)
 (*c={k1 -> 1.0, k2 -> 0.75, s1 -> 0.3,s2 -> 0.3};*)
c={k1 -> 1.02, k2 -> 0.683, s1 -> 0.21,s2 -> 0.21};

func=T+Vne+Vee+Vnn /. c;

t=FindMinimum[func,{P,0.32},{Q,1.985},{Method -> "Newton", MaxIterations -> 1000}]

u=t[[2]];

Vne /.u

Vee /.u  /. c

3*Vex /. u

Vnn /. u

-(Vne+Vee+Vnn)/T /.u/.c
(R+S)*0.529177 /.u

plot1=Graphics[{Circle[{0,0},P],Disk[{0,0},0.1],
			Circle[{R,0},Q],Circle[{-R/2,w3*R/2},Q],Circle[{-R/2,-w3*R/2},Q],Disk[{R+S,0},0.1],Disk[{-(R+S)/2,w3*(R+S)/2},0.1],Disk[{-(R+S)/2,-w3*(R+S)/2},0.1],
			{Thickness[0.012],Line[{{0,0},{R+S,0}}],Line[{{0,0},{-(R+S)/2,w3*(R+S)/2}}],Line[{{0,0},{-(R+S)/2,-w3*(R+S)/2}}]}}]/. u;

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

plot2=Graphics3D[
{Opacity[0.3],{Sphere[{R,0,0},Q],Sphere[{-R/2,w3*R/2,0},Q],Sphere[{-R/2,-w3*R/2,0},Q]},{Darker[Red],Sphere[{R+S,0,0},0.1],Sphere[{-(R+S)/2,w3*(R+S)/2,0},0.1],Sphere[{-(R+S)/2,-w3*(R+S)/2,0},0.1],Sphere[{0,0,0},P]},{Darker[Red,1],Sphere[{0,0,0},0.1]}}]/. u;

Show[plot2, {AspectRatio->Automatic,Boxed -> False}]