(* SiO2: Ne ions after Kimball,  7.05.2012 *)
Clear[k1,k2,k3,k4,sig1,sig2,sig3,sig4,R1,R2,R3,R4,c,z,z1];
c = {k1 -> 0.902, k2 -> 1.45, k3 -> 0.98, k4 -> 2.46, 
sig1 -> 0.3, sig2 -> 0.3, sig3 -> 0.3, sig4 -> 0.3};  
z=14.0; z1=8.0; M=4.4394; (* Madelung const. of beta-Cristobalite *)
ad=Sqrt[3/8]; bohr=0.529177; rf = 4^(1/3);

(* He and Neon shell of Si+4 *)
Tsi = 2.25*k1/R1^2+9.0*k2/R2^2 /. c;
Veesi = 3.0*sig1/R1+12.0*sig2/R2+(16+24*ad)/(R1+R2) /. c;
Vnesi = -3.0*z/R1-8.0*z/(R1+R2) /. c;
S2 = R2*rf;
Esi=Tsi+Veesi+Vnesi;

(* 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) /. c;
S3 = R4*rf;
Eom=To2+Veeo+Vneo;

Vee=Veesi+2*Veeo;
Vne=Vnesi+2*Vneo;
Ekin=Tsi+2*To2;

(* Si+4 O-2 ionpair and crystal electrostatic interaction *)
d=S2+S3;
Epot=Vee+Vne-M*(z-10)/d;
Etot=Ekin+Epot  /. c;

t = FindMinimum[Etot, {R1,0.09385}, {R2,0.48347}, {R3,0.188}, {R4,1.6625}, 
    {MaxIterations -> 500}]

Esi /. t[[2]]
Eom /. t[[2]]
(t[[1]] - Esi - 2*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/Ekin /. c /. t[[2]]