(* Al atom Kimball, Ne-tbp core 27.12.2012 *)

Clear[k1,k2,k3,k4, sig1,sig2,c,z,R1,R2,R3,R4,a,b,f,g,d];
c = {k1 -> 1.0, k2 -> 1.0, k3 -> 1.0, k4 -> 1.0, sig1 -> 0.28, sig2 -> 0.3, sig3 -> 0.3};
z=13.0; 
ad = Sqrt[3]/2;
T = 2.25*k1/R1^2+6.75*k2/R2^2+3.375*k3/R3^2+2.25*k4/R4^2 /. c;  
b = (R1+R2)/2+Sqrt[(R2+R3)^2-3*(R1+R2)^2/4];
a = (R1+R2)*Sqrt[3];d = b*Sqrt[3];
f = Sqrt[(R1+R4)^2+(R1+R2)^2];g=Sqrt[b^2+(R1+R4)^2];Vee=3.0*sig1/R1+3.0*3.0*sig2/R2 /. c;
Vee = Vee + 12/(R1+R2)+6/b+12/a+4/(R1+R4)+12/f+3/d+12/(R2+R3)+6/(b+R1+R2)+6/g+0.5/(R1+R4);
Vne=-3.0*z/R1-6*z/(R1+R2)-2*z/(R1+R4)-3*z/b;

func = T + Vee + Vne /. c;

t = FindMinimum[func, {R1,0.109}, {R2,0.4198}, {R3,1.668},{R4,0.416},{MaxIterations -> 500},{Method-> "Newton"}]

-(Vne+Vee)/T /. t[[2]]

plot1=Graphics3D[{Opacity[0.5],{Lighter[Green,0.2],
Sphere[{-(R1+R2)/2,ad*(R1+R2),0},R2],
Sphere[{-(R1+R2)/2,-ad*(R1+R2),0},R2],
Sphere[{R1+R2,0,0},R2]},
{Darker[Red,1],Sphere[{0,0,0},R1]},
{Lighter[Blue,0.2],
Sphere[{0,0,R1+R4},R4],
Sphere[{0,0,-(R1+R4)},R4]},
Sphere[{-b,0,0},R3],
Sphere[{b/2,ad*b,0},R3], 
Sphere[{b/2,-ad*b,0},R3]}] /. t[[2]];

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