(* CH4   28.05.2013; all interactions with matrix operations; e.g. -4*(3-(p/R2)^2)/R2 for H atoms
 is included in dwnr  *)

Clear[z,sig,k,nc,R1,R2,p,Ekin,Epot,vne,vnn,vee,ow,cn,n,w,c,c2,ec,dw,dn,dwn,dwnr,dwr,dnr,d1,d2,pi,i,j,t,λ];

(*R1=0.26537992;R2=1.2343054;p=0.5334356688274884;*)

z=6.0;
nc=5; bohr=0.529177;sq=Sqrt[3];
k={1.04605075,1.20889711,1.20889711,1.20889711,1.20889711}; 
sig={0.2650207,0.30451569,0.30451569,0.30451569,0.30451569};

dw=Table[0,{i,1,nc},{j,1,nc}];
dn=dw; dwn=dw; dwr=dw; dnr=dw;dwnr=dw;
c2=Table[0,{i,1,nc}];

ow={-2,-2,-2,-2,-2};
cn={z,1,1,1,1};

c = {R1, R2, R2, R2, R2}; 
ec = {0, p, p, p, p}; 
d1 = (R1 + R2 + p)/sq; d2 = (R1 + R2)/sq; 
n = {{0, 0, 0}, {d1, d1, d1}, {-d1, -d1, d1}, {d1, -d1, -d1}, {-d1, d1, -d1}}; 
w = {{0, 0, 0}, {d2, d2, d2}, {-d2, -d2, d2}, {d2, -d2, -d2}, {-d2, d2, -d2}}; 
For[i = 1, i <= nc, i++, c2[[i]] = k[[i]]/c[[i]]^2;
For[j = 1, j <= nc, j++, 
    dw[[i,j]] = Norm[w[[i]] - w[[j]]]; dn[[i,j]] = Norm[n[[i]] - n[[j]]]; 
    dwn[[i,j]] = Norm[w[[i]] - n[[j]]]; dwr[[i,j]] = If[i != j, 1/dw[[i,j]], 
        1.5*(sig[[i]]/c[[i]])]; dnr[[i,j]] = If[i != j, 1/dn[[i,j]], 0]; 
    dwnr[[i,j]] = If[i != j, 1/dwn[[i,j]], (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;

t =FindMinimum[func,{R1,0.26},{R2,1.24},{p,0.53},{Method->Automatic,MaxIterations->500}]

Ekin    /. t[[2]]  
vne     /. t[[2]]  
vee     /. t[[2]]  
vnn     /. t[[2]]  
vir=-Epot/Ekin     /. t[[2]]  
d1*sq*bohr     /. t[[2]]  
d1/. t[[2]]

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

wd=w.RotationMatrix[π/4,{0,1,0}];
nd=n.RotationMatrix[π/4,{0,1,0}];
plot2=Table[Graphics[{Circle[{wd[[i,1]],wd[[i,2]]},c[[i]]],Disk[{nd[[i,1]],nd[[i,2]]},0.08]} ],{i,1,nc}]   /. t[[2]]  ;  
Show[plot2,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-3,3},{-3,3}}, Frame -> True}]