(* 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[      7409,        227]
NotebookOptionsPosition[      6535,        196]
NotebookOutlinePosition[      7047,        215]
CellTagsIndexPosition[      7004,        212]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
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.2650207; sig2=0.30451569; sig4=sig2; nc=5; pi=0.0;
k1=1.04605075; k2=1.20889711; k4=k2; bohr=0.529177; rad=57.29578;

(* C He shell, i=j each *)
Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
vne=-3.0*z/R1;

(* 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};        
xw={d1,-d4,-d4,d4,d4};
yw={0,d5,-d5,0,0};
zw={0,0,0,d5,-d5};
vne=vne-4*(3-(p/R2)^2)/R2;
vnn = 0.0;
(*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]]]

(*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]]]

(*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]]]]

Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.26},{R2,1.24},{p,0.54},{Method -> Automatic}, \
{MaxIterations -> 500}]
(* 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.5784886449106407`*^9, 3.5784887143151627`*^9}}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "40.56089226557408`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.2653799183950385`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.2343053903057881`"}], ",", 
     RowBox[{"p", "\[Rule]", "0.5334356394795194`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398212075`*^9}],

Cell[BoxData[
 RowBox[{"-", "121.15886714387797`"}]], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}],

Cell[BoxData["26.425381110909388`"], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}],

Cell[BoxData["13.611701083963231`"], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}],

Cell[BoxData["1.999999989698028`"], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}],

Cell[BoxData["1.075880843995231`"], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}],

Cell[BoxData["109.4712215648118`"], "Output",
 CellChangeTimes->{{3.578488427415059*^9, 3.578488454746307*^9}, 
   3.5784887398368073`*^9}]
}, Open  ]],

Cell[CellGroupData[{

Cell["\<\
plot1=Graphics[{Circle[{xw[[1]],yw[[1]]},R1],      \
Circle[{xw[[2]],yw[[2]]},R2],Circle[{xw[[3]],yw[[3]]},R2],Circle[{xw[[4]],yw[[\
4]]},R2],Circle[{xw[[5]],yw[[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}]

plot2=Graphics[{Circle[{xw[[1]],zw[[1]]},R1],      \
Circle[{xw[[2]],zw[[2]]},R2],Circle[{xw[[3]],zw[[3]]},R2],Circle[{xw[[4]],zw[[\
4]]},R2],Circle[{xw[[5]],zw[[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.578488444762289*^9, 3.578488449660698*^9}}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0., 0}, 0.2653799183950385], 
   CircleBox[{-0.8658437166781492, 1.2244879270217663`}, 1.2343053903057881`],
    CircleBox[{-0.8658437166781492, -1.2244879270217663`}, 
    1.2343053903057881`], 
   CircleBox[{0.8658437166781492, 0}, 1.2343053903057881`], 
   CircleBox[{0.8658437166781492, 0}, 1.2343053903057881`], 
   DiskBox[{0., 0}, 0.08], 
   DiskBox[{-1.17382292672699, 1.660036302801789}, 0.08], 
   DiskBox[{-1.17382292672699, -1.660036302801789}, 0.08], 
   DiskBox[{1.17382292672699, 0}, 0.08], DiskBox[{1.17382292672699, 0}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{{3.5784884274618587`*^9, 3.5784884548087068`*^9}, 
   3.578488739883608*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0., 0}, 0.2653799183950385], 
   CircleBox[{-0.8658437166781492, 0}, 1.2343053903057881`], 
   CircleBox[{-0.8658437166781492, 0}, 1.2343053903057881`], 
   CircleBox[{0.8658437166781492, 1.2244879270217663`}, 1.2343053903057881`], 
   CircleBox[{0.8658437166781492, -1.2244879270217663`}, 1.2343053903057881`],
    DiskBox[{0., 0}, 0.08], DiskBox[{-1.17382292672699, 0}, 0.08], 
   DiskBox[{-1.17382292672699, 0}, 0.08], 
   DiskBox[{1.17382292672699, 1.660036302801789}, 0.08], 
   DiskBox[{1.17382292672699, -1.660036302801789}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-3, 3}, {-3, 3}}]], "Output",
 CellChangeTimes->{{3.5784884274618587`*^9, 3.5784884548087068`*^9}, 
   3.578488739883608*^9}]
}, Open  ]]
},
WindowToolbars->"EditBar",
WindowSize->{882, 964},
WindowMargins->{{Automatic, 109}, {Automatic, -3}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
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, 1946, 69, 1034, "Input"],
Cell[2528, 93, 417, 11, 31, "Output"],
Cell[2948, 106, 157, 3, 31, "Output"],
Cell[3108, 111, 140, 2, 31, "Output"],
Cell[3251, 115, 140, 2, 31, "Output"],
Cell[3394, 119, 139, 2, 31, "Output"],
Cell[3536, 123, 139, 2, 31, "Output"],
Cell[3678, 127, 139, 2, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[3854, 134, 1058, 22, 337, "Input"],
Cell[4915, 158, 803, 17, 379, "Output"],
Cell[5721, 177, 798, 16, 411, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)