(* CH4 molecule with C[1s2 4sp], overlap, parametrized ES 30.12.2012 *)
(* Variables: P,Q,R,S *)
(* with exchange *)

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

Z=6.0; w3=Sqrt[3]; wd=Sin[ArcCos[-1/3]];   ad=Sqrt[8/3];                                               
R=1.046; S= 1.012; 
(* core *)
Tc=2.25*k1/P^2;
Vnec=-3.0*Z/P;
Veec=3.0*s1/P;

(* 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=(2/Q)*(3.0-(R/Q)^2-0.6*(P/Q)^2);
Vex  = -1.2*(P/Q)^3/Q;
Veecb = Veecb + Vex;
Vnncb=Z/(R+S);

(* bond-bond *)
d1= (R+S)*ad;
d2= Sqrt[R^2+(R+S)^2+(2/3)*R*(R+S)];
d3= R*ad;
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+4*Tb;
Vne = Vnec+4*Vneb+4*Vnecb+6*Vnebb;
Vee = Veec+4*Veeb+4*Veecb+6*Veebb;
Vnn = 4*Vnncb+6*Vnnbb;

(* c={k1 -> 1.0, k2-> 0.625, s1 -> 0.3, s2 -> 0.3};  *)c={k1 -> 1.0525, k2 -> 0.6, s1 -> 0.25742249,s2 -> 0.22958899};

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

t=FindMinimum[func,{P,0.27},{Q,1.2},{Method -> "Newton"},{MaxIterations -> 500}]

u=t[[2]];

Vne /.u

Vee /.u  /. c

4*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[{0,0},Q],
	Circle[{-wd*R/2,wd*w3*R/2},Q],Circle[{-wd*R/2,-wd*w3*R/2},Q],Circle[{wd*R,0},Q],Disk[{-wd*(R+S)/2,wd*w3*(R+S)/2},0.1],Disk[{-wd*(R+S)/2,-wd*w3*(R+S)/2},0.1],Disk[{wd*(R+S),0},0.1],
	{Thickness[0.012],Line[{{0,0},{wd*(R+S),0}}],Line[{{0,0},{-wd*(R+S)/2,wd*w3*(R+S)/2}}],Line[{{0,0},{-wd*(R+S)/2,-wd*w3*(R+S)/2}}]}}]/. u;
Show[plot1, {AspectRatio->Automatic, PlotRange->{{-3,3},{-3,3}},Axes->True,GridLines->Automatic,Frame->True}]

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

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