C2H6 Molecule after Kimball, non overlapping, parameterized with MP2/6-31G(d)

Clear[z,σ,k,nc,R1,R2,p,Ekin,Epot,vne,vnn,vee,vir,ow,cn,n,w,c,c2,ec,
dwnr,dwr,dnr,d,d1,d2,pi,i,j,t,λ,cs,ss,wi];
bohr=0.529177; 
 (*Parameterlist c *)
(*c={k1->1.02669238,k2->1.36467050,k3->1.20088548,s1->0.30102681,s2->0.34863843,s3->0.34858036};*)
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; n=2;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={0,0,0,p,p,p,p,p,p};
d1=R1+R2; d2=d1+(R1+R3)*cs; d3=(R1+R3)*ss; d4= d3/2; d5= d3*Sqrt[3]/2;
d6=d1+(R1+R3+p)*cs; d7=(R1+R3+p)*ss; d8=d7/2; d9=d7*Sqrt[3]/2;
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 = {{d1, 0, 0}, {-d1, 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},{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*d1  /.u;
dcc=cc*bohr

(*Virial Theorem *)
vir/.u

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

plot1=Table[Graphics3D[{Opacity[0.6],Sphere[{w[[i]]},c[[i]]],{Darker[Red,1.0],If[i≠3,Sphere[{n[[i]]},0.08]]}} ],{i,1,nc}]   /. u  ;
Show[plot1,{AspectRatio -> Automatic,Axes -> True}]