(* CH4 molecule with C[1s2 4sp], no overlap, parametrized ES 06.11.2012 *)
(* Variables: P,Q,R *)
Clear[k1,k2,s1,s2, z, P,Q,R,Vee,Vne,Vnn,T,d1,d2,d3,par,w3,ad];
Clear[Vnec,Veec,Vneb,Veeb,Veecb,Vnecb,Vnncb,Vnebb,Veebb,Vnnbb];
Clear[Tc,Tb,wd];

z=6.0; wd=0.94281; ad=Sqrt[8/3]; w3=Sqrt[3];

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

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

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

(* bond-bond *)
d1=(P+Q+R)*ad;
d2=Sqrt[(P+Q)^2+(P+Q+R)^2+(2/3)*(P+Q)*(P+Q+R)];
d3=(P+Q)*ad;
Vnebb=-4.0/d2;
Veebb=4.0/d3;
Vnnbb=1.0/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;

(* par={k1 -> 1.0, k2-> 1.0, s1 -> 0.4, s2 -> 0.4}; *) 
(* par={k1 -> 1.01613950, k2 -> 0.65858068,s1-> 0.25742249
s2-> 0.22958899};*)
par={k1->1.05518276,k2-> 1.21738169,s1-> 0.24406300,s2-> 0.30451569};

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

t=FindMinimum[func,{P,0.26},{Q,1.11},{R,0.476}]

u=t[[2]];

Vne /.u

Vee /.u  /. par

Vnn /. u

-(Vne+Vee+Vnn)/T /. u /.par

(P+Q+R)*0.529177 /. u

plot1=Graphics[{Circle[{0,0},P],Disk[{0,0},0.1],Circle[{0,0},Q],
	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+R)/2,wd*w3*(P+Q+R)/2},0.1],Disk[{-wd*(P+Q+R)/2,-wd*w3*(P+Q+R)/2},0.1],Disk[{wd*(P+Q+R),0},0.1],
	{Thickness[0.012],Line[{{0,0},{wd*(P+Q+R),0}}],Line[{{0,0},{-wd*(P+Q+R)/2,wd*w3*(P+Q+R)/2}}],Line[{{0,0},{-wd*(P+Q+R)/2,-wd*w3*(P+Q+R)/2}}]}}]/. u

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