(* 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[     16558,        460]
NotebookOptionsPosition[     15029,        411]
NotebookOutlinePosition[     15569,        431]
CellTagsIndexPosition[     15526,        428]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* H3C-CH3 with  two fused methanes at vertex with CH4 data CH4Mt(min).nb  \
01.07.2012 *)
Clear[z,sig1,sig2,sig4,k1,k2,k4,nc,R1,R2,R3,R4,w,p,vee,vne,vnn,
xc,yc,zc,xn,yn,zn,oc,ch,rr,cs,ss,d1,d2,d3,d4,d5,pi,i,j,d,t];
z=6.0;
nc=9;   (* number of clouds *)
sig1=0.3; sig2=0.3; sig4=sig2; (* screening const. from e-e interaction in \
doubly occ. clouds *)
k1=1.027; k2=1.23; k4=k2; (* parameters for kinetic energy of clouds; k=1.0 \
Kimball's lowest value *)
bohr=0.529177; rad=57.29578;
\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, 
   3.5657594185941954`*^9},
 FontSize->16],

Cell["\<\
(* C He-shells *)
Ekin = 2*(2.25*k1/R1^2);
vee=2*(3.0*sig1/R1);
vne=-2*(3.0*z/R1);\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, 
   3.5657594185941954`*^9},
 FontSize->16],

Cell["\<\
(* this is the common vertex assumption *)
(* R4=R2; *)
\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594297326155`*^9}, {3.56575972282593*^9, 
  3.5657597299551425`*^9}},
 FontSize->16],

Cell["\<\
(* bonding pairs *)
Ekin = Ekin + 2.25*(6*k2/R2^2+k4/R4^2);
vee=vee+3.0*(6*sig2/R2+sig4/R4);

(* cloud occupation *)
oc={-2,-2,-2,-2,-2,-2,-2,-2,-2};

(* nuclear charges for C1,C2,H3,H4,H5,H6,H7,H8,vertex9 *)
ch={6,6,1,1,1,1,1,1,0};

(* cloud radii in the same order *)
rr={R1,R1,R2,R2,R2,R2,R2,R2,R4};

(* w is half angle between two C-H of CH4, i.e. 109.47\[Degree]/2 *)
w=ArcCos[-1/3]/2;
cs=Cos[w];  
ss=Sin[w];

(* edge length of tetrahedron of 4 equal clouds *)
a=4*(R1+R2)/Sqrt[6]; (* 4/Sqrt[6] is also Sqrt[8/3] *)

(* x is C-C bond axis, xz one mirrorplane of D3d molecule *)
(* R1+R2 is radius of outer sphere for each tetrahedron of equal clouds *)
d1=R1+R4; 
d2=R1+R2;
  
(* cloud coordinates in terms of radii *)        
d3=d1+d2/3; 
d4=d2*Sqrt[2/3];  
d5=d2*Sqrt[8]/3;
xc={-d1,d1,-d3,-d3,-d3,d3,d3,d3,0};
yc={0,0,0,-d4,d4,0,d4,-d4,0};
zc={0,0,d5,-d5/2,-d5/2,-d5,d5/2,d5/2,0}; (* staggered conformation *)

(* nuclear coordinates in terms of radii; C nuclei assumed in center of C(1s) \
cloud *)
(* CH3 units span tetrahedron of circumsphere radius R1+R2+p, if tetrahedral! \
*)
d6=R1+R2+p; 
d7=d1+d6/3; 
d8=d6*Sqrt[2/3]; 
d9=d6*Sqrt[8]/3;
xn={-d1,d1,-d7,-d7,-d7,d7,d7,d7,0};
yn={0,0,0,d8,-d8,0,d8,-d8,0};
zn={0,0,d9,-d9/2,-d9/2,-d9,d9/2,d9/2,0};

(* potential energy of protons in CH-clouds with eccentricity p *)
vne=vne-6*(3-(p/R2)^2)/R2;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594297326155`*^9}, {3.5657595016331415`*^9, 
  3.565759544658017*^9}, {3.5657596661354303`*^9, 3.565759695307482*^9}},
 FontSize->16],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594509486523`*^9}},
 FontSize->16],

