(* 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[     10185,        283]
NotebookOptionsPosition[      9202,        249]
NotebookOutlinePosition[      9742,        269]
CellTagsIndexPosition[      9699,        266]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
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;

(* C He-shells *)
Ekin = 2*(2.25*k1/R1^2);
vee=2*(3.0*sig1/R1);
vne=-2*(3.0*z/R1);
(* this is the common vertex assumption *)
R4=R2;

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

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

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.26},{R2,1.24},{p,0.54},{Method -> \"Newton\"}, \
{MaxIterations -> 500}]*)

(* func *)
vne 
vee 
vnn
-Epot/Ekin
2*d1*bohr
(R1+R2+p)*bohr
(R1+R2)*bohr
2*w*rad
4*(R1*(R1+2*R4))^(3/2)/(R1+R4)^3\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657586454256372`*^9, 3.5657586490916443`*^9}, {
  3.5657589497977724`*^9, 3.5657589703274083`*^9}, {3.565759013399084*^9, 
  3.5657590346619215`*^9}},
 FontSize->16],

Cell[BoxData[
 RowBox[{"-", "275.9265225899106`"}]], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539747553`*^9}],

Cell[BoxData["68.32355993618081`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["36.13992668814448`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["2.153628525695231`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["1.596523833938`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["1.0839446081290198`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["0.798261916969`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["109.4712215648118`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.5657590539903555`*^9}],

Cell[BoxData["0.7159310103245544`"], "Output",
 CellChangeTimes->{3.5657586759080915`*^9, 3.5657588697540317`*^9, 
  3.565759054005955*^9}]
}, Open  ]],

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]} ];

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]} ];

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

\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5657587051269426`*^9, 3.565758711039353*^9}, {
  3.565758836869174*^9, 3.5657588421731834`*^9}, {3.565759042742736*^9, 
  3.5657590473291435`*^9}},
 FontSize->16],

Cell[BoxData[
 GraphicsBox[{CircleBox[{-1.508497, 0}, 0.262361], 
   CircleBox[{1.508497, 0}, 0.262361], 
   CircleBox[{-2.0113293333333333`, 0}, 1.246136], 
   CircleBox[{-2.0113293333333333`, -1.2316826428397318`}, 1.246136], 
   CircleBox[{-2.0113293333333333`, 1.2316826428397318`}, 1.246136], 
   CircleBox[{2.0113293333333333`, 0}, 1.246136], 
   CircleBox[{2.0113293333333333`, 1.2316826428397318`}, 1.246136], 
   CircleBox[{2.0113293333333333`, -1.2316826428397318`}, 1.246136], 
   CircleBox[{0, 0}, 1.246136], DiskBox[{-1.508497, 0}, 0.08], 
   DiskBox[{1.508497, 0}, 0.08], DiskBox[{-2.1912834199999995`, 0}, 0.08], 
   DiskBox[{-2.1912834199999995`, 1.6724783323016468`}, 0.08], 
   DiskBox[{-2.1912834199999995`, -1.6724783323016468`}, 0.08], 
   DiskBox[{2.19128342, 0}, 0.08], 
   DiskBox[{2.19128342, 1.6724783323016468`}, 0.08], 
   DiskBox[{2.19128342, -1.6724783323016468`}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.565758676017291*^9, 3.565758869816432*^9, 
  3.565759054037155*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{-1.508497, 0}, 0.262361], 
   CircleBox[{1.508497, 0}, 0.262361], 
   CircleBox[{-2.0113293333333333`, 1.4222246107994179`}, 1.246136], 
   CircleBox[{-2.0113293333333333`, -0.7111123053997089}, 1.246136], 
   CircleBox[{-2.0113293333333333`, -0.7111123053997089}, 1.246136], 
   CircleBox[{2.0113293333333333`, -1.4222246107994179`}, 1.246136], 
   CircleBox[{2.0113293333333333`, 0.7111123053997089}, 1.246136], 
   CircleBox[{2.0113293333333333`, 0.7111123053997089}, 1.246136], 
   CircleBox[{0, 0}, 1.246136], DiskBox[{-1.508497, 0}, 0.08], 
   DiskBox[{1.508497, 0}, 0.08], 
   DiskBox[{-2.1912834199999995`, 1.9312116307363443`}, 0.08], 
   DiskBox[{-2.1912834199999995`, -0.9656058153681721}, 0.08], 
   DiskBox[{-2.1912834199999995`, -0.9656058153681721}, 0.08], 
   DiskBox[{2.19128342, -1.9312116307363443`}, 0.08], 
   DiskBox[{2.19128342, 0.9656058153681721}, 0.08], 
   DiskBox[{2.19128342, 0.9656058153681721}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.565758676017291*^9, 3.565758869816432*^9, 
  3.5657590540527554`*^9}]
}, Open  ]]
},
WindowToolbars->"EditBar",
WindowSize->{1131, 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[CellGroupData[{
Cell[579, 22, 3285, 100, 1553, "Input"],
Cell[3867, 124, 157, 3, 31, "Output"],
Cell[4027, 129, 140, 2, 31, "Output"],
Cell[4170, 133, 140, 2, 31, "Output"],
Cell[4313, 137, 140, 2, 31, "Output"],
Cell[4456, 141, 137, 2, 31, "Output"],
Cell[4596, 145, 141, 2, 31, "Output"],
Cell[4740, 149, 137, 2, 31, "Output"],
Cell[4880, 153, 140, 2, 31, "Output"],
Cell[5023, 157, 139, 2, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[5199, 164, 1679, 33, 489, "Input"],
Cell[6881, 199, 1115, 22, 293, "Output"],
Cell[7999, 223, 1187, 23, 325, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)