(* 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[     32220,        789]
NotebookOptionsPosition[     29564,        705]
NotebookOutlinePosition[     30104,        725]
CellTagsIndexPosition[     30061,        722]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* H2C=CH2 with spherical bananas: Two fused methanes with CH4 data \
CH4Mt(min).nb  20.06.2012 TEST! *)
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;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, 
   3.5668006247245502`*^9}],

Cell["\<\
pi=0.0; (* distance of bananas from C-C axis, pi=0 is regular tetrahedric *)
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.27; k4=1.656; (* parameters for kinetic energy of clouds; \
k=1.0 Kimball's lowest value *)
bohr=0.529177; rad=57.29578;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006303405604`*^9}, {3.5668170547514086`*^9, 
  3.5668170590102158`*^9}, {3.5668172256185083`*^9, 3.566817256584563*^9}, {
  3.566817851008007*^9, 3.566817851554008*^9}, {3.56681792693334*^9, 
  3.5668179459965734`*^9}, {3.566818031469124*^9, 3.566818055025165*^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",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.56680063598777*^9}, {3.566802700978997*^9, 
  3.56680270817061*^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[ArcCos[-1/3]/2];  ss=Sin[ArcCos[-1/3]/2];
(* edge length of tetrahedron of 4 equal clouds *)
a=R4; (*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,-pi-a};
(* potential energy of protons in CH-clouds with eccentricity p *)
vne=vne-4*(3-(p/R2)^2)/R2;\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.566800645940588*^9}, {3.566803037549588*^9, 
  3.5668030837880697`*^9}, {3.5668036983355494`*^9, 3.5668037058079624`*^9}, {
  3.5668042399061003`*^9, 3.566804247784114*^9}, {3.5668167100220027`*^9, 
  3.5668167195380197`*^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]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.566800657812209*^9}}],

Cell["\<\
(* 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]]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006625858173`*^9}}],

Cell["\<\
(* 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.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006669226246`*^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.265},{d1,1.24},{R2,1.20},{R4,1.68},{p,0.54},{w,1.3}\
,{Method -> \"Newton\", MaxIterations -> 800}]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006719614334`*^9}, {3.566802680886162*^9, 
  3.5668026869077725`*^9}, {3.5668027663743124`*^9, 3.5668027861551466`*^9}, {
  3.5668031435049744`*^9, 3.5668031532549915`*^9}, {3.566803746040433*^9, 
  3.5668037659304676`*^9}, {3.5668040820650234`*^9, 3.5668040971502495`*^9}, {
  3.56680426798615*^9, 3.5668042822601748`*^9}, {3.566804358279108*^9, 
  3.5668043613679132`*^9}, {3.5668048168263135`*^9, 3.566804817497115*^9}, {
  3.566817988179048*^9, 3.5668180095978856`*^9}}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "76.17212418843475`"}], ",", 
   RowBox[{"{", 
    RowBox[{
     RowBox[{"R1", "\[Rule]", "0.26521676639954245`"}], ",", 
     RowBox[{"d1", "\[Rule]", "1.2547215686010675`"}], ",", 
     RowBox[{"R2", "\[Rule]", "1.2036067366105956`"}], ",", 
     RowBox[{"R4", "\[Rule]", "1.699600670922859`"}], ",", 
     RowBox[{"p", "\[Rule]", "0.5823516907180097`"}], ",", 
     RowBox[{"w", "\[Rule]", "1.3`"}]}], "}"}]}], "}"}]], "Output",
 CellChangeTimes->{
  3.566800788930439*^9, 3.5668025580671463`*^9, 3.566802714644621*^9, 
   3.566802794204761*^9, 3.5668030940996876`*^9, 3.566803167139016*^9, 
   3.566803799189726*^9, 3.5668041083510695`*^9, 3.566804368528326*^9, 
   3.5668050178326664`*^9, 3.5668167378056517`*^9, 3.566817071084637*^9, {
   3.566817234026923*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933500952*^9, 3.566817951955784*^9}, {3.5668180177411*^9, 
   3.566818060984376*^9}}]
}, Open  ]],

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.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006905878663`*^9}}],