Cell["\<\
(* cc: sum of cloud-cloud potential energies *)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + \
oc[[i]]*oc[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])\
^2]]]

(* nn: sum of nuclei-nuclei potential energies *)
vnn = 0.0;
For[i = 1, i < nc-2, 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]]]

(* cn: sum of cloud-nuclei potential energies *)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < nc+1, j++,
    If[i != j,
      vne = vne + \
oc[[i]]*ch[[j]]/Sqrt[(xc[[i]]-xn[[j]])^2+(yc[[i]]-yn[[j]])^2+(zc[[i]]-zn[[j]])\
^2]]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594509486523`*^9}},
 FontSize->16],

Cell[CellGroupData[{

Cell["\<\
Epot=vne+vee+vnn;
func=Ekin+Epot;
(* results of CH4 computation; if this is not available, decomment the \
minimize function *)
(* R1=0.2623610; R2=1.2461360; p=0.53986226; *)
(* minimization function for R1, R2, p *)
t = FindMinimum[func,{R1,0.2623610},{R2,1.2461360},{R4,1.2461360},{p,0.\
53986226},{Method -> \"Newton\"}, {MaxIterations -> 500}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.565759460339869*^9}, {3.565759795990059*^9, 
  3.5657598023080697`*^9}},
 FontSize->16],

Cell[BoxData[
 RowBox[{
  StyleBox[
   RowBox[{"FindMinimum", "::", "lstol"}], "MessageName"], 
  RowBox[{
  ":", " "}], "\<\"The line search decreased the step size to within the \
tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find \
a sufficient decrease in the function. You may need more than \
\[NoBreak]\\!\\(MachinePrecision\\)\[NoBreak] digits of working precision to \
meet these tolerances. \\!\\(\\*ButtonBox[\\\"\[RightSkeleton]\\\", \
ButtonStyle->\\\"Link\\\", ButtonFrame->None, \
ButtonData:>\\\"paclet:ref/message/FindMinimum/lstol\\\", ButtonNote -> \
\\\"FindMinimum::lstol\\\"]\\)\"\>"}]], "Message", "MSG",
 CellChangeTimes->{3.565759486875516*^9, 3.5657598141484904`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "4.720649941855578`*^15"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.2537825146433071`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.0096318340753785`"}], ",", 
     RowBox[{"R4", "\[Rule]", 
      RowBox[{"-", "0.5583462051844958`"}]}], ",", 
     RowBox[{"p", "\[Rule]", "7.771561172376096`*^-16"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{3.5657594868911157`*^9, 3.5657598141484904`*^9}]
}, Open  ]],

Cell["\<\

(* func *)\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594804639044`*^9}},
 FontSize->16],

Cell[CellGroupData[{

Cell["\<\
vne  /. t[[2]] 
vee  /. t[[2]] 
vnn  /. t[[2]]
-Epot/Ekin  /. t[[2]]
2*d1*bohr  /. t[[2]]
(R1+R2+p)*bohr  /. t[[2]]
(R1+R2)*bohr  /. t[[2]]
2*w*rad  /. t[[2]]
4*(R1*(R1+2*R4))^(3/2)/(R1+R4)^3 /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594804639044`*^9}},
 FontSize->16],

