GeometryOptimizations - Maple Help

Home : Support : Online Help : Toolboxes : Quantum Chemistry : Courses : Lessons : GeometryOptimizations

Geometry Optimization and Normal Modes

 Overview In this exercise, you will determine the optimum geometry for gaseous hypochlorous acid, HOCl. Namely, you will determine the two bond lengths and one bond angle that minimize the potential energy due to all nuclear and electronic interactions.   Furthermore, you will determine the vibrational frequencies of HOCl using a normal mode analysis, which treats the potential surface in the vicinity of the minimum as being harmonic and calculates the corresponding 3N - 6 (for non-linear molecule) or 3N - 5 (for a linear molecule) vibrational modes, each corresponding to a concerted motion of all atoms.  For a non-linear triatomic, this corresponds to a symmetric and antisymmetric stretch and a bending motion. For a linear triatomic, this corresponds to symmetric and antisymmetric stretch and two degenerate bending motions.  Once you have calculated the geometry and vibrational frequencies, they can be compared with experimental measurements.   While computational chemistry can be a powerful tool, the quality of results depends on the accuracy of the electronic structure method employed and the atomic orbital basis used.   In this activity, you can compare geometries and vibrational frequencies calculated using Hartree-Fock, Density Functional Theory, or any method supported by the QuantumChemistry package.  Note that geometry optimizations for more sophisticated methods can be computationally expensive!



Initialize

Here you can choose the electronic structure method and atomic orbital basis for the geometry optimization and frequency calculation.

 > $\mathrm{with}\left(\mathrm{QuantumChemistry}\right);\mathrm{Digits}≔15;$
 $\left[{\mathrm{AOLabels}}{,}{\mathrm{ActiveSpaceCI}}{,}{\mathrm{ActiveSpaceSCF}}{,}{\mathrm{AtomicData}}{,}{\mathrm{BondAngles}}{,}{\mathrm{BondDistances}}{,}{\mathrm{Charges}}{,}{\mathrm{ChargesPlot}}{,}{\mathrm{CorrelationEnergy}}{,}{\mathrm{CoupledCluster}}{,}{\mathrm{DensityFunctional}}{,}{\mathrm{DensityPlot3D}}{,}{\mathrm{Dipole}}{,}{\mathrm{DipolePlot}}{,}{\mathrm{Energy}}{,}{\mathrm{ExcitationEnergies}}{,}{\mathrm{ExcitationSpectra}}{,}{\mathrm{ExcitationSpectraPlot}}{,}{\mathrm{ExcitedStateEnergies}}{,}{\mathrm{ExcitedStateSpins}}{,}{\mathrm{FullCI}}{,}{\mathrm{GeometryOptimization}}{,}{\mathrm{HartreeFock}}{,}{\mathrm{Interactive}}{,}{\mathrm{Isotopes}}{,}{\mathrm{MOCoefficients}}{,}{\mathrm{MODiagram}}{,}{\mathrm{MOEnergies}}{,}{\mathrm{MOIntegrals}}{,}{\mathrm{MOOccupations}}{,}{\mathrm{MOOccupationsPlot}}{,}{\mathrm{MOSymmetries}}{,}{\mathrm{MP2}}{,}{\mathrm{MolecularData}}{,}{\mathrm{MolecularGeometry}}{,}{\mathrm{NuclearEnergy}}{,}{\mathrm{NuclearGradient}}{,}{\mathrm{OscillatorStrengths}}{,}{\mathrm{Parametric2RDM}}{,}{\mathrm{PlotMolecule}}{,}{\mathrm{Populations}}{,}{\mathrm{RDM1}}{,}{\mathrm{RDM2}}{,}{\mathrm{RTM1}}{,}{\mathrm{ReadXYZ}}{,}{\mathrm{Restore}}{,}{\mathrm{Save}}{,}{\mathrm{SaveXYZ}}{,}{\mathrm{SearchBasisSets}}{,}{\mathrm{SearchFunctionals}}{,}{\mathrm{SkeletalStructure}}{,}{\mathrm{Thermodynamics}}{,}{\mathrm{TransitionDipolePlot}}{,}{\mathrm{TransitionDipoles}}{,}{\mathrm{TransitionOrbitalPlot}}{,}{\mathrm{TransitionOrbitals}}{,}{\mathrm{Variational2RDM}}{,}{\mathrm{VibrationalModeAnimation}}{,}{\mathrm{VibrationalModes}}{,}{\mathrm{Video}}\right]$
 ${\mathrm{Digits}}{≔}{15}$ (2.1)
 > $\mathrm{ESmethod}≔\mathrm{HartreeFock};$
 ${\mathrm{ESmethod}}{≔}{\mathrm{HartreeFock}}$ (2.2)
 > $\mathrm{AObasis}≔"sto-3g";$
 ${\mathrm{AObasis}}{≔}{"sto-3g"}$ (2.3)



