(* 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[     12976,        419]
NotebookOptionsPosition[     11125,        357]
NotebookOutlinePosition[     11637,        376]
CellTagsIndexPosition[     11594,        373]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* CH4 methane.  with k,sig of CH41401.pas and minim. ident. num. results \
17.06.2012 *)
Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,
a1,a2,a3,b1,b2,b3,c1,c2,c3,cs,ss,d1,d2,d3,d4,d5,pi,i,j,d,t];
z=6.0; sig1=0.3; sig2=0.3; sig4=sig2; nc=5; pi=0.0;
k1=1.027; k2=1.23; k4=k2; bohr=0.529177; rad=57.29578;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{3.5668390630515456`*^9}],

Cell["\<\
(* C He shell, i=j each *)
Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
vne=-3.0*z/R1;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390682151546`*^9}}],

Cell["\<\
(* R1=0.262361; R2=1.246136; p=0.539862264; *)
R4=R2;
(* bonding pairs *)
Ekin = Ekin + 2.25*(2*k2/R2^2+2*k4/R4^2);
vee=vee+3.0*(2*sig2/R2+2*sig4/R4);
oc={-2,-2,-2,-2,-2};
ch={z,1,1,1,1};
rr={R1,R2,R2,R2,R2};
w=ArcCos[-1/3]/2; a=4*(R1+R2)/Sqrt[6];
cs=Cos[w];  ss=Sin[w];  (* w is half tetraeder angle *)
(* d1=R1+R3;  d1=(R1+R2)*cs; *)  (* R1+R2+p is radius of outer sphere *) \
d1=0.0;
d2=d1+(R1+R2+p)*cs; d3=(R1+R2+p)*ss; d4=d1+(R1+R2)*cs;  d5=(R1+R2)*ss;
xn={d1,-d2,-d2,d2,d2};
yn={0,d3,-d3,0,0};
zn={0,0,0,d3,-d3};        
xc={d1,-d4,-d4,d4,d4};
yc={0,d5,-d5,0,0};
zc={0,0,0,d5,-d5};
vne=vne-4*(3-(p/R2)^2)/R2;
vnn = 0.0;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839078760773*^9}, {
  3.566839318626795*^9, 3.5668393238528037`*^9}}],

Cell["\<\
(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + \
oc[[i]]*oc[[j]]/Sqrt[(xw[[i]]-xw[[j]])^2+(yw[[i]]-yw[[j]])^2+(zw[[i]]-zw[[j]])\
^2]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839083924382*^9}}],

Cell["\<\
(*nn*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vnn = vnn + \
ch[[i]]*ch[[j]]/Sqrt[(xn[[i]]-xn[[j]])^2+(yn[[i]]-yn[[j]])^2+(zn[[i]]-zn[[j]])\
^2]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390890567913`*^9}}],

Cell["\<\
(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nc+1, j++,
    If[i != j,
      vne = vne + \
oc[[i]]*ch[[j]]/Sqrt[(xw[[i]]-xn[[j]])^2+(yw[[i]]-yn[[j]])^2+(zw[[i]]-zn[[j]])\
^2]]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668390979332066`*^9}}],

Cell[CellGroupData[{

Cell["\<\
Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.26},{R2,1.24},{p,0.54},{Method -> Automatic}, \
{MaxIterations -> 500}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839104968819*^9}}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "40.69904362506576`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.2623609956342706`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.246135949525892`"}], ",", 
     RowBox[{"p", "\[Rule]", "0.5398621903909542`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{3.5668391985689836`*^9, 3.5668393324172187`*^9, 
  3.566840202945548*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* N[func, 14] *)
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
(R1+R2+p)*bohr /. t[[2]]
2*w*rad /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839112846833*^9}}],

Cell[BoxData[
 RowBox[{"-", "121.57720036929429`"}]], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202945548*^9}],

Cell[BoxData["26.668672410019127`"], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202961148*^9}],

Cell[BoxData["13.510440690728258`"], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202961148*^9}],

Cell[BoxData["1.9999999995475228`"], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202961148*^9}],

Cell[BoxData["1.0839445422735332`"], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202961148*^9}],

Cell[BoxData["109.4712215648118`"], "Output",
 CellChangeTimes->{3.5668391985845833`*^9, 3.5668393324172187`*^9, 
  3.566840202961148*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Hvnfx: sum of H2-nuclei forces *)
Hvnfx = 0.0;
  For[j = 1, j < 6, j++,
    If[j != 2,
    Hvnfx = Hvnfx - \
(xn[[j]]-xn[[2]])*ch[[3]]*ch[[j]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(\
zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefx: sum of H2 - clouds forces *)
Hvefx = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Hvefx = Hvefx- \
(xc[[j]]-xn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(\
zc[[j]]-zn[[2]])^2)^(3/2)]]
Hforcetotx = Hvnfx+Hvefx/. t[[2]]
Hforcennx=Hvnfx /. t[[2]]
HForcenex=Hvefx /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839137994077*^9}}],

