(* CH4   25.05.2013; all interactions with matrix operations; e.g. -4*(3-(p/R2)^2)/R2 for H atoms is included in dwnr  *)
(* does not function with FindMinimum! however, with given R1,R2,p from old calc, result ok *)

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; 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};

(* cloud charge, nuclear charge, radii,excentricity *)

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

c={R1,R2,R2,R2,R2};
ec={0,p,p,p,p};
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}};
w={{0,0,0},{d4,d4,d4},{-d4,-d4,d4},{d4,-d4,-d4},{-d4,d4,-d4}};
dw=Table[Norm[w[[i]]-w[[j]]],{i,1,nc},{j,1,nc}];
dn =Table[Norm[n[[i]]-n[[j]]],{i,1,nc},{j,1,nc}];
dwn=Table[Norm[w[[i]]-n[[j]]],{i,1,nc},{j,1,nc}];
dwr=Table[If[i≠j,1/dw[[i,j]],1.5*sig[[i]]/c[[i]]],{i,1,nc},{j,1,nc}];
dnr=Table[If[dn[[i,j]]> 0,1/dn[[i,j]],0],{i,1,nc},{j,1,nc}];
dwnr=Table[If[i≠j,1/dwn[[i,j]],(1.5-0.5*(ec[[i]]/c[[i]])^2)/c[[i]]],{i,1,nc},{j,1,nc}];
c2=Table[k[[i]]/c[[i]]^2,{i,1,nc}];
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.26538},{R2,1.2343},{p,0.5334},{Method -> "Newton"}, {MaxIterations -> 500}]

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[{w[[i]]},c[[i]]],{Darker[Red,1.0],Sphere[{n[[i]]},0.08]}} ],{i,1,nc}]   /. t[[2]]  ;
Show[plot1,{AspectRatio -> Automatic,Axes -> True}]

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