(* 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[     14507,        400]
NotebookOptionsPosition[     13019,        350]
NotebookOutlinePosition[     13532,        369]
CellTagsIndexPosition[     13489,        366]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* H2C=CH2 with spherical bananas: Two fused methanes with CH4 data \
CH4Mt(min).nb  17.06.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;
pi=0.0; (* distance of bananas from C-C axis, pi =0 is regular tetraedric *)
nc=8;   (* 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",
 CellChangeTimes->{3.566705957935094*^9}],

Cell["\<\
(* 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 edge assumption *)
R4=R2;\
\>", "Input",
 CellChangeTimes->{{3.566705957935094*^9, 3.5667059682779117`*^9}}],

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

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

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

(* cloud radii in the same order *)
rr={R1,R1,R2,R2,R2,R2,R4,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, xy plane of molecule *)
(* nuclear coordinates in terms of radii; C nucleus assumed in center of \
C(1s) cloud *)
d1=(R1+R2)*cs;  (* R1+R2 is radius of outer sphere for tetrahedron of equal \
clouds *)
d2=d1+(R1+R2+p)*cs; d3=(R1+R2+p)*ss; 
xn={-d1,d1,-d2,-d2,d2,d2,0,0};
yn={0,0,d3,-d3,-d3,d3,0,0};
zn={0,0,0,0,0,0,0,0};

(* cloud coordinates in terms of radii *)        
d4=d1+(R1+R2)*cs;  d5=(R1+R2)*ss;
xc={-d1,d1,-d4,-d4,d4,d4,0,0};
yc={0,0,d5,-d5,-d5,d5,0,0};
zc={0,0,0,0,0,0,pi+a/2,-pi-a/2};

(* potential energy of protons in CH-clouds with eccentricity p *)
vne=vne-4*(3-(p/R2)^2)/R2;\
\>", "Input",
 CellChangeTimes->{{3.566705957935094*^9, 3.5667059817719355`*^9}, {
  3.56670638868305*^9, 3.5667064556383677`*^9}}],

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

(* 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",
 CellChangeTimes->{{3.566705957935094*^9, 3.566705993019555*^9}}],

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.26},{R2,1.24},{p,0.54},{Method -> Automatic}, \
{MaxIterations -> 500}]\
\>", "Input",
 CellChangeTimes->{{3.566705957935094*^9, 3.566706016871997*^9}}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "77.46830379455953`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.26273679950137213`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.2563671190699615`"}], ",", 
     RowBox[{"p", "\[Rule]", "0.5812385024419693`"}]}], "}"}]}], 
  "}"}]], "Output",
 CellChangeTimes->{3.5667063204953303`*^9, 3.566707142413974*^9, 
  3.566707239258944*^9, 3.566707308398266*^9, 3.5667073673663692`*^9}]
}, Open  ]],

Cell[CellGroupData[{

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

Cell[BoxData[
 RowBox[{"-", "257.1960678700093`"}]], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["61.87420023182061`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["40.385260070844495`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["2.000000000281081`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["0.9282347270703322`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["1.1114529013245567`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}],

Cell[BoxData["109.4712215648118`"], "Output",
 CellChangeTimes->{3.56670632051093*^9, 3.566707142429574*^9, 
  3.566707239274544*^9, 3.5667073084138656`*^9, 3.5667073673819695`*^9}]
}, Open  ]],

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]]},R4],Circle[{xc[[8]],yc[[8]]},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]} ] /. t[[2]];
\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{
  3.5667060446088457`*^9, {3.566706282602864*^9, 3.5667062931484823`*^9}}],