Cell[BoxData[
 RowBox[{"-", "0.3221484194569486`"}]], "Output",
 CellChangeTimes->{3.566839198787384*^9, 3.5668393324484186`*^9, 
  3.5668402029767475`*^9}],

Cell[BoxData[
 RowBox[{"-", "0.9520128128756769`"}]], "Output",
 CellChangeTimes->{3.566839198787384*^9, 3.5668393324484186`*^9, 
  3.5668402029767475`*^9}],

Cell[BoxData["0.6298643934187282`"], "Output",
 CellChangeTimes->{3.566839198787384*^9, 3.5668393324484186`*^9, 
  3.5668402029767475`*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Hvnfy: sum of H2-nuclei forces *)
Hvnfy = 0.0;
  For[j = 1, j < 6, j++,
    If[j != 2,
    Hvnfy = Hvnfy - \
(yn[[j]]-yn[[2]])*ch[[2]]*ch[[j]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(\
zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefy: sum of H2 - clouds forces *)
Hvefy = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Hvefy = Hvefy - \
(yc[[j]]-yn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(\
zc[[j]]-zn[[2]])^2)^(3/2)]]
Hforcetoty = Hvnfy+Hvefy /. t[[2]]
Hforcenny=Hvnfy /. t[[2]]
HForceney=Hvefy /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.566839148118495*^9}}],

Cell[BoxData["0.45558666389307323`"], "Output",
 CellChangeTimes->{3.5668391989589844`*^9, 3.5668393324484186`*^9, 
  3.5668402029923477`*^9}],

Cell[BoxData["1.3463494315217415`"], "Output",
 CellChangeTimes->{3.5668391989589844`*^9, 3.5668393324484186`*^9, 
  3.5668402029923477`*^9}],

Cell[BoxData[
 RowBox[{"-", "0.8907627676286681`"}]], "Output",
 CellChangeTimes->{3.5668391989589844`*^9, 3.5668393324484186`*^9, 
  3.566840203007948*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
(* Wvnfy: sum of w2-nuclei forces *)
Wvnfy = 0.0;
  For[j = 1, j < 6, j++,
    Wvnfy = Wvnfy - \
(yn[[j]]-zc[[2]])*oc[[2]]*ch[[j]]/((xc[[2]]-xn[[j]])^2+(yc[[2]]-yn[[j]])^2+(\
zc[[2]]-zn[[j]])^2)^(3/2)]

(* Wvefy: sum of W2 - clouds forces *)
Wvefy = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j != 2,
    Wvefy = Wvefy - \
(yc[[j]]-zc[[2]])*oc[[j]]*oc[[2]]/((xc[[j]]-xc[[2]])^2+(yc[[j]]-yc[[2]])^2+(\
zc[[j]]-zc[[2]])^2)^(3/2)]]
Wforcetoty = Wvnfy+Wvefy /. t[[2]]
Wforcenny=Wvnfy /. t[[2]]
WForceney=Wvefy /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668390630515456`*^9, 3.5668391587889137`*^9}}],

Cell[BoxData["21.454299348040824`"], "Output",
 CellChangeTimes->{3.566839199286585*^9, 3.566839332464019*^9, 
  3.566840203007948*^9}],

Cell[BoxData["21.12471058739574`"], "Output",
 CellChangeTimes->{3.566839199286585*^9, 3.566839332464019*^9, 
  3.566840203023548*^9}],

Cell[BoxData["0.32958876064508486`"], "Output",
 CellChangeTimes->{3.566839199286585*^9, 3.566839332464019*^9, 
  3.566840203023548*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},R1],      \
Circle[{xc[[2]],yc[[2]]},R2],Circle[{xc[[3]],yc[[3]]},R2],Circle[{xc[[4]],yc[[\
4]]},R2],Circle[{xc[[5]],yc[[5]]},R2],Disk[{xn[[1]],yn[[1]]},0.08], \
Disk[{xn[[2]],yn[[2]]},0.08], Disk[{xn[[3]],yn[[3]]},0.08], \
Disk[{xn[[4]],yn[[4]]},0.08],Disk[{xn[[5]],yn[[5]]},0.08]} ] /. t[[2]];