Cell[BoxData[
 RowBox[{"-", "239.92461929387977`"}]], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["54.27116667389774`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["33.309204242379685`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["1.999999999990379`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["1.3279395910152143`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["1.0854347354914802`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["0.7772676148523958`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["148.969028`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}],

Cell[BoxData["0.5052454074208601`"], "Output",
 CellChangeTimes->{
  3.566800788946039*^9, 3.5668025580827465`*^9, 3.566802714660221*^9, 
   3.566802794220361*^9, 3.566803094115288*^9, 3.5668031671546164`*^9, 
   3.566803799205326*^9, 3.5668041083510695`*^9, 3.5668043685439262`*^9, 
   3.5668050178482666`*^9, 3.5668167378056517`*^9, 3.566817071100237*^9, {
   3.5668172340425234`*^9, 3.5668172614985714`*^9}, 3.5668178580904193`*^9, {
   3.566817933516552*^9, 3.566817951955784*^9}, {3.5668180177566996`*^9, 
   3.5668180609999757`*^9}}]
}, Open  ]],

Cell["\<\
(* Cvnf: sum of C1-nuclei forces *)
Cvnf = 0.0;
  For[j = 2, j < 7, j++,
    Cvnf = Cvnf - \
(xn[[j]]-xn[[1]])*ch[[1]]*ch[[j]]/((xn[[1]]-xn[[j]])^2+(yn[[1]]-yn[[j]])^2+(\
zn[[1]]-zn[[j]])^2)^(3/2)]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668006962662764`*^9}}],

Cell[CellGroupData[{

Cell["\<\
(* Cvef: sum of C1 - clouds forces *)
Cvef = 0.0;
  For[j = 2, j < nc+1, j++,
    Cvef = Cvef - \
(xc[[j]]-xn[[1]])*oc[[j]]*ch[[1]]/((xc[[1]]-xn[[j]])^2+(yc[[1]]-yn[[j]])^2+(\
zc[[1]]-zn[[j]])^2)^(3/2)]
Cforcetot = Cvnf+Cvef /. t[[2]]
Cforcenn=Cvnf /. t[[2]]
CForcene=Cvef /. t[[2]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.566800704503091*^9}}],

Cell[BoxData["11.265281864281395`"], "Output",
 CellChangeTimes->{
  3.566800788961639*^9, 3.5668025580983467`*^9, 3.5668027146758213`*^9, 
   3.566802794251561*^9, 3.566803094146488*^9, 3.566803167170216*^9, 
   3.5668037992209263`*^9, 3.566804108366669*^9, 3.5668043685595264`*^9, 
   3.566805017879467*^9, 3.566816737821252*^9, 3.566817071115837*^9, {
   3.566817234058123*^9, 3.566817261529772*^9}, 3.566817858121619*^9, {
   3.566817933532152*^9, 3.5668179519869843`*^9}, {3.5668180177879*^9, 
   3.566818061015576*^9}}],

Cell[BoxData[
 RowBox[{"-", "4.734480287054216`"}]], "Output",
 CellChangeTimes->{
  3.566800788961639*^9, 3.5668025580983467`*^9, 3.5668027146758213`*^9, 
   3.566802794251561*^9, 3.566803094146488*^9, 3.566803167170216*^9, 
   3.5668037992209263`*^9, 3.566804108366669*^9, 3.5668043685595264`*^9, 
   3.566805017879467*^9, 3.566816737821252*^9, 3.566817071115837*^9, {
   3.566817234058123*^9, 3.566817261529772*^9}, 3.566817858121619*^9, {
   3.566817933532152*^9, 3.5668179519869843`*^9}, {3.5668180177879*^9, 
   3.566818061015576*^9}}],

Cell[BoxData["15.999762151335615`"], "Output",
 CellChangeTimes->{
  3.566800788961639*^9, 3.5668025580983467`*^9, 3.5668027146758213`*^9, 
   3.566802794251561*^9, 3.566803094146488*^9, 3.566803167170216*^9, 
   3.5668037992209263`*^9, 3.566804108366669*^9, 3.5668043685595264`*^9, 
   3.566805017879467*^9, 3.566816737821252*^9, 3.566817071115837*^9, {
   3.566817234058123*^9, 3.566817261529772*^9}, 3.566817858121619*^9, {
   3.566817933532152*^9, 3.5668179519869843`*^9}, {3.5668180177879*^9, 
   3.566818061015576*^9}}]
}, Open  ]],

Cell["\<\
(* Hvnfx: sum of H3-nuclei forces *)
Hvnfx = 0.0;
  For[j = 1, j < 7, j++,
    If[j != 3,
    Hvnfx = Hvnfx - \
(xn[[j]]-xn[[3]])*ch[[3]]*ch[[j]]/((xn[[3]]-xn[[j]])^2+(yn[[3]]-yn[[j]])^2+(\
zn[[3]]-zn[[j]])^2)^(3/2)]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668007097759*^9}}],

Cell[CellGroupData[{

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

Cell[BoxData[
 RowBox[{"-", "0.4087162158475007`"}]], "Output",
 CellChangeTimes->{
  3.566800788977239*^9, 3.5668025581139464`*^9, 3.5668027146914215`*^9, 
   3.5668027942671614`*^9, 3.566803094162088*^9, 3.5668031671858163`*^9, 
   3.5668037992365265`*^9, 3.566804108382269*^9, 3.5668043685751266`*^9, 
   3.5668050178950667`*^9, 3.5668167378368516`*^9, 3.5668170711314373`*^9, {
   3.5668172340737233`*^9, 3.5668172615453715`*^9}, 3.5668178581372194`*^9, {
   3.566817933547752*^9, 3.5668179520025845`*^9}, {3.5668180178035*^9, 
   3.5668180610311756`*^9}}],

Cell[BoxData[
 RowBox[{"-", "1.2211476306870983`"}]], "Output",
 CellChangeTimes->{
  3.566800788977239*^9, 3.5668025581139464`*^9, 3.5668027146914215`*^9, 
   3.5668027942671614`*^9, 3.566803094162088*^9, 3.5668031671858163`*^9, 
   3.5668037992365265`*^9, 3.566804108382269*^9, 3.5668043685751266`*^9, 
   3.5668050178950667`*^9, 3.5668167378368516`*^9, 3.5668170711314373`*^9, {
   3.5668172340737233`*^9, 3.5668172615453715`*^9}, 3.5668178581372194`*^9, {
   3.566817933547752*^9, 3.5668179520025845`*^9}, {3.5668180178035*^9, 
   3.5668180610311756`*^9}}],

Cell[BoxData["0.8124314148395975`"], "Output",
 CellChangeTimes->{
  3.566800788977239*^9, 3.5668025581139464`*^9, 3.5668027146914215`*^9, 
   3.5668027942671614`*^9, 3.566803094162088*^9, 3.5668031671858163`*^9, 
   3.5668037992365265`*^9, 3.566804108382269*^9, 3.5668043685751266`*^9, 
   3.5668050178950667`*^9, 3.5668167378368516`*^9, 3.5668170711314373`*^9, {
   3.5668172340737233`*^9, 3.5668172615453715`*^9}, 3.5668178581372194`*^9, {
   3.566817933547752*^9, 3.5668179520025845`*^9}, {3.5668180178035*^9, 
   3.5668180610311756`*^9}}]
}, Open  ]],

Cell["\<\
(* Hvnfy: sum of H3-nuclei forces *)
Hvnfy = 0.0;
  For[j = 1, j < 7, j++,
    If[j != 3,
    Hvnfy = Hvnfy - \
(yn[[j]]-yn[[3]])*ch[[3]]*ch[[j]]/((xn[[3]]-xn[[j]])^2+(yn[[3]]-yn[[j]])^2+(\
zn[[3]]-zn[[j]])^2)^(3/2)]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.566800727622331*^9}}],

Cell[CellGroupData[{

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

Cell[BoxData["0.5290953932405441`"], "Output",
 CellChangeTimes->{
  3.566800788992839*^9, 3.566802558129546*^9, 3.5668027147070217`*^9, 
   3.566802794282761*^9, 3.5668030941776876`*^9, 3.566803167217016*^9, 
   3.5668037992365265`*^9, 3.566804108397869*^9, 3.5668043685907264`*^9, 
   3.5668050179106665`*^9, 3.566816737852452*^9, 3.5668170711470375`*^9, {
   3.566817234089323*^9, 3.5668172615609713`*^9}, 3.5668178581528196`*^9, {
   3.566817933563352*^9, 3.5668179520181847`*^9}, {3.5668180178191*^9, 
   3.566818061046776*^9}}],

Cell[BoxData["1.4203281382404769`"], "Output",
 CellChangeTimes->{
  3.566800788992839*^9, 3.566802558129546*^9, 3.5668027147070217`*^9, 
   3.566802794282761*^9, 3.5668030941776876`*^9, 3.566803167217016*^9, 
   3.5668037992365265`*^9, 3.566804108397869*^9, 3.5668043685907264`*^9, 
   3.5668050179106665`*^9, 3.566816737852452*^9, 3.5668170711470375`*^9, {
   3.566817234089323*^9, 3.5668172615609713`*^9}, 3.5668178581528196`*^9, {
   3.566817933563352*^9, 3.5668179520181847`*^9}, {3.5668180178191*^9, 
   3.566818061046776*^9}}],

Cell[BoxData[
 RowBox[{"-", "0.891232744999933`"}]], "Output",
 CellChangeTimes->{
  3.566800788992839*^9, 3.566802558129546*^9, 3.5668027147070217`*^9, 
   3.566802794282761*^9, 3.5668030941776876`*^9, 3.566803167217016*^9, 
   3.5668037992365265`*^9, 3.566804108397869*^9, 3.5668043685907264`*^9, 
   3.5668050179106665`*^9, 3.566816737852452*^9, 3.5668170711470375`*^9, {
   3.566817234089323*^9, 3.5668172615609713`*^9}, 3.5668178581528196`*^9, {
   3.566817933563352*^9, 3.5668179520181847`*^9}, {3.5668180178191*^9, 
   3.566818061046776*^9}}]
}, Open  ]],

Cell["\<\
(* Hvnfz: sum of H3-nuclei forces *)
Hvnfz = 0.0;
  For[j = 7, j < 9, j++,
    If[j != 3,
    Hvnfz = Hvnfz - \
(zn[[j]]-zn[[3]])*ch[[3]]*ch[[j]]/((xn[[3]]-xn[[j]])^2+(yn[[3]]-yn[[j]])^2+(\
zn[[3]]-zn[[j]])^2)^(3/2)]]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668007404299536`*^9}, {3.5668020045937743`*^9, 
  3.566802043765443*^9}, {3.566802206068128*^9, 3.566802232900175*^9}, {
  3.5668022898090754`*^9, 3.566802295066284*^9}, {3.56680237223962*^9, 
  3.56680239537446*^9}}],

Cell[CellGroupData[{

Cell["\<\
(* Hvefz: sum of H3 - clouds forces *)
Hvefz = 0.0;
  For[j = 7, j < nc+1, j++,
    If[j != 3,
    Hvefz = Hvefz - \
(zc[[j]]-zn[[3]])*oc[[j]]*ch[[3]]/((xc[[j]]-xn[[3]])^2+(yc[[j]]-yn[[3]])^2+(\
zc[[j]]-zn[[3]])^2)^(3/2)]]
Hforcetotz = Hvnfz+Hvefz /. t[[2]]
Hforcennz=Hvnfz /. t[[2]]
HForcenez=Hvefz /. t[[2]]        \
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668007404299536`*^9}, {3.5668020045937743`*^9, 
  3.566802043765443*^9}, {3.566802206068128*^9, 3.566802232900175*^9}, {
  3.5668022898090754`*^9, 3.566802295066284*^9}, {3.56680237223962*^9, 
  3.566802503638651*^9}}],

Cell[BoxData["0.`"], "Output",
 CellChangeTimes->{
  3.5668025581607466`*^9, 3.5668027147226214`*^9, 3.5668027942983613`*^9, 
   3.566803094193288*^9, 3.5668031672326164`*^9, 3.5668037992521267`*^9, 
   3.5668041084134693`*^9, 3.5668043686219263`*^9, 3.566805017941867*^9, 
   3.5668167378680515`*^9, 3.566817071162637*^9, {3.566817234089323*^9, 
   3.5668172615765715`*^9}, 3.56681785816842*^9, {3.566817933578952*^9, 
   3.5668179520337844`*^9}, {3.5668180178347*^9, 3.566818061046776*^9}}],

Cell[BoxData["0.`"], "Output",
 CellChangeTimes->{
  3.5668025581607466`*^9, 3.5668027147226214`*^9, 3.5668027942983613`*^9, 
   3.566803094193288*^9, 3.5668031672326164`*^9, 3.5668037992521267`*^9, 
   3.5668041084134693`*^9, 3.5668043686219263`*^9, 3.566805017941867*^9, 
   3.5668167378680515`*^9, 3.566817071162637*^9, {3.566817234089323*^9, 
   3.5668172615765715`*^9}, 3.56681785816842*^9, {3.566817933578952*^9, 
   3.5668179520337844`*^9}, {3.5668180178347*^9, 3.566818061062376*^9}}],

Cell[BoxData["0.`"], "Output",
 CellChangeTimes->{
  3.5668025581607466`*^9, 3.5668027147226214`*^9, 3.5668027942983613`*^9, 
   3.566803094193288*^9, 3.5668031672326164`*^9, 3.5668037992521267`*^9, 
   3.5668041084134693`*^9, 3.5668043686219263`*^9, 3.566805017941867*^9, 
   3.5668167378680515`*^9, 3.566817071162637*^9, {3.566817234089323*^9, 
   3.5668172615765715`*^9}, 3.56681785816842*^9, {3.566817933578952*^9, 
   3.5668179520337844`*^9}, {3.5668180178347*^9, 3.566818061062376*^9}}]
}, Open  ]],

Cell["\<\
(* Wvnfy: sum of w7-nuclei forces *)
Wvnfy = 0.0;
  For[j = 1, j < 7, j++,
    Wvnfy = Wvnfy - \
(yn[[j]]-zc[[7]])*oc[[7]]*ch[[j]]/((xc[[7]]-xn[[j]])^2+(yc[[7]]-yn[[j]])^2+(\
zc[[7]]-zn[[j]])^2)^(3/2)]\
\>", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.566799684588899*^9, 3.56679969039211*^9}, {
  3.5668006247245502`*^9, 3.5668007404299536`*^9}, {3.5668020045937743`*^9, 
  3.566802043765443*^9}, {3.566802512827067*^9, 3.5668025183962765`*^9}}],

Cell[CellGroupData[{

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

Cell[BoxData[
 RowBox[{"-", "2.0025182975471365`"}]], "Output",
 CellChangeTimes->{
  3.566800789008439*^9, 3.5668025581763463`*^9, 3.5668027147382216`*^9, 
   3.566802794329561*^9, 3.566803094208888*^9, 3.5668031672482166`*^9, 
   3.5668037992677264`*^9, 3.5668041084290695`*^9, 3.5668043686375265`*^9, 
   3.566805017957467*^9, 3.5668167378836517`*^9, 3.566817071178237*^9, {
   3.5668172341049232`*^9, 3.5668172615921717`*^9}, 3.5668178581840196`*^9, {
   3.566817933594552*^9, 3.566817952049384*^9}, {3.5668180178503*^9, 
   3.566818061062376*^9}}],

Cell[BoxData[
 RowBox[{"-", "4.668644408233097`"}]], "Output",
 CellChangeTimes->{
  3.566800789008439*^9, 3.5668025581763463`*^9, 3.5668027147382216`*^9, 
   3.566802794329561*^9, 3.566803094208888*^9, 3.5668031672482166`*^9, 
   3.5668037992677264`*^9, 3.5668041084290695`*^9, 3.5668043686375265`*^9, 
   3.566805017957467*^9, 3.5668167378836517`*^9, 3.566817071178237*^9, {
   3.5668172341049232`*^9, 3.5668172615921717`*^9}, 3.5668178581840196`*^9, {
   3.566817933594552*^9, 3.566817952049384*^9}, {3.5668180178503*^9, 
   3.566818061077976*^9}}],

Cell[BoxData["2.6661261106859606`"], "Output",
 CellChangeTimes->{
  3.566800789008439*^9, 3.5668025581763463`*^9, 3.5668027147382216`*^9, 
   3.566802794329561*^9, 3.566803094208888*^9, 3.5668031672482166`*^9, 
   3.5668037992677264`*^9, 3.5668041084290695`*^9, 3.5668043686375265`*^9, 
   3.566805017957467*^9, 3.5668167378836517`*^9, 3.566817071178237*^9, {
   3.5668172341049232`*^9, 3.5668172615921717`*^9}, 3.5668178581840196`*^9, {
   3.566817933594552*^9, 3.566817952049384*^9}, {3.5668180178503*^9, 
   3.566818061077976*^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]]},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]];

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

Cell[BoxData[
 GraphicsBox[{CircleBox[{-1.2547215686010675`, 0}, 0.26521676639954245`], 
   CircleBox[{1.2547215686010675`, 0}, 0.26521676639954245`], 
   CircleBox[{-2.10274721345602, 1.199289368194063}, 1.2036067366105956`], 
   CircleBox[{-2.10274721345602, -1.199289368194063}, 1.2036067366105956`], 
   CircleBox[{2.10274721345602, -1.199289368194063}, 1.2036067366105956`], 
   CircleBox[{2.10274721345602, 1.199289368194063}, 1.2036067366105956`], 
   CircleBox[{0, 0}, 1.699600670922859], CircleBox[{0, 0}, 1.699600670922859],
    DiskBox[{-1.2547215686010675`, 0}, 0.08], 
   DiskBox[{1.2547215686010675`, 0}, 0.08], 
   DiskBox[{-2.4389681188550965`, 1.674777532562799}, 0.08], 
   DiskBox[{-2.4389681188550965`, -1.674777532562799}, 0.08], 
   DiskBox[{2.4389681188550965`, -1.674777532562799}, 0.08], 
   DiskBox[{2.4389681188550965`, 1.674777532562799}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{
  3.5668007890240393`*^9, 3.5668025581919465`*^9, 3.5668027147538214`*^9, 
   3.5668027943451614`*^9, 3.5668030942244883`*^9, 3.5668031672638164`*^9, 
   3.566803799283326*^9, 3.5668041084446697`*^9, 3.5668043686375265`*^9, 
   3.566805017957467*^9, 3.566816737899252*^9, 3.566817071178237*^9, {
   3.5668172341205235`*^9, 3.566817261607772*^9}, 3.5668178581996193`*^9, {
   3.5668179336101522`*^9, 3.5668179520649843`*^9}, {3.5668180178659*^9, 
   3.566818061077976*^9}}]
}, Open  ]],

Cell[CellGroupData[{

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

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

Cell[BoxData[
 GraphicsBox[{CircleBox[{-1.2547215686010675`, 0}, 0.26521676639954245`], 
   CircleBox[{1.2547215686010675`, 0}, 0.26521676639954245`], 
   CircleBox[{-2.10274721345602, 0}, 1.2036067366105956`], 
   CircleBox[{-2.10274721345602, 0}, 1.2036067366105956`], 
   CircleBox[{2.10274721345602, 0}, 1.2036067366105956`], 
   CircleBox[{2.10274721345602, 0}, 1.2036067366105956`], 
   CircleBox[{0, 1.699600670922859}, 1.699600670922859], 
   CircleBox[{0, -1.699600670922859}, 1.699600670922859], 
   DiskBox[{-1.2547215686010675`, 0}, 0.08], 
   DiskBox[{1.2547215686010675`, 0}, 0.08], 
   DiskBox[{-2.4389681188550965`, 0}, 0.08], 
   DiskBox[{-2.4389681188550965`, 0}, 0.08], 
   DiskBox[{2.4389681188550965`, 0}, 0.08], 
   DiskBox[{2.4389681188550965`, 0}, 0.08]},
  AspectRatio->Automatic,
  Axes->True,
  Frame->True,
  GridLines->Automatic,
  PlotRange->{{-4, 4}, {-3, 3}}]], "Output",
 CellChangeTimes->{
  3.5668007890396395`*^9, 3.5668025581919465`*^9, 3.5668027147694216`*^9, 
   3.5668027943451614`*^9, 3.566803094240088*^9, 3.566803167279416*^9, 
   3.5668037992989264`*^9, 3.5668041084602695`*^9, 3.5668043686531267`*^9, 
   3.566805017973067*^9, 3.566816737899252*^9, 3.566817071193837*^9, {
   3.5668172341361237`*^9, 3.566817261607772*^9}, 3.5668178581996193`*^9, {
   3.566817933625752*^9, 3.5668179520805845`*^9}, {3.5668180178815002`*^9, 
   3.566818061093576*^9}}]
}, Open  ]],

Cell["", "Input",
 PageWidth->WindowWidth,
 CellChangeTimes->{{3.5668007580735846`*^9, 3.5668007629251933`*^9}}]
},
WindowToolbars->"EditBar",
WindowSize->{843, 964},
WindowMargins->{{Automatic, 129}, {-7, Automatic}},
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, 379, 9, 99, "Input"],
Cell[939, 31, 756, 14, 133, "Input"],
Cell[1698, 47, 347, 11, 116, "Input"],
Cell[2048, 60, 1541, 37, 490, "Input"],
Cell[3592, 99, 371, 10, 99, "Input"],
Cell[3966, 111, 366, 10, 99, "Input"],
Cell[4335, 123, 392, 11, 116, "Input"],
Cell[CellGroupData[{
Cell[4752, 138, 967, 18, 133, "Input"],
Cell[5722, 158, 970, 19, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[6729, 182, 369, 13, 167, "Input"],
Cell[7101, 197, 556, 9, 31, "Output"],
Cell[7660, 208, 538, 8, 31, "Output"],
Cell[8201, 218, 539, 8, 31, "Output"],
Cell[8743, 228, 538, 8, 31, "Output"],
Cell[9284, 238, 539, 8, 31, "Output"],
Cell[9826, 248, 539, 8, 31, "Output"],
Cell[10368, 258, 539, 8, 31, "Output"],
Cell[10910, 268, 531, 8, 31, "Output"],
Cell[11444, 278, 539, 8, 31, "Output"]
}, Open  ]],
Cell[11998, 289, 365, 10, 99, "Input"],
Cell[CellGroupData[{
Cell[12388, 303, 448, 13, 150, "Input"],
Cell[12839, 318, 525, 8, 31, "Output"],
Cell[13367, 328, 541, 9, 31, "Output"],
Cell[13911, 339, 525, 8, 31, "Output"]
}, Open  ]],
Cell[14451, 350, 381, 11, 116, "Input"],
Cell[CellGroupData[{
Cell[14857, 365, 473, 14, 167, "Input"],
Cell[15333, 381, 560, 9, 31, "Output"],
Cell[15896, 392, 560, 9, 31, "Output"],
Cell[16459, 403, 543, 8, 31, "Output"]
}, Open  ]],
Cell[17017, 414, 383, 11, 116, "Input"],
Cell[CellGroupData[{
Cell[17425, 429, 477, 14, 167, "Input"],
Cell[17905, 445, 533, 8, 31, "Output"],
Cell[18441, 455, 533, 8, 31, "Output"],
Cell[18977, 465, 549, 9, 31, "Output"]
}, Open  ]],
Cell[19541, 477, 580, 14, 116, "Input"],
Cell[CellGroupData[{
Cell[20146, 495, 681, 17, 167, "Input"],
Cell[20830, 514, 492, 7, 31, "Output"],
Cell[21325, 523, 492, 7, 31, "Output"],
Cell[21820, 532, 492, 7, 31, "Output"]
}, Open  ]],
Cell[22327, 542, 468, 11, 99, "Input"],
Cell[CellGroupData[{
Cell[22820, 557, 475, 14, 167, "Input"],
Cell[23298, 573, 552, 9, 31, "Output"],
Cell[23853, 584, 551, 9, 31, "Output"],
Cell[24407, 595, 535, 8, 31, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[24979, 608, 751, 14, 201, "Input"],
Cell[25733, 624, 1491, 26, 293, "Output"]
}, Open  ]],
Cell[CellGroupData[{
Cell[27261, 655, 773, 14, 201, "Input"],
Cell[28037, 671, 1396, 27, 293, "Output"]
}, Open  ]],
Cell[29448, 701, 112, 2, 31, "Input"]
}
]
*)

(* End of internal cache information *)