Cell[BoxData[
 RowBox[{"-", "4.720649941855927`*^15"}]], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["139.0331601306798`"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["113.79362455446977`"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["4.870505494199379`*^13"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData[
 RowBox[{"-", "0.3223362001390293`"}]], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["0.6685698148119084`"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["0.6685698148119078`"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData["109.4712215648118`"], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141796904`*^9}],

Cell[BoxData[
 RowBox[{"0.`", "\[VeryThinSpace]", "+", 
  RowBox[{"14.509963863387917`", " ", "\[ImaginaryI]"}]}]], "Output",
 CellChangeTimes->{3.5657594869067154`*^9, 3.5657598141952906`*^9}]
}, Open  ]],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592721567383`*^9, 3.565759315743215*^9}, {
  3.5657594185941954`*^9, 3.5657594804639044`*^9}},
 FontSize->16],

Cell[CellGroupData[{

Cell["\<\
(* projection on xy-plane of molecule *)
plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},R1],      \
Circle[{xc[[2]],yc[[2]]},R1],Circle[{xc[[3]],yc[[3]]},R2],Circle[{xc[[4]],yc[[\
4]]},R2],Circle[{xc[[5]],yc[[5]]},R2],Circle[{xc[[6]],yc[[6]]},R2],Circle[{xc[\
[7]],yc[[7]]},R2],Circle[{xc[[8]],yc[[8]]},R2],Circle[{xc[[9]],yc[[9]]},R4],\
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],Disk[{xn[[6]],yn[[6]\
]},0.08],Disk[{xn[[7]],yn[[7]]},0.08],Disk[{xn[[8]],yn[[8]]},0.08]} ]  /. \
t[[2]]

Show[plot1,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-4,4},{-3,3}}, Frame -> True}]

(* projection on xz-plane, perpendicular to molecular plane *)
plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},R1],      \
Circle[{xc[[2]],zc[[2]]},R1],Circle[{xc[[3]],zc[[3]]},R2],Circle[{xc[[4]],zc[[\
4]]},R2],Circle[{xc[[5]],zc[[5]]},R2],Circle[{xc[[6]],zc[[6]]},R2],Circle[{xc[\
[7]],zc[[7]]},R2],Circle[{xc[[8]],zc[[8]]},R2],Circle[{xc[[9]],zc[[9]]},R4],\
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],Disk[{xn[[6]],zn[[6]\
]},0.08],Disk[{xn[[7]],zn[[7]]},0.08],Disk[{xn[[8]],zn[[8]]},0.08]} ]  /. \
t[[2]]

Show[plot2,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-4,4},{-3,3}}, Frame -> True}]

\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657592523291035`*^9, 3.5657592557923098`*^9}},
 FontSize->16],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.11657442569837306`, 0}, 1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, 1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 0}, 1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, -1.0315734960238365`}, 
    1.0096318340753785`], CircleBox[{0, 0}, -0.5583462051844958], 
   DiskBox[{0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 1.0315734960238372`}, 0.08], 
   DiskBox[{-0.1165744256983734, -1.0315734960238372`}, 0.08], 
   DiskBox[{0.1165744256983734, 0}, 0.08], 
   DiskBox[{0.1165744256983734, 1.0315734960238372`}, 0.08], 
   DiskBox[{0.1165744256983734, -1.0315734960238372`}, 0.08]}]], "Output",
 CellChangeTimes->{3.565759332482044*^9, 3.565759486937916*^9, 
  3.565759814210891*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.11657442569837306`, 0}, 1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, 1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 0}, 1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 1.0315734960238365`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, -1.0315734960238365`}, 
    1.0096318340753785`], CircleBox[{0, 0}, -0.5583462051844958], 
   DiskBox[{0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 1.0315734960238372`}, 0.08], 
   DiskBox[{-0.1165744256983734, -1.0315734960238372`}, 0.08], 
   DiskBox[{0.1165744256983734, 0}, 0.08], 
   DiskBox[{0.1165744256983734, 1.0315734960238372`}, 0.08], 
   DiskBox[{0.1165744256983734, -1.0315734960238372`}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.565759332482044*^9, 3.565759486937916*^9, 
  3.565759814210891*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.11657442569837306`, 1.1911584712364907`}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -0.5955792356182453}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -0.5955792356182453}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, -1.1911584712364907`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 0.5955792356182453}, 1.0096318340753785`],
    CircleBox[{0.11657442569837306`, 0.5955792356182453}, 
    1.0096318340753785`], CircleBox[{0, 0}, -0.5583462051844958], 
   DiskBox[{0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 1.1911584712364915`}, 0.08], 
   DiskBox[{-0.1165744256983734, -0.5955792356182458}, 0.08], 
   DiskBox[{-0.1165744256983734, -0.5955792356182458}, 0.08], 
   DiskBox[{0.1165744256983734, -1.1911584712364915`}, 0.08], 
   DiskBox[{0.1165744256983734, 0.5955792356182458}, 0.08], 
   DiskBox[{0.1165744256983734, 0.5955792356182458}, 0.08]}]], "Output",
 CellChangeTimes->{3.565759332482044*^9, 3.565759486937916*^9, 
  3.565759814210891*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.30456369054118876`, 0}, 0.2537825146433071], 
   CircleBox[{-0.11657442569837306`, 1.1911584712364907`}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -0.5955792356182453}, 
    1.0096318340753785`], 
   CircleBox[{-0.11657442569837306`, -0.5955792356182453}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, -1.1911584712364907`}, 
    1.0096318340753785`], 
   CircleBox[{0.11657442569837306`, 0.5955792356182453}, 1.0096318340753785`],
    CircleBox[{0.11657442569837306`, 0.5955792356182453}, 
    1.0096318340753785`], CircleBox[{0, 0}, -0.5583462051844958], 
   DiskBox[{0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.30456369054118876`, 0}, 0.08], 
   DiskBox[{-0.1165744256983734, 1.1911584712364915`}, 0.08], 
   DiskBox[{-0.1165744256983734, -0.5955792356182458}, 0.08], 
   DiskBox[{-0.1165744256983734, -0.5955792356182458}, 0.08], 
   DiskBox[{0.1165744256983734, -1.1911584712364915`}, 0.08], 
   DiskBox[{0.1165744256983734, 0.5955792356182458}, 0.08], 
   DiskBox[{0.1165744256983734, 0.5955792356182458}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.565759332482044*^9, 3.565759486937916*^9, 
  3.565759814210891*^9}]
}, Open  ]]
},
WindowToolbars->"EditBar",
WindowSize->{1147, 964},
WindowMargins->{{313, Automatic}, {Automatic, 0}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
PrivateNotebookOptions->{"VersionedStylesheet"->{"Default.nb"[8.] -> False}},
CellLabelAutoDelete->True,
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, 647, 17, 185, "Input"],
Cell[1207, 39, 243, 9, 90, "Input"],
Cell[1453, 50, 292, 9, 71, "Input"],
Cell[1748, 61, 1638, 54, 888, "Input"],
Cell[3389, 117, 178, 4, 33, "Input"],
Cell[3570, 123, 839, 27, 318, "Input"],
Cell[CellGroupData[{
Cell[4434, 154, 583, 14, 147, "Input"],
Cell[5020, 170, 719, 13, 39, "Message"],
Cell[5742, 185, 482, 12, 33, "Output"]
}, Open  ]],
Cell[6239, 200, 197, 7, 52, "Input"],
Cell[CellGroupData[{
Cell[6461, 211, 387, 14, 185, "Input"],
Cell[6851, 227, 134, 2, 31, "Output"],
Cell[6988, 231, 113, 1, 31, "Output"],
Cell[7104, 234, 114, 1, 31, "Output"],
Cell[7221, 237, 117, 1, 31, "Output"],
Cell[7341, 240, 131, 2, 31, "Output"],
Cell[7475, 244, 114, 1, 31, "Output"],
Cell[7592, 247, 114, 1, 31, "Output"],
Cell[7709, 250, 113, 1, 31, "Output"],
Cell[7825, 253, 193, 3, 31, "Output"]
}, Open  ]],
Cell[8033, 259, 178, 4, 33, "Input"],
Cell[CellGroupData[{
Cell[8236, 267, 1603, 33, 489, "Input"],
Cell[9842, 302, 1199, 22, 447, "Output"],
Cell[11044, 326, 1311, 27, 293, "Output"],
Cell[12358, 355, 1270, 23, 447, "Output"],
Cell[13631, 380, 1382, 28, 293, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)