(* C2H6 Molecule after Kimball, non overlapping,parameterized with G3MP2//6-31g(d) *)

Clear[k1,k2,k3,s1,s2,s3];
 n=2.; cw=1/3; sw=Sqrt[2]/3; (*cos and sin of 180-109.5 *)
(* Parameterlist c *)
c={k1->1.01684982,k2->1.36707348,k3-> 1.19949452,s1-> 0.33365014,s2-> 0.33365040,s3-> 0.33448109};
(*c = {k1 ->1.01682967(*4905*),k2 ->1.36655617(*37426345*),k3 -> 1.19993090,s1 -> 0.33353563,s2 -> 0.33353647(*30677632*),s3 -> 0.33459553};*)
(*c = {k1 -> 1.02246687,k2 -> 1.37426345,k3 -> 1.20537762,s1 -> 0.30582536,s2 -> 0.30677632,s3 -> 0.35441063};*)
(* Terms from ChemEdu/Kimball/Alkan.sys *)

T=2.25*n*k1/P^2 +(n - 1.0)*2.25*k2/Q^2 + 4.5*(n + 1.0)*k3/R^2 /.c;
p =X; p2= p*p; z= 6.; ad= 8./3.; zd= 2./ 3.; vd =4./ 3.; q= P + Q; r= P + R; q2= q*q; r2= r*r; r2p2=r2*p2; r2p= r2*p; qr= q*r; pr=p*r;qrp=qr*p;

Vee=3.0*n*s1/P+(n - 1)*3.0*s2/ Q+(n + 1)*6.0*s3/ R + 8.0/ q+ 24.0/ Sqrt[zd*qr+q2 + r2]+ 24.0/ r + 4.0/(2*q) + 24.0/ Sqrt[vd*qr + 4.0*q2 + r2]+ 24.0/ Sqrt[ad*r2] + 12.0/ Sqrt[ad*qr + 4.0*q2 +4.0*r2]+ 24.0/ Sqrt[ad*qr + 4.0*q2 + vd*r2] /.c;

Vne =-2.0*(n + 1.0)*(3.0 -((p - 1.0)*(1.0 + P/ R))^2)/R - 3.0*n*z/ P- 4.0*z/ q - 12.0/ Sqrt[zd*qrp + r2p2 + q2]- 12.0*z/ Sqrt[r2] - 12.0/ pr - 4.0*z/ (2*q)- 12.0*z/ Sqrt[vd*qr + 4.0*q2 + r2] - 12.0/ Sqrt[vd*qrp + r2p2 + 4.0*q2]- 24.0/ Sqrt[zd*r2p + r2p2 + r2]-12.0/ Sqrt[vd*qr + vd*qrp + 2.0*r2p + r2p2 + 4.0*q2 + r2]- 24.0/ Sqrt[vd*qr + vd*qrp - zd*r2p + r2p2 + 4.0*q2 + r2] ;

Vnn =6.0*z/ (r*p) + z*z/ (2*q)+ 6.0*z/ Sqrt[vd*qrp + r2p2 + 4.0*q2] + 6.0/ Sqrt[ad*r2p2]+ 3.0/ Sqrt[ad*qrp + 4.0*r2p2 + 4.0*q2]+ 6.0/ Sqrt[ad*qrp + vd*r2p2 + 4.0*q2] ;

func = T + Vne + Vee + Vnn;
t = FindMinimum[func,{P,0.25,0.27},{Q,1.1,1.2},{R,1.2,1.3},{X,1.3,1.4},{Method->Automatic,MaxIterations->500}]

u = t[[2]];
T /.u
Vne /.u
Vee /.u /.c
Vnn /.u
ch = X*(P + R) /.u;
dch = ch*0.529177                                             (* C-H distance *)
cc = 2.*(P + Q) /.u;
dcc=cc*0.529177                                               (* C-C distance *)
-(Vne + Vee + Vnn)/T /.u /.c              (* Virial Theorem *)
(*atomization energy at 0 K in kcal/mol *)(-t[[1]] - n*37.78899 -( 2*n + 2)*0.50177)*627.50956

plot1 = Graphics[{Circle[{-q,0},P],Disk[{-q,0},0.1],
Circle[{q,0},P],Disk[{q,0},0.1],
Circle[{0,0},Q],
Circle[{-q-r*cw,-r*2*sw},R],Disk[{-q- pr*cw,-pr*2*sw},0.1],Circle[{q +r*cw,r*2*sw},R],Disk[{q+ pr*cw,pr*2*sw},0.1],
{Dashing[{0.02,0.02}],
Circle[{q+r*cw,-r*sw},R],Disk[{q + pr*cw,-pr*sw},0.1],Circle[{-q-r*cw,r*sw},R],Disk[{-q-pr*cw,pr*sw},0.1]},{Thickness[0.01],Line[{{-q,0},{q,0}}],Line[{{q,0},{q + pr*cw,pr*2*sw}}],Line[{{-q,0},{-q- pr/3.,-pr*2*sw}}],
{Dashing[{0.01,0.01}],
Line[{{q,0},{q + pr*cw,-pr*sw}}],
Line[{{-q,0},{-q - pr*cw,pr*sw}}]}}}] /.u;

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