(* [KCl6/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.99, k2 -> 1.2, k3 -> 3.1, k4 -> 0.99, k5 -> 1.218, k6 -> 3.26, k6->3.5, sig1 -> 0.3, sig2 -> 0.3, sig3 -> 0.3, sig4 -> 0.3, sig5 -> 0.3, sig6->0.3}; z=19.0; z1=17.0;
ad = Sqrt[3./8.]; M = 1.7476; (* Madelung constant for rocksalt structure *)

(* He and Neon shell of K+ *)
T = 2.25*k1/R1^2+9.0*k2/R2^2 /. c;
Vee = 3.0*sig1/R1+12.0*sig2/R2+16/(R1+R2)+24*ad/(R1+R2) /. c;
Vne = -3.0*z/R1-8.0*z/(R1+R2);
S2 = R2*4^(1/3);

(* Argon shell, for K+ *)
T = T + 9.0*k3/R3^2 /. c;
Vee = Vee+12.0*sig3/R3+80.0/(S2+R3)+24.0*ad/(S2+R3) /. c;
Vne = Vne-8.35*z/(S2+R3);
S3 = R3*4^(1/3);

(* He and Neon shell of Cl- *)
T = T + 2.25*k4/R4^2+9.0*k5/R5^2 /. c;
Vee = Vee + 3.0*sig4/R4+12.0*sig5/R5+16/(R4+R5)+24*ad/(R4+R5) /. c;
Vne = Vne -3.0*z1/R4-8.0*z1/(R4+R5);
S4 = R5*4^(1/3);

(* Argon shell, full for Cl- *)
T = T + 9.0*k6/R6^2 /. c;
Vee = Vee+12.0*sig6/R6+80.0/(S4+R6)+24.0*ad/(S4+R6) /. c;
Vne = Vne-8.35*z1/(S4+R6);
S5 = R6*4^(1/3);

(* K+Cl- ionpair; these terms sum up to -1/d *)
d=S3+S5 /. c;
(* Vnn=z*z1/d;
Veeg=(z-1)*(z1+1)/d;
Vneg=-(z*(z1+1)+(z-1)*z1)/d;
Vee:=Vee+Veeg;
Vne:=Vne+Vneg;*)
func = T + Vee + Vne - M/d;

(* the start values have been determined by several attempts
going back and forth between FindMinimum and its results. Criterium
used is the virial ratio, now 1.999998824. Can be improved by
one or two cycles more. perhaps, there is an automatism for this. *)

t = FindMinimum[func, {R1,0.07291497005}, {R2,0.2862678358},{R3,1.584372841},{R4,0.08221147752}, {R5,0.3349643624}, {R6,2.159407892}, {Method -> Automatic},
{MaxIterations -> 500}]

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