(* CH4 methane.  with k,sig of CH41401.pas and minim. ident. num. results 17.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=sig2; nc=5; pi=0.0;
k1=1.027; k2=1.23; k4=k2; bohr=0.529177; rad=57.29578;

(* C He shell, i=j each *)
Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
vne=-3.0*z/R1;

(* R1=0.262361; R2=1.246136; p=0.539862264; *)
R4=R2;
(* bonding pairs *)
Ekin = Ekin + 2.25*(2*k2/R2^2+2*k4/R4^2);
vee=vee+3.0*(2*sig2/R2+2*sig4/R4);
oc={-2,-2,-2,-2,-2};
ch={z,1,1,1,1};
rr={R1,R2,R2,R2,R2};
w=ArcCos[-1/3]/2; a=4*(R1+R2)/Sqrt[6];
cs=Cos[w];  ss=Sin[w];  (* w is half tetraeder angle *)
(* d1=R1+R3;  d1=(R1+R2)*cs; *)  (* R1+R2+p is radius of outer sphere *) d1=0.0;
d2=d1+(R1+R2+p)*cs; d3=(R1+R2+p)*ss; d4=d1+(R1+R2)*cs;  d5=(R1+R2)*ss;
xn={d1,-d2,-d2,d2,d2};
yn={0,d3,-d3,0,0};
zn={0,0,0,d3,-d3};        
xw={d1,-d4,-d4,d4,d4};
yw={0,d5,-d5,0,0};
zw={0,0,0,d5,-d5};
vne=vne-4*(3-(p/R2)^2)/R2;
vnn = 0.0;
(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + oc[[i]]*oc[[j]]/Sqrt[(xw[[i]]-xw[[j]])^2+(yw[[i]]-yw[[j]])^2+(zw[[i]]-zw[[j]])^2]]]

(*nn*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vnn = vnn + ch[[i]]*ch[[j]]/Sqrt[(xn[[i]]-xn[[j]])^2+(yn[[i]]-yn[[j]])^2+(zn[[i]]-zn[[j]])^2]]]

(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nc+1, j++,
    If[i != j,
      vne = vne + oc[[i]]*ch[[j]]/Sqrt[(xw[[i]]-xn[[j]])^2+(yw[[i]]-yn[[j]])^2+(zw[[i]]-zn[[j]])^2]]]]

Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.26},{R2,1.24},{p,0.54},{Method -> Automatic}, {MaxIterations -> 500}]
(* N[func, 14] *)
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
(R1+R2+p)*bohr /. t[[2]]
2*w*rad /. t[[2]]

plot1=Graphics[{Circle[{xw[[1]],yw[[1]]},R1],      Circle[{xw[[2]],yw[[2]]},R2],Circle[{xw[[3]],yw[[3]]},R2],Circle[{xw[[4]],yw[[4]]},R2],Circle[{xw[[5]],yw[[5]]},R2],Disk[{xn[[1]],yn[[1]]},0.08], Disk[{xn[[2]],yn[[2]]},0.08], Disk[{xn[[3]],yn[[3]]},0.08], Disk[{xn[[4]],yn[[4]]},0.08],Disk[{xn[[5]],yn[[5]]},0.08]} ] /. t[[2]]

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

plot2=Graphics[{Circle[{xw[[1]],zw[[1]]},R1],      Circle[{xw[[2]],zw[[2]]},R2],Circle[{xw[[3]],zw[[3]]},R2],Circle[{xw[[4]],zw[[4]]},R2],Circle[{xw[[5]],zw[[5]]},R2],Disk[{xn[[1]],zn[[1]]},0.08], Disk[{xn[[2]],zn[[2]]},0.08], Disk[{xn[[3]],zn[[3]]},0.08], Disk[{xn[[4]],zn[[4]]},0.08],Disk[{xn[[5]],zn[[5]]},0.08]} ] /. t[[2]]

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