Define Molecule

In this exercise, you will determine the geometry of the hypochlorous acid compound.  In order to generate a starting structure, begin by estimating bond lengths and angle.

 > $\mathrm{r1}≔1.0;$
 ${\mathrm{r1}}{≔}{1.00000000}$ (3.1)
 > $\mathrm{r2}≔1.7;$
 ${\mathrm{r2}}{≔}{1.70000000}$ (3.2)
 > $\mathrm{theta}≔104.5;$
 ${\mathrm{\theta }}{≔}{104.50000000}$ (3.3)

Now, convert the bond-angle coordinates to a three-dimensional Cartesian coordinate.  For simplicity, allow the molecule to occupy the 2D xy plane and the H-O bond along the x-axis: H (0,0,0) and O (r1,0,0).  Then, the coordinates of the Cl atom becomes (r1+ r2 cos (180 - theta), r2 sin(180 - theta), 0).  (Verify for yourself this is correct!)

 > $\mathrm{theta}≔\frac{\mathrm{theta}\cdot \mathrm{Pi}}{180.0};$
 ${\mathrm{\theta }}{≔}{1.82386907}$ (3.4)
 >
 ${\mathrm{molec}}{≔}\left[\left[{"H"}{,}{0}{,}{0}{,}{0}\right]{,}\left[{"O"}{,}{1.00000000}{,}{0}{,}{0}\right]{,}\left[{"Cl"}{,}{1.42564601}{,}{1.64585099}{,}{0}\right]\right]$ (3.5)
 > $\mathrm{PlotMolecule}\left(\mathrm{molec}\right);$

Geometry Optimization

To get an idea of how good (or bad) your initial guess was, first calculate the energy of the initial geometry.

 > $\mathrm{data0}≔\mathrm{ESmethod}\left(\mathrm{molec},\mathrm{basis}=\mathrm{AObasis}\right);$
 ${\mathrm{table}}{}\left({\mathrm{%id}}{=}{18446744322342913598}\right)$ (4.1)