Show[plot1,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-3,3},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.5668391835773573`*^9}, {
  3.566840085851742*^9, 3.566840129656619*^9}}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0., 0}, 0.2623609956342706], 
   CircleBox[{-0.870931117359948, 1.2316825980631922`}, 1.246135949525892], 
   CircleBox[{-0.870931117359948, -1.2316825980631922`}, 1.246135949525892], 
   CircleBox[{0.870931117359948, 0}, 1.246135949525892], 
   CircleBox[{0.870931117359948, 0}, 1.246135949525892], 
   DiskBox[{0., 0}, 0.08], 
   DiskBox[{-1.1826206983074663`, 1.6724782306895591`}, 0.08], 
   DiskBox[{-1.1826206983074663`, -1.6724782306895591`}, 0.08], 
   DiskBox[{1.1826206983074663`, 0}, 0.08], 
   DiskBox[{1.1826206983074663`, 0}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.5668391993177853`*^9, 3.566839332479619*^9, 
  3.566840203023548*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},R1],      \
Circle[{xc[[2]],zc[[2]]},R2],Circle[{xc[[3]],zc[[3]]},R2],Circle[{xc[[4]],zc[[\
4]]},R2],Circle[{xc[[5]],zc[[5]]},R2],Disk[{xn[[1]],zn[[1]]},0.08], \
Disk[{xn[[2]],zn[[2]]},0.08], Disk[{xn[[3]],zn[[3]]},0.08], \
Disk[{xn[[4]],zn[[4]]},0.08],Disk[{xn[[5]],zn[[5]]},0.08]} ] /. t[[2]];

Show[plot2,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-3,3},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.5668391881949654`*^9}, {
  3.566840146036648*^9, 3.566840172307094*^9}}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0., 0}, 0.2623609956342706], 
   CircleBox[{-0.870931117359948, 0}, 1.246135949525892], 
   CircleBox[{-0.870931117359948, 0}, 1.246135949525892], 
   CircleBox[{0.870931117359948, 1.2316825980631922`}, 1.246135949525892], 
   CircleBox[{0.870931117359948, -1.2316825980631922`}, 1.246135949525892], 
   DiskBox[{0., 0}, 0.08], DiskBox[{-1.1826206983074663`, 0}, 0.08], 
   DiskBox[{-1.1826206983074663`, 0}, 0.08], 
   DiskBox[{1.1826206983074663`, 1.6724782306895591`}, 0.08], 
   DiskBox[{1.1826206983074663`, -1.6724782306895591`}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.566839199364585*^9, 3.5668393324952188`*^9, 
  3.566840203054748*^9}]
}, Open  ]],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668391703953342`*^9, 3.566839177774147*^9}}]
},
WindowToolbars->"EditBar",
WindowSize->{993, 964},
WindowMargins->{{Automatic, 392}, {Automatic, 0}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
PrivateNotebookOptions->{"VersionedStylesheet"->{"Default.nb"[8.] -> False}},
FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (November 20, 2012)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 415, 9, 99, "Input"],
Cell[975, 31, 199, 7, 82, "Input"],
Cell[1177, 40, 794, 25, 354, "Input"],
Cell[1974, 67, 282, 9, 82, "Input"],
Cell[2259, 78, 284, 9, 82, "Input"],
Cell[2546, 89, 302, 10, 99, "Input"],
Cell[CellGroupData[{
Cell[2873, 103, 250, 7, 65, "Input"],
Cell[3126, 112, 415, 11, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[3578, 128, 241, 10, 133, "Input"],
Cell[3822, 140, 156, 3, 31, "Output"],
Cell[3981, 145, 139, 2, 31, "Output"],
Cell[4123, 149, 139, 2, 31, "Output"],
Cell[4265, 153, 139, 2, 31, "Output"],
Cell[4407, 157, 139, 2, 31, "Output"],
Cell[4549, 161, 138, 2, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[4724, 168, 644, 21, 286, "Input"],
Cell[5371, 191, 156, 3, 31, "Output"],
Cell[5530, 196, 156, 3, 31, "Output"],
Cell[5689, 201, 139, 2, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[5865, 208, 646, 21, 286, "Input"],
Cell[6514, 231, 142, 2, 31, "Output"],
Cell[6659, 235, 141, 2, 31, "Output"],
Cell[6803, 239, 156, 3, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[6996, 247, 632, 20, 269, "Input"],
Cell[7631, 269, 135, 2, 31, "Output"],
Cell[7769, 273, 134, 2, 31, "Output"],
Cell[7906, 277, 136, 2, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[8079, 284, 638, 12, 133, "Input"],
Cell[8720, 298, 801, 17, 379, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[9558, 320, 638, 12, 133, "Input"],
Cell[10199, 334, 797, 16, 379, "Output"]
}, Open  ]],
Cell[11011, 353, 110, 2, 31, "Input"]
}
]
*)

(* End of internal cache information *)