Cell[CellGroupData[{

Cell["\<\
Show[plot1,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-4,4},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{3.5667060446088457`*^9, 3.566706282602864*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{-0.8770550563141747, 0}, 0.26273679950137213`], 
   CircleBox[{0.8770550563141747, 0}, 0.26273679950137213`], 
   CircleBox[{-1.7541101126283494`, 1.2403431555874045`}, 
    1.2563671190699615`], 
   CircleBox[{-1.7541101126283494`, -1.2403431555874045`}, 
    1.2563671190699615`], 
   CircleBox[{1.7541101126283494`, -1.2403431555874045`}, 
    1.2563671190699615`], 
   CircleBox[{1.7541101126283494`, 1.2403431555874045`}, 1.2563671190699615`],
    CircleBox[{0, 0}, 1.2563671190699615`], 
   CircleBox[{0, 0}, 1.2563671190699615`], 
   DiskBox[{-0.8770550563141747, 0}, 0.08], 
   DiskBox[{0.8770550563141747, 0}, 0.08], 
   DiskBox[{-2.089688318476595, 1.7149224055348242`}, 0.08], 
   DiskBox[{-2.089688318476595, -1.7149224055348242`}, 0.08], 
   DiskBox[{2.089688318476595, -1.7149224055348242`}, 0.08], 
   DiskBox[{2.089688318476595, 1.7149224055348242`}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.5667063205265303`*^9, 3.566707142445174*^9, 
  3.566707239290144*^9, 3.5667073084294653`*^9, 3.566707367397569*^9}]
}, Open  ]],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{3.5667060446088457`*^9, 3.566706282602864*^9}],

Cell["\<\
(* 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]]},R4],Circle[{xc[[8]],zc[[8]]},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]} ] /. t[[2]];
\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5667060446088457`*^9, 3.566706052174859*^9}, {
   3.566706258844022*^9, 3.5667062634928303`*^9}, 3.566706299138893*^9}],

Cell[CellGroupData[{

Cell["\<\
Show[plot2,{AspectRatio \[Rule] Automatic,Axes -> True,GridLines -> \
Automatic, PlotRange \[Rule] {{-4,4},{-3,3}}, Frame -> True}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5667060446088457`*^9, 3.566706052174859*^9}, 
   3.566706258844022*^9}],

Cell[BoxData[
 GraphicsBox[{CircleBox[{-0.8770550563141747, 0}, 0.26273679950137213`], 
   CircleBox[{0.8770550563141747, 0}, 0.26273679950137213`], 
   CircleBox[{-1.7541101126283494`, 0}, 1.2563671190699615`], 
   CircleBox[{-1.7541101126283494`, 0}, 1.2563671190699615`], 
   CircleBox[{1.7541101126283494`, 0}, 1.2563671190699615`], 
   CircleBox[{1.7541101126283494`, 0}, 1.2563671190699615`], 
   CircleBox[{0, 1.2403431555874045`}, 1.2563671190699615`], 
   CircleBox[{0, -1.2403431555874045`}, 1.2563671190699615`], 
   DiskBox[{-0.8770550563141747, 0}, 0.08], 
   DiskBox[{0.8770550563141747, 0}, 0.08], 
   DiskBox[{-2.089688318476595, 0}, 0.08], 
   DiskBox[{-2.089688318476595, 0}, 0.08], 
   DiskBox[{2.089688318476595, 0}, 0.08], 
   DiskBox[{2.089688318476595, 0}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{3.5667063205265303`*^9, 3.5667071424763737`*^9, 
  3.566707239305744*^9, 3.5667073084450655`*^9, 3.566707367428769*^9}]
}, Open  ]],

Cell["\<\
(* 3D of molecule *)
plot3=Graphics3D[{Sphere[{xc[[1]],yc[[1]],0},R1],      \
Sphere[{xc[[2]],yc[[2]],0},R1],{Opacity[0.4],Sphere[{xc[[3]],yc[[3]],0},R2],\
Sphere[{xc[[4]],yc[[4]],0},R2],Sphere[{xc[[5]],yc[[5]],0},R2],Sphere[{xc[[6]],\
yc[[6]],0},R2]},Sphere[{xc[[7]],yc[[7]],zc[[7]]},R4],Sphere[{xc[[8]],yc[[8]],\
zc[[8]]},R4],Sphere[{xn[[1]],yn[[1]],0},0.08], \
Sphere[{xn[[2]],yn[[2]],0},0.08],Sphere[{xn[[3]],yn[[3]],0},0.08],Sphere[{xn[[\
4]],yn[[4]],0},0.08],
Sphere[{xn[[5]],yn[[5]],0},0.08],Sphere[{xn[[6]],yn[[6]],0},0.08]} ] /. t[[2]];
\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5667060446088457`*^9, 3.566706052174859*^9}, 
   3.566706258844022*^9, 3.5667066813395643`*^9, {3.566706717126027*^9, 
   3.566706723568838*^9}, {3.5667068030509777`*^9, 3.5667068293838243`*^9}, {
   3.566706865763088*^9, 3.566706888554728*^9}, {3.566707002668928*^9, 
   3.5667071319151554`*^9}, {3.566707194814466*^9, 3.566707232784933*^9}, {
   3.5667072976654468`*^9, 3.566707303577857*^9}, {3.5667073620155597`*^9, 
   3.5667073629827614`*^9}}],

Cell[CellGroupData[{

Cell["\<\
Show[plot3,{AspectRatio \[Rule] Automatic, Boxed -> False}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5667060446088457`*^9, 3.566706052174859*^9}, {
  3.566706701385599*^9, 3.5667067467972794`*^9}}],

Cell[BoxData[
 Graphics3DBox[{SphereBox[{-0.8770550563141747, 0, 0}, 0.26273679950137213`], 
   SphereBox[{0.8770550563141747, 0, 0}, 0.26273679950137213`], 
   {Opacity[0.4], 
    SphereBox[{-1.7541101126283494`, 1.2403431555874045`, 0}, 
     1.2563671190699615`], 
    SphereBox[{-1.7541101126283494`, -1.2403431555874045`, 0}, 
     1.2563671190699615`], 
    SphereBox[{1.7541101126283494`, -1.2403431555874045`, 0}, 
     1.2563671190699615`], 
    SphereBox[{1.7541101126283494`, 1.2403431555874045`, 0}, 
     1.2563671190699615`]}, 
   SphereBox[{0, 0, 1.2403431555874045`}, 1.2563671190699615`], 
   SphereBox[{0, 0, -1.2403431555874045`}, 1.2563671190699615`], 
   SphereBox[{-0.8770550563141747, 0, 0}, 0.08], 
   SphereBox[{0.8770550563141747, 0, 0}, 0.08], 
   SphereBox[{-2.089688318476595, 1.7149224055348242`, 0}, 0.08], 
   SphereBox[{-2.089688318476595, -1.7149224055348242`, 0}, 0.08], 
   SphereBox[{2.089688318476595, -1.7149224055348242`, 0}, 0.08], 
   SphereBox[{2.089688318476595, 1.7149224055348242`, 0}, 0.08]},
  AspectRatio->Automatic,
  AutomaticImageSize->True,
  Boxed->False,
  ImageSize->{350.13976243111085`, 355.95825889070466`},
  ViewPoint->{0.40543748195694873`, -2.9198445666533663`, 1.6613633422009855`},
  ViewVertical->{0.006552929448527978, -0.14364686199059146`, 
   1.1971644331112477`}]], "Output",
 CellChangeTimes->{3.5667071427727747`*^9, 3.566707239336944*^9, 
  3.566707308476266*^9, 3.5667073674443693`*^9}]
}, Open  ]]
},
WindowToolbars->"EditBar",
WindowSize->{911, 964},
WindowMargins->{{493, Automatic}, {-5, Automatic}},
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, 637, 14, 167, "Input"],
Cell[1197, 36, 223, 8, 116, "Input"],
Cell[1423, 46, 1280, 41, 626, "Input"],
Cell[2706, 89, 720, 22, 286, "Input"],
Cell[CellGroupData[{
Cell[3451, 115, 407, 10, 116, "Input"],
Cell[3861, 127, 461, 11, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[4359, 143, 238, 10, 150, "Input"],
Cell[4600, 155, 198, 3, 31, "Output"],
Cell[4801, 160, 181, 2, 31, "Output"],
Cell[4985, 164, 182, 2, 31, "Output"],
Cell[5170, 168, 181, 2, 31, "Output"],
Cell[5354, 172, 182, 2, 31, "Output"],
Cell[5539, 176, 182, 2, 31, "Output"],
Cell[5724, 180, 181, 2, 31, "Output"]
}, Open  ]],
Cell[5920, 185, 643, 13, 133, "Input"],
Cell[CellGroupData[{
Cell[6588, 202, 247, 5, 48, "Input"],
Cell[6838, 209, 1174, 24, 293, "Output"]
}, Open  ]],
Cell[8027, 236, 108, 2, 31, "Input"],
Cell[8138, 240, 712, 13, 133, "Input"],
Cell[CellGroupData[{
Cell[8875, 257, 275, 6, 48, "Input"],
Cell[9153, 265, 1050, 21, 293, "Output"]
}, Open  ]],
Cell[10218, 289, 1068, 19, 150, "Input"],
Cell[CellGroupData[{
Cell[11311, 312, 228, 5, 31, "Input"],
Cell[11542, 319, 1461, 28, 403, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)