Now optimize the geometry using the GeometryOptimization function.  This will not only determine the optimum geometry for the electronic structure method  and basis of interest but will also output the new energy:

 > $\mathrm{molec},\mathrm{data}≔\mathrm{GeometryOptimization}\left(\mathrm{molec},\mathrm{ESmethod},\mathrm{basis}=\mathrm{AObasis}\right);$
 ${\mathrm{molec}}{,}{\mathrm{data}}{≔}\left[\left[{"H"}{,}{0}{,}{0}{,}{0}\right]{,}\left[{"O"}{,}{1.00034967}{,}{-0.04822483}{,}{-1.81396399}{}{{10}}^{{-10}}\right]{,}\left[{"Cl"}{,}{1.39489050}{,}{1.64397799}{,}{1.04972729}{}{{10}}^{{-8}}\right]\right]{,}{table}{}\left(\left[{\mathrm{rdm1}}{=}\begin{array}{c}\left[\begin{array}{ccccccc}{0.55984726}& {-0.02000254}& {-0.07255045}& {-0.72400953}& {0.08298590}& {4.33352152}{}{{10}}^{{-10}}& {\dots }\\ {-0.02000254}& {2.11270990}& {-0.47981347}& {-0.06815357}& {0.06747518}& {4.15609446}{}{{10}}^{{-10}}& {\dots }\\ {-0.07255045}& {-0.47981347}& {2.09780294}& {0.38938444}& {-0.34442754}& {-2.11643907}{}{{10}}^{{-9}}& {\dots }\\ {-0.72400953}& {-0.06815357}& {0.38938444}& {0.98274723}& {-0.19112565}& {-1.32575252}{}{{10}}^{{-9}}& {\dots }\\ {0.08298590}& {0.06747518}& {-0.34442754}& {-0.19112565}& {0.80725903}& {-7.59129612}{}{{10}}^{{-9}}& {\dots }\\ {4.33352152}{}{{10}}^{{-10}}& {4.15609446}{}{{10}}^{{-10}}& {-2.11643907}{}{{10}}^{{-9}}& {-1.32575252}{}{{10}}^{{-9}}& {-7.59129612}{}{{10}}^{{-9}}& {2.01194752}& {\dots }\\ {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {}\end{array}\right]\\ \hfill {\text{15 Ã— 15 Matrix}}\end{array}{,}{\mathrm{converged}}{=}{1}{,}{\mathrm{dipole}}{=}\left[\begin{array}{c}{-1.81875578}\\ {-0.84341703}\\ {-5.51951083}{}{{10}}^{{-9}}\end{array}\right]{,}{\mathrm{aolabels}}{=}\begin{array}{c}\left[\begin{array}{c}{"0 H 1s"}\\ {"1 O 1s"}\\ {"1 O 2s"}\\ {"1 O 2px"}\\ {"1 O 2py"}\\ {"1 O 2pz"}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{mo_occ}}{=}\begin{array}{c}\left[\begin{array}{c}{2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{charges}}{=}\left[\begin{array}{c}{0.22918508}\\ {-0.15562330}\\ {-0.07356178}\end{array}\right]{,}{\mathrm{mo_symmetry}}{=}\begin{array}{c}\left[\begin{array}{c}{"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{mo_coeff}}{=}\begin{array}{c}\left[\begin{array}{ccccccc}{0.00008252}& {-0.00566888}& {-0.00183497}& {0.00199296}& {0.00151695}& {0.}& {\dots }\\ {-0.00003468}& {0.99446887}& {0.00076590}& {-0.00099367}& {0.00013630}& {0.}& {\dots }\\ {0.00022856}& {0.02399435}& {-0.00492321}& {0.00739233}& {-0.00101079}& {0.}& {\dots }\\ {0.00008282}& {-0.00280176}& {-0.00162710}& {0.00195446}& {-0.00146988}& {0.}& {\dots }\\ {0.00022205}& {0.00202907}& {-0.00421377}& {0.00589088}& {0.00017504}& {1.24664798}{}{{10}}^{{-11}}& {\dots }\\ {0.}& {1.24065356}{}{{10}}^{{-11}}& {-2.66694167}{}{{10}}^{{-11}}& {3.72449266}{}{{10}}^{{-11}}& {0.}& {-0.00198439}& {\dots }\\ {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {}\end{array}\right]\\ \hfill {\text{15 Ã— 15 Matrix}}\end{array}{,}{\mathrm{mo_energy}}{=}\begin{array}{c}\left[\begin{array}{c}{-103.72539368}\\ {-20.39140764}\\ {-10.40721342}\\ {-7.85176462}\\ {-7.84568972}\\ {-7.84548880}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{populations}}{=}\begin{array}{c}\left[\begin{array}{c}{0.77081492}\\ {1.99849901}\\ {1.90005347}\\ {1.25545769}\\ {1.00161313}\\ {2.00000000}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{group}}{=}{"C1"}{,}{\mathrm{e_tot}}{=}{-528.93124657}\right]\right)$ (4.2)

How well did your initial geometry match the optimized geometry?

 > $\mathrm{PlotMolecule}\left(\mathrm{molec}\right);$
 > $\mathrm{r1}≔\mathrm{sqrt}\left(\mathrm{molec}{\left[2,2\right]}^{2}+\mathrm{molec}{\left[2,3\right]}^{2}+\mathrm{molec}{\left[2,4\right]}^{2}\right);$
 ${\mathrm{r1}}{≔}{1.00151141}$ (4.3)
 > $\mathrm{r2}≔\mathrm{sqrt}\left({\left(\mathrm{molec}\left[2,2\right]-\mathrm{molec}\left[3,2\right]\right)}^{2}+{\left(\mathrm{molec}\left[2,3\right]-\mathrm{molec}\left[3,3\right]\right)}^{2}+{\left(\mathrm{molec}\left[2,4\right]-\mathrm{molec}\left[3,4\right]\right)}^{2}\right);$
 ${\mathrm{r2}}{≔}{1.73758823}$ (4.4)
 > $\mathrm{r3}≔\mathrm{sqrt}\left({\left(\mathrm{molec}\left[1,2\right]-\mathrm{molec}\left[3,2\right]\right)}^{2}+{\left(\mathrm{molec}\left[1,3\right]-\mathrm{molec}\left[3,3\right]\right)}^{2}+{\left(\mathrm{molec}\left[1,4\right]-\mathrm{molec}\left[3,4\right]\right)}^{2}\right);$
 ${\mathrm{r3}}{≔}{2.15601093}$ (4.5)
 > $\mathrm{theta}≔\frac{\mathrm{arccos}\left(\frac{{\mathrm{r1}}^{2}+{\mathrm{r2}}^{2}-{\mathrm{r3}}^{2}}{2\cdot \mathrm{r1}\cdot \mathrm{r2}}\right)\cdot 180}{\mathrm{Pi}};$
 ${\mathrm{\theta }}{≔}{100.36420347}$ (4.6)

Is the compound more or less bent than water?

 Answer For Hartree-Fock and STO-3G, the bond angle is less than that of water, so HOCl is more bent!



Vibrational Frequencies

Now that you have determined the optimum geometry, calculate the vibrational frequencies of the triatomic using a normal mode analysis. You may animate each normal mode using the VibrationalModeAnimation function.

 > $\mathrm{emodes},\mathrm{vmodes}≔\mathrm{VibrationalModes}\left(\mathrm{molec},\mathrm{ESmethod},\mathrm{basis}=\mathrm{AObasis}\right);$
 ${\mathrm{emodes}}{,}{\mathrm{vmodes}}{≔}\left[\begin{array}{c}{4165.35322686}\\ {1534.92740704}\\ {983.37877259}\\ {13.36540956}\end{array}\right]{,}\begin{array}{c}\left[\begin{array}{cccc}{-0.96735856}& {0.08846150}& {-0.06034427}& {0.00189287}\\ {0.05577983}& {0.95895925}& {-0.12095658}& {-0.00039001}\\ {1.03071482}{}{{10}}^{{-6}}& {-4.60627261}{}{{10}}^{{-7}}& {4.05301840}{}{{10}}^{{-8}}& {0.98654617}\\ {0.24703509}& {0.12670831}& {-0.19414483}& {0.00749708}\\ {-0.00164650}& {-0.21496833}& {-0.78847543}& {-0.00239470}\\ {-3.35259246}{}{{10}}^{{-7}}& {-5.40508443}{}{{10}}^{{-8}}& {-6.49938152}{}{{10}}^{{-7}}& {0.15982978}\\ {⋮}& {⋮}& {⋮}& {⋮}\end{array}\right]\\ \hfill {\text{9 Ã— 4 Matrix}}\end{array}$ (5.1)
 > $\mathrm{VibrationalModeAnimation}\left(\mathrm{molec},\mathrm{Re}\left(\mathrm{emodes}\right),\mathrm{vmodes},1\right);$

How well do the calculated values compare to experimental values: 3609 cm-1 (OH stretch), 1239 cm-1 (bend), and 724 cm-1 (Cl-O stretch)?

 Answer For the Hartree-Fock and STO-3G basis, we see that the motions are consistent with the assignments but the calculated frequencies are too high.

In the overview, it was mentioned that there should be only 3 normal modes for a nonlinear triatomic. If using Hartree-Fock and the STO-3G basis, note you get 4 modes, not 3.  Can you reconcile the 4th mode?

 Answer The 4th mode is a rotation, not a vibration.  This is obvious from the magnitude of the frequency, which is on the order of 20 cm-1.
 References 1. https://webbook.nist.gov/cgi/cbook.cgi?ID=7790-92-3