(* MgO: Ne ions nach Kimball, Aussen-sp-Kugeln 28.04.2012 *)
Clear[k1,k2,k3,k4,sig1,sig2,sig3,sig4,c,z,z1];
c = {k1 -> 1.0, k2 -> 1.0, k3 -> 1.1, k4 -> 1.0, sig1 -> 0.3, sig2 -> 0.3, sig3 -> 0.3, sig4 -> 0.3};  z=12.0; z1=8.0; M=1.74756; fmg=2.85; fo=2.85;
ad=Sqrt[3/8]; bohr=0.529177;

(* He and Neon shell of Mg+2 *)
Tmg = 2.25*k1/R1^2+9.0*k2/R2^2 /. c;
Veemg = 3.0*sig1/R1+12.0*sig2/R2+(16+24*ad)/(R1+R2) /. c;
Vnemg = -3.0*z/R1-8.0*z/(R1+R2);
S2 = R2*fmg;
Emg = Tmg+Veemg+Vnemg;

(* He and Neon shell of O-2 *)
To2 = 2.25*k3/R3^2+9.0*k4/R4^2 /. c;
Veeo = 3.0*sig3/R3+12.0*sig4/R4+(16+24*ad)/(R3+R4) /. c;
Vneo = -3.0*z1/R3-8.0*z1/(R3+R4);
S3 = R4*fo;
Eom = To2+Veeo+Vneo;

(* Mg+2 O-2 ionpair *)
d=S2+S3;
T=Tmg+To2;
Vee=Veemg+Veeo;
Vne=Vnemg+Vneo;
Epot=Vee+Vne-M*(z-10)*(10-z1)/d+0.372*2.0/S3; (*Hard Sphere model ++*)

func = T + Epot /. c;

t = FindMinimum[func, {R1,0.119}, {R2,0.4705}, {R3,0.2045}, {R4,0.922}, {MaxIterations -> 500}]

Emg /. t[[2]]
Eom /. t[[2]]
(t[[1]] - Emg - Eom /.t[[2]])*627.5095*4.184
bohr*S2 /. c /. t[[2]]
bohr*S3 /. c /. t[[2]]
S2/S3 /. t[[2]]
bohr*(S2+S3) /. c /. t[[2]]
Epot/T /. c /. t[[2]]