C2H6, non overlapping, polarized, parameterized with MP2/6-31G(d) ES 04 June 2013

Clear[z,σ,k,nc,R1,R2,R3,p,Ekin,Epot,vne,vnn,vee,vir,ow,cn,wi,c2,ec,
dwnr,dwr,dnr,d,d1,d2,d3,d4,d5,d6,d7,d8,d9,nn,n,q,i,j,t,λ,cs,ss,wi];
bohr=0.529177;  rad=57.2957795;
 (*Parameterlists *)
k={1.02669238,1.02669238,1.36467050,1.20088548,1.20088548,1.20088548,1.20088548,1.20088548,1.20088548}; 
σ={0.30102681,0.30102681,0.34863843,0.34858036,0.34858036,0.34858036,0.34858036,0.34858036,0.34858036};

(*Terms analogous to CH4_mat *)
z=6.0; nc=9; nn=2;(*wi=π-ArcCos[-1/3];*) cs=Cos[wi]; ss=Sin[wi];
dwr=Table[0,{i,1,nc},{j,1,nc}];
dnr=dwr;dwnr=dwr;
c2=Table[0,{i,1,nc}];

ow={-2,-2,-2,-2,-2,-2,-2,-2,-2};
cn={z,z,0,1,1,1,1,1,1};
c={R1,R1,R2,R3,R3,R3,R3,R3,R3};
ec={q,q,0,p,p,p,p,p,p};
d0=R1+R2+q; d1=R1+R2; 
d2=d1+(R1+R3)*cs;         d3=(R1+R3)*ss;        d4= d3/2; d5= d4*Sqrt[3];
d6=d0+(R1+R3+p)*cs; d7=(R1+R3+p)*ss; d8=d7/2; d9=d8*Sqrt[3];
w = {{d1, 0, 0}, {-d1, 0,0}, {0,0,0},
   {d2, d3, 0},    {d2, -d4, d5}, {d2, -d4, -d5},
{-d2,-d3,0},{-d2,d4,-d5},   {-d2,d4,d5}}; 
n = {{d0, 0, 0}, {-d0, 0,0}, {0,0,0},
   {d6, d7, 0},   {d6, -d8, d9}, {d6, -d8, -d9},
{-d6,-d7,0},{-d6,d8,-d9},  {-d6,d8,d9}}; 
For[i = 1, i ≤ nc, i++,
c2[[i]] = k[[i]]/c[[i]]^2;
For[j = 1, j ≤ nc, j++,
d = Norm[w[[i]] - w[[j]]];
dwr[[i,j]] = If[i ≠ j, 1/d,1.5*(σ[[i]]/c[[i]])];
d= Norm[n[[i]] - n[[j]]];
dnr[[i,j]] = If[i ≠ j, 1/d, 0];
d= Norm[w[[i]] - n[[j]]];
dwnr[[i,j]] = If[i ≠ j, 1/d,(1.5 - 0.5*(ec[[i]]/c[[i]])^2)/c[[i]]]]]; 
Ekin = -(9/8)*ow . c2; 
vee = 0.5*ow . (ow . dwr); 
vnn = 0.5*cn . (cn . dnr); 
vne = cn . (ow . dwnr); 
Epot = vne + vee + vnn; 
func = Ekin + Epot;

vir=-Epot/Ekin;

func =Ekin+Epot;

t = NMinimize[{func},{{R1, 0.25, 0.27}, {R2, 1.19, 1.22}, 
   {R3, 1.24, 1.28}, {p, 0.53,0.56},{q, 0.0001,0.02},{wi,1.2,1.25}}]

u = t[[2]];

(*Total Energy *)
Etot=t[[1]]

(*Kinetic energy*)
Ekin /.u

(*NUCLEUS-ELECTRON POTENTIAL ENERGY  *)
vne /.u

(*ELECTRON-ELECTRON POTENTIAL ENERGY *)
vee /.u

(*NUCLEUS-NUCLEUS POTENTIAL ENERGY *)
vnn /.u

(*C-H distance  Å*)
ch =R1+R3+p/.u;
dch = ch*bohr

(*C-C distance  Å *)
cc =  2*d0 /.u;
dcc=cc*bohr

(* angle CCH *)
(π-wi)*rad /. u

(*Virial Theorem *)
vir/.u

(*Atomization Energy at 0 K in kcal/mol *)
(-Etot - nn*37.78899  -( 2*nn + 2)*0.50177)*627.5096  /.u

plot0=Graphics3D[{Opacity[0.6],{Darker[Orange,0.5],Specularity[1,20],Sphere[{w[[1]]},c[[1]]],Sphere[{w[[2]]},c[[2]]]},{Lighter[Orange,0.4],Specularity[1,20],Sphere[{w[[3]]},c[[3]]]},{Darker[Red,1.0],Specularity[1,20],Sphere[{n[[1]]},0.08],Sphere[{n[[2]]},0.08]}}] /.u;
plot1=Table[Graphics3D[{Opacity[0.6],Specularity[1,20],Sphere[{w[[i]]},c[[i]]],{Darker[Red,1.0],Specularity[1,20],Sphere[{n[[i]]},0.08]}}],{i,4,nc}]   /. u  ;
Show[plot0,plot1,{AspectRatio -> Automatic,Axes -> True,SphericalRegion->True}]