(* CH4   15.06.2012 *)

Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,
a1,a2,a3,b1,b2,b3,c1,c2,c3,cs,ss,d1,d2,d3,d4,d5,pi,i,j,d,t];

(*z=6.0; sig1=0.3; sig2=0.3; sig4=0.3; nc=5; pi=0.0;
k1=1.0; k2=1.2; k4=1.66; bohr=0.529177;*)

z=6.0; 
k1=1.04605075;  k2=1.20889711;  (*sig1=0.25273353;  sig2=0.32960771;*)
sig1=0.2650207; sig2=0.30451569;
nc=5; pi=0.0; bohr=0.529177;sq=Sqrt[3];

(* C He shell, i=j each *)

Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
vne=-3.0*z/R1;
(* bonding pairs *)
Ekin = Ekin + 2.25*(4*k2/R2^2);
vee=vee+3.0*(4*sig2/R2);
c1={-2,-2,-2,-2,-2};
c2={z,1,1,1,1};
c3={R1,R2,R2,R2,R2};
d2=(R1+R2+p)/sq;  d4=(R1+R2)/sq;
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};
vne=vne-4*(3-(p/R2)^2)/R2;
(*vne=vne-4*(3-((p-1)*(1+R1/R2))^2)/R2;*)
vnn = 0.0;
(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
vee =vee + c1[[i]]*c1[[j]]/Sqrt[(a1[[i]]-a1[[j]])^2+(a2[[i]]-a2[[j]])^2+(a3[[i]]-a3[[j]])^2]]]

(*nn*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,vnn = vnn + c2[[i]]*c2[[j]]/Sqrt[(b1[[i]]-b1[[j]])^2+(b2[[i]]-b2[[j]])^2+(b3[[i]]-b3[[j]])^2]]]
(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nc+1, j++,d=Sqrt[(a1[[i]]-b1[[j]])^2+(a2[[i]]-b2[[j]])^2+(a3[[i]]-b3[[j]])^2];
    If[d≠0,vne = vne + c1[[i]]*c2[[j]]/d]]]      

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

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

plot1=Graphics[{Circle[{a1[[1]],a2[[1]]},R1],      Circle[{a1[[2]],a2[[2]]},R1],Circle[{a1[[3]],a2[[3]]},R2],Circle[{a1[[4]],a2[[4]]},R2],Circle[{a1[[5]],a2[[5]]},R2],Circle[{a1[[6]],a2[[6]]},R2],Circle[{a1[[7]],a2[[7]]},R4],Circle[{a1[[8]],a2[[8]]},R4],Disk[{b1[[1]],b2[[1]]},0.08], Disk[{b1[[2]],b2[[2]]},0.08], Disk[{b1[[3]],b2[[3]]},0.08], Disk[{b1[[4]],b2[[4]]},0.08],Disk[{b1[[5]],b2[[5]]},0.08],Disk[{b1[[6]],b2[[6]]},0.08]} ] /. t[[2]]

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

plot2=Graphics[{Circle[{a1[[1]],a2[[1]]},R1],      Circle[{a1[[2]],a3[[2]]},R1],Circle[{a1[[3]],a3[[3]]},R2],Circle[{a1[[4]],a3[[4]]},R2],Circle[{a1[[5]],a3[[5]]},R2],Circle[{a1[[6]],a3[[6]]},R2],Circle[{a1[[7]],a3[[7]]},R4],Circle[{a1[[8]],a3[[8]]},R4],Disk[{b1[[1]],b3[[1]]},0.08], Disk[{b1[[2]],b3[[2]]},0.08], Disk[{b1[[3]],b3[[3]]},0.08], Disk[{b1[[4]],b3[[4]]},0.08],Disk[{b1[[5]],b3[[5]]},0.08],Disk[{b1[[6]],b3[[6]]},0.08]} ] /. t[[2]]

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