(* 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[      3784,        123]
NotebookOptionsPosition[      2968,         95]
NotebookOutlinePosition[      3665,        118]
CellTagsIndexPosition[      3622,        115]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["\<\
(* [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.99, k2 -> 1.98, k3 -> 0.99, k4 -> 1.1, k5 -> 2.83, sig1 -> 0.3, \
sig2 -> 0.3, sig3 -> 0.3, sig4 -> 0.3, sig5 -> 0.3};
(* c = {k1 -> 0.99, k2 -> 2.65, k3 -> 0.99, k4 -> 1.22, k5 -> 3.42, sig1 -> \
0.3, sig2 -> 0.3, sig3 -> 0.3, sig4 -> 0.3, sig5 -> 0.3};*)
z=11.0; z1=17.0;
ad = Sqrt[3./8.]; M = 1.7476; (* 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/(R1+R2)+24*ad/(R1+R2) /. c;
Vne = -3.0*z/R1-8.0*z/(R1+R2);
S2 = R2*4^(1/3);

(* 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/(R3+R4)+24*ad/(R3+R4) /. c;
Vne = Vne -3.0*z1/R3-8.0*z1/(R3+R4);
S3 = R4*4^(1/3);

(* 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*4^(1/3);

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

func = T + Vee + Vne - M/d;

t = FindMinimum[func, {R1,0.134369285}, {R2,0.8846053088},{R3,0.08158192127}, \
{R4,0.3116682854}, {R5,2.149202772}, {Method -> Automatic},
{MaxIterations -> 500}]

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]\
\>", "Input"],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "589.8332663222693`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.1346744461330032`"}], ",", 
     RowBox[{"R2", "\[Rule]", "0.8952646740677106`"}], ",", 
     RowBox[{"R3", "\[Rule]", "0.08151922191546108`"}], ",", 
     RowBox[{"R4", "\[Rule]", "0.3106646614080302`"}], ",", 
     RowBox[{"R5", "\[Rule]", "1.9056499555167568`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{3.573143175548428*^9}],

Cell[BoxData["2.0000000002175264`"], "Output",
 CellChangeTimes->{3.573143175564028*^9}],

Cell[BoxData["2.352813457737358`"], "Output",
 CellChangeTimes->{3.573143175564028*^9}],

Cell[BoxData["0.7520367636823887`"], "Output",
 CellChangeTimes->{3.573143175564028*^9}],

Cell[BoxData["1.6007766940549697`"], "Output",
 CellChangeTimes->{3.573143175564028*^9}]
}, Open  ]]
},
PrintingStyleEnvironment->"Printout",
WindowSize->{802, 940},
WindowMargins->{{Automatic, 437}, {Automatic, -9}},
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, 1518, 44, 694, "Input"],
Cell[2100, 68, 489, 12, 52, "Output"],
Cell[2592, 82, 88, 1, 31, "Output"],
Cell[2683, 85, 87, 1, 31, "Output"],
Cell[2773, 88, 88, 1, 31, "Output"],
Cell[2864, 91, 88, 1, 63, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)