(* [NaCl6/6]infG Kimball, Ne/Ar-centers,  17.2.2012 *)
Clear[k1,k2,k3,k4,k5,sig1,sig2,sig3,sig4,sig5,c,z,R1,R2,R3,R4,R5,S2,S3,S4];
c = {k1 -> 0.99000000, k2 -> 2.65000000, k3 -> 0.99000000, k4 -> 1.22000000, k5 -> 3.42000000, 
sig1 -> 0.3000000, sig2 -> 0.3000000, sig3 -> 0.3000000, sig4 -> 0.3000000, sig5 -> 0.3000000};
z=11.0000000000; z1=17.0000000000; sq23=1.5874010519682;
ad = Sqrt[3./8.]; M = 1.74756000000000000; (* Madelung constant for rocksalt structure *)

(* He and Neon shell of Na+ *)
T = 2.25*k1/R1^2.+9.0*k2/R2^2. /. c;
Vee = 3.0*sig1/R1+12.0*sig2/R2+16.0/(R1+R2)+24.0*ad/(R1+R2) /. c;
Vne = -3.0*z/R1-8.0*z/(R1+R2);
S2 = R2*sq23;

(* He and Neon shell of Cl- *)
T = T + 2.25*k3/R3^2.+9.0*k4/R4^2. /. c;
Vee = Vee + 3.0*sig3/R3+12.0*sig4/R4+16.0/(R3+R4)+24.0*ad/(R3+R4) /. c;
Vne = Vne -3.0*z1/R3-8.0*z1/(R3+R4);
S3 = R4*sq23;

(* Argon shell, full for Cl- *)
T = T + 9.0*k5/R5^2. /. c;
Vee = Vee+12.0*sig5/R5+80.0/(S3+R5)+24.0*ad/(S3+R5) /. c;
Vne = Vne-8.35*z1/(S3+R5);
S4 = R5*sq23;

(* Na+Cl- ionpair; these terms sum up to -1/d *)
d=S2+S4 /. c;
(* Vnn=z*z1/d;
Veeg=(z-1.0)*(z1+1.0)/d;
Vneg=-(z*(z1+1.0)+(z-1.0)*z1)/d;
Vee=Vee+Veeg;
Vne=Vne+Vneg; *)
func = T + Vee + Vne - M/d;

t = FindMinimum[func, {R1,0.136}, {R2,1.131},{R3,0.082228}, {R4,0.33557}, {R5,2.2276785},
{Method -> Automatic}, {MaxIterations -> 500}, {WorkingPrecision -> 14}]

N[Vne /. c /. t[[2]],10]
N[Vee /. c /. t[[2]],10]
N[Vnn /. c /. t[[2]],10]
N[-(Vee + Vne - M/d)/T /. c /. t[[2]],10]
N[d*0.529177 /. c /. t[[2]],10]
N[S2*0.529177 /. c /. t[[2]],10]
N[S4*0.529177 /. c /. t[[2]],10]