(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[      4304,        141]
NotebookOptionsPosition[      3480,        113]
NotebookOutlinePosition[      4178,        136]
CellTagsIndexPosition[      4135,        133]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["\<\
(* [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]\
\>", "Input",
 FontSize->16],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "995.4501913154753`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.07291461118526976`"}], ",", 
     RowBox[{"R2", "\[Rule]", "0.2862675276818523`"}], ",", 
     RowBox[{"R3", "\[Rule]", "1.5843786484186706`"}], ",", 
     RowBox[{"R4", "\[Rule]", "0.0822115068150589`"}], ",", 
     RowBox[{"R5", "\[Rule]", "0.33496350885566767`"}], ",", 
     RowBox[{"R6", "\[Rule]", "2.158472108171909`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{3.5731429969593143`*^9},
 FontSize->16],

Cell[BoxData["1.9999999999662796`"], "Output",
 CellChangeTimes->{3.5731429969593143`*^9},
 FontSize->16],

Cell[BoxData["3.144054994534134`"], "Output",
 CellChangeTimes->{3.5731429969593143`*^9},
 FontSize->16],

Cell[BoxData["1.330903615118112`"], "Output",
 CellChangeTimes->{3.5731429969749146`*^9},
 FontSize->16],

Cell[BoxData["1.8131513794160223`"], "Output",
 CellChangeTimes->{3.5731429969749146`*^9},
 FontSize->16]
}, Open  ]]
},
PrintingStyleEnvironment->"Printout",
WindowSize->{1558, 940},
WindowMargins->{{100, Automatic}, {-7, Automatic}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
PrintingOptions->{"PrintCellBrackets"->False,
"PrintMultipleHorizontalPages"->False,
"PrintRegistrationMarks"->False,
"PrintingMargins"->{{34, 14.125}, {56.6875, 56.6875}}},
PrivateNotebookOptions->{"VersionedStylesheet"->{"Default.nb"[8.] -> False}},
FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (January 25, 2013)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[579, 22, 1885, 56, 964, "Input"],
Cell[2467, 80, 567, 14, 33, "Output"],
Cell[3037, 96, 105, 2, 33, "Output"],
Cell[3145, 100, 104, 2, 33, "Output"],
Cell[3252, 104, 104, 2, 33, "Output"],
Cell[3359, 108, 105, 2, 65, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)