(* CH4   24.05.2013; all with matrix operations; e.g. -4*(3-(p/R2)^2)/R2 for H atoms is included in cn.(ow.dwnr), see below *)

Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,n,c,dw,dn,dwn,dwnr,dwr,dnr,a1,a2,a3,b1,b2,b3,d1,d2,d3,d4,d5,pi,i,j,d,t];

R1=0.26537992;R2=1.2343054;p=0.5334356688274884;

z=6.0;  vne=0.0;
nc=5; pi=0.0; bohr=0.529177;sq=Sqrt[3];
(*k=Table[1.0,{nc}]*)
k={1.04605075,1.20889711,1.20889711,1.20889711,1.20889711}; 
(*sig=Table[0.4,{nc}]*)
sig={0.2650207,0.30451569,0.30451569,0.30451569,0.30451569};

(* C He shell, i=j *)

ow={-2,-2,-2,-2,-2};
cn={z,1,1,1,1};
c={R1,R2,R2,R2,R2};
Ekin = -(9/8)*Sum[k[[i]]*ow[[i]]/c[[i]]^2,{i,1,nc}]
vee   = -(3/2)*Sum[sig[[i]]*ow[[i]]/c[[i]],{i,1,nc}]

d2=(R1+R2+p)/sq;  d4=(R1+R2)/sq;

n={{0,0,0},{d2,d2,d2},{-d2,-d2,d2},{d2,-d2,-d2},{-d2,d2,-d2}};
MatrixForm[n]

w={{0,0,0},{d4,d4,d4},{-d4,-d4,d4},{d4,-d4,-d4},{-d4,d4,-d4}};
MatrixForm[w]

dw=Table[Norm[w[[i]]-w[[j]]],{i,1,nc},{j,1,nc}];
MatrixForm[dw]

dn =Table[Norm[n[[i]]-n[[j]]],{i,1,nc},{j,1,nc}];
MatrixForm[dn]

dwn=Table[Norm[w[[i]]-n[[j]]],{i,1,nc},{j,1,nc}];
(*dwn[[1,1]]=1/1.5*R1;*)
MatrixForm[dwn]

dwr=Table[If[dw[[i,j]]> 0,1/dw[[i,j]],0],{i,1,nc},{j,1,nc}];
MatrixForm[dwr]

dnr=Table[If[dn[[i,j]]> 0,1/dn[[i,j]],0],{i,1,nc},{j,1,nc}];
MatrixForm[dnr]

dwnr=Table[If[dwn[[i,j]]> 0,1/dwn[[i,j]],1.5/c[[i]]],{i,1,nc},{j,1,nc}];
MatrixForm[dwnr]

vee=vee+0.5*ow.(ow.dwr);
vnn=0.5*cn.(cn.dnr);
vne=(*vne+*)cn.(ow.dwnr);  (*0.90069**)
MatrixForm[ow.dwnr]
cn.%

b1={0,d2,-d2,d2,-d2};
b2={0,d2,-d2,-d2,d2};
b3={0,d2,d2,-d2,-d2};        
a1={0,d4,-d4,d4,-d4};
a2={0,d4,-d4,-d4,d4};
a3={0,d4,d4,-d4,-d4};

Epot=vne+vee+vnn;
func=Ekin+Epot;
(*t = FindMinimum[func,{R1,0.267},{R2,1.244},{p,0.65},{Method -> Automatic}, {MaxIterations -> 500}]*)

func
Ekin
vne  (*/. t[[2]]*)
vee  (*/. t[[2]]*)
vnn  (*/. t[[2]]*)
vir=-Epot/Ekin  (*/. t[[2]]*)
(R1+R2+p)*bohr  (*/. t[[2]]*)

plot1=Table[Graphics3D[{Opacity[0.6],Sphere[{a1[[i]],a2[[i]],a3[[i]]},c[[i]]],{Darker[Red,1.0],Sphere[{b1[[i]],b2[[i]],b3[[i]]},0.08]}} ],{i,1,nc}](*/. t[[2]]*);

Show[plot1,{AspectRatio -> Automatic,Axes -> True}]

plot2=Table[Graphics[{Circle[{a1[[i]],a3[[i]]},c[[i]]],Disk[{b1[[i]],b3[[i]]},0.08]} ],{i,1,nc}](*/. t[[2]]*);(*/. t[[2]]*);

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

(* from old program *)       
Ekin=2.25*k[[1]]/R1^2+9.0*k[[2]]/R2^2;

vne=-3*z/R1

vee=3*sig[[1]]/R1;
       p=p/(R1+R2)+1; wd=sq; f=Sqrt[2/3];
       p2=p*p; r=R1+R2; r2=r*r;

vne=vne-4*(3-((p-1)*(1+R1/R2))^2)/R2

vee=vee+4*3*sig[[2]]/R2;
       vee=vee+8*2/r;

vne=vne-8*z/r

vne=vne-8/(r*p)

vnn=4*z/(r*p);
       vee=vee+12/(f*r);

vne=vne-24*wd/Sqrt[2*r2*p + 3*r2*p2 + 3*r2]

vnn=vnn+3/(f*r*p);



Ekin
Epot=vne+vnn+vee
Epot+Ekin
vne  (*/. t[[2]]*)
vee  (*/. t[[2]]*)
vnn  (*/. t[[2]]*)
vir=-Epot/Ekin  (*/. t[[2]]*)
(R1+R2)*p*bohr  (*/. t[[2]]*)