Contents

02 Harmonic Vibrational Analysis

The purpose of this project is to extend your fundamental Python language programming techniques through a normal coordinate/harmonic vibrational frequency calculation. The theoretical background and a concise set of instructions for this project may be found here.

Original authors (Crawford, et al.) thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project.

# Following os.chdir code is only for thebe (live code), since only in thebe default directory is /home/jovyan
import os
if os.getcwd().split("/")[-1] != "Project_02":
    os.chdir("source/Project_02")
from solution_02 import Molecule as SolMol
# Following code is called for numpy array pretty printing
import numpy as np
np.set_printoptions(precision=7, linewidth=120, suppress=True)

Step 1: Read the Coordinate Data

The coordinate data are given in a format identical to that for Project #1. The test case for the remainder of this project is the water molecule (input/h2o_geom.txt), optimized at the SCF/DZP level of theory. You can find the coordinates (in bohr) in the input directory.

In this project, Molecule object can be initialized as the following toggled code.

class Molecule:

    def __init__(self):
        self.atom_charges = NotImplemented  # type: np.ndarray
        self.atom_coords = NotImplemented  # type: np.ndarray
        self.natm = NotImplemented  # type: int
        self.hess = NotImplemented  # type: np.ndarray

    def construct_from_dat_file(self, file_path: str):
        # Same to Project 01
        with open(file_path, "r") as f:
            dat = np.array([line.split() for line in f.readlines()][1:])
            self.atom_charges = np.array(dat[:, 0], dtype=float).astype(int)
            self.atom_coords = np.array(dat[:, 1:4], dtype=float)
            self.natm = self.atom_charges.shape[0]

Step 2: Read the Cartesian Hessian Data

The primary input data for the harmonic vibrational calculation is the Hessian matrix, which consists of second derivatives of the energy with respect to atomic positions.

\[ F_{r_i s_j} = \frac{\partial^2 E}{\partial q_{r_i} \partial q_{s_j}} \]

Notations:

  • \(E\): Total energy of molecule

  • \(i, j\): refer to index of atom (0, 1, 2, …)

  • \(r, s\): refer to coordinate components (\(x\), \(y\), \(z\))

  • \(q_{x_0}\): \(x\) coordinate component of atom 0

The Hessian matrix (in units of \(E_\mathrm{h} / a_0^2\), where \(E_\mathrm{h}\) stands for Hartree energy, and \(a_0\) Bohr radius) can be downloaded here (input/h2o_hessian.txt) for the H2O test case. The first integer in the file is the number of atoms (which you may compare to the corresponding value from the geometry file as a test of consistency), while the remaining values have the following format:

\[\begin{split} \begin{matrix} F_{x_1, x_1} & F_{x_1, y_1} & F_{x_1, z_1} \\ F_{x_1, x_2} & F_{x_1, y_2} & F_{x_1, z_2} \\ \vdots & \vdots & \vdots \\ F_{x_2, x_1} & F_{x_2, y_1} & F_{x_2, z_1} \\ \vdots & \vdots & \vdots \end{matrix} \end{split}\]

Hint 1: Reading and storing Hessian

The Hessian stored in memory should be a symmetric matrix, while the format of the input file is rectangular. Understanding the translation between the two takes a bit of thinking.

One may take advantage of numpy.reshape (NumPy API) to compactly convert rectangular matrix to symmetric matrix.

Implementation

Reader should fill all NotImplementedError in the following code:

def obtain_hessian(mole: Molecule, file_path: str):
    # Input: Read Hessian file from `file_path`
    # Attribute modification: Obtain raw Hessian to `mole.hess`
    raise NotImplementedError("About 2~15 lines of code")

Molecule.obtain_hessian = obtain_hessian

Solution

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o_geom.txt")
sol_mole.obtain_hessian("input/h2o_hessian.txt")
sol_mole.hess

array([[ 0.0927643,  0.       ,  0.       , -0.0463822,  0.       ,  0.       , -0.0463822,  0.       ,  0.       ],
       [ 0.       ,  0.3171327,  0.       ,  0.       , -0.1585664,  0.0800202,  0.       , -0.1585664, -0.0800202],
       [ 0.       ,  0.       ,  0.2800907,  0.       ,  0.0347766, -0.1400454,  0.       , -0.0347766, -0.1400454],
       [-0.0463822,  0.       ,  0.       ,  0.0514668,  0.       ,  0.       , -0.0050847,  0.       ,  0.       ],
       [ 0.       , -0.1585664,  0.0347766,  0.       ,  0.1730076, -0.0573984,  0.       , -0.0144412,  0.0226218],
       [ 0.       ,  0.0800202, -0.1400454,  0.       , -0.0573984,  0.1268373,  0.       , -0.0226218,  0.013208 ],
       [-0.0463822,  0.       ,  0.       , -0.0050847,  0.       ,  0.       ,  0.0514668,  0.       ,  0.       ],
       [ 0.       , -0.1585664, -0.0347766,  0.       , -0.0144412, -0.0226218,  0.       ,  0.1730076,  0.0573984],
       [ 0.       , -0.0800202, -0.1400454,  0.       ,  0.0226218,  0.013208 ,  0.       ,  0.0573984,  0.1268373]])

Step 3: Mass-Weight the Hessian Matrix

Divide each element of the Hessian matrix by the product of square-roots of the masses of the atoms associated with the given coordinates:

\[ F_{r_i s_j}^\mathrm{M} = \frac{F_{r_i s_i}}{\sqrt{m_i m_j}} \]

where \(m_i\) represents the mass of the atom corresponding to atom \(i\). Use atomic mass units (\(\mathsf{amu}\)) for the masses, just as for Project 01.

Implementation

Reader should fill all NotImplementedError in the following code:

def mass_weighted_hess(mole: Molecule) -> np.ndarray or list:
    # Output: Mass-weighted Hessian matrix (in unit Eh/(amu*a0^2))
    raise NotImplementedError("About 2~10 lines of code")

Molecule.mass_weighted_hess = mass_weighted_hess

Solution

sol_mole.mass_weighted_hess()

array([[ 0.005798 ,  0.       ,  0.       , -0.01155  ,  0.       ,  0.       , -0.01155  ,  0.       ,  0.       ],
       [ 0.       ,  0.0198215,  0.       ,  0.       , -0.0394859,  0.0199265,  0.       , -0.0394859, -0.0199265],
       [ 0.       ,  0.       ,  0.0175063,  0.       ,  0.00866  , -0.0348738,  0.       , -0.00866  , -0.0348738],
       [-0.01155  ,  0.       ,  0.       ,  0.0510613,  0.       ,  0.       , -0.0050446,  0.       ,  0.       ],
       [ 0.       , -0.0394859,  0.00866  ,  0.       ,  0.1716445, -0.0569462,  0.       , -0.0143274,  0.0224436],
       [ 0.       ,  0.0199265, -0.0348738,  0.       , -0.0569462,  0.1258381,  0.       , -0.0224436,  0.013104 ],
       [-0.01155  ,  0.       ,  0.       , -0.0050446,  0.       ,  0.       ,  0.0510613,  0.       ,  0.       ],
       [ 0.       , -0.0394859, -0.00866  ,  0.       , -0.0143274, -0.0224436,  0.       ,  0.1716445,  0.0569462],
       [ 0.       , -0.0199265, -0.0348738,  0.       ,  0.0224436,  0.013104 ,  0.       ,  0.0569462,  0.1258381]])

Step 4: Diagonalize the Mass-Weighted Hessian Matrix

Compute the eigenvalues of the mass-weighted Hessian:

\[ \mathbf{F}^\mathrm{M} \mathbf{L} = \mathbf{L} \mathbf{\Lambda} \]

where \(\mathbf{\Lambda}\) is a diagonal matrix containing eigenvalues, and \(\mathbf{L}\) contains eigenvectors. You should consider using the same canned diagonalization function you used in Project 01.

Implementation

Reader should fill all NotImplementedError in the following code:

def eig_mass_weight_hess(mole: Molecule) -> np.ndarray or list:
    # Output: Eigenvalue of mass-weighted Hessian matrix (in unit Eh/(amu*a0^2))
    raise NotImplementedError("Exactly 1 line of code using numpy")

Molecule.eig_mass_weight_hess = eig_mass_weight_hess

Solution

sol_mole.eig_mass_weight_hess()

array([-0.       ,  0.       ,  0.       ,  0.0518147,  0.0547485,  0.0561059,  0.1317344,  0.2106858,  0.2351242])

Step 5: Compute the Harmonic Vibrational Frequencies

The vibrational frequencies are proportional to the squareroot of the eigenvalues of the mass-weighted Hessian:

\[ \tilde \omega_i = \mathrm{constant} \times \sqrt{\lambda_i} \]

The most common units to use for vibrational frequencies is \(\mathsf{cm}^{-1}\) (spectroscopy wavenumber).

Hint 1: Unit dimension analysis

It should be relatively appearant that unit of eigenvalues \(\sqrt{\lambda_i}\) is the same to \(F_{r_i s_j}^\mathrm{M}\), i.e. \(E_\mathrm{h} / (\mathsf{amu} \cdot a_0^2)\).

If we regard \(E_\mathrm{h}\), \(\mathrm{amu}\), \(a_0\) as unit conversion constant, then we may have

\[ \frac{E_\mathrm{h}}{\mathrm{amu} \, a_0^2} \frac{\mathsf{J}}{\mathsf{kg} \cdot \mathsf{m}^2} = \frac{E_\mathrm{h}}{\mathrm{amu} \, a_0^2} \, \mathsf{s}^{-2} \]

So unit of \(\sqrt{\lambda_i}\) is exactly \(\mathsf{s}^{-1} = \mathsf{Hz}\). In spectroscopy, wavenumber refers to a frequency which is divided by speed of light in vacuum:

\[ \tilde \omega_i = \frac{\sqrt{\lambda_i}}{2 \pi c} \]

So finally, the unit conversion should be

\[ \frac{\mathrm{centi}}{2 \pi c} \sqrt{\frac{E_\mathrm{h}}{\mathrm{amu} \, a_0^2}} \, \mathsf{cm}^{-1} \]

All the unit conversion constants in the fomular above could be simply obtained from scipy.constants (SciPy API).

Hint 2: Imaginary frequency

In some cases, you may find some \(\lambda_i\) could be smaller than zero, making \(\sqrt{\lambda_i}\) becoming an imaginary value. This “frequency” is called imaginary frequency. If these values are far from zero, they can be some indication that the molecule is far from optimized geometry or near transition state.

We do not discuss imaginary frequency in detail here. Usually, these imaginary values is expressed in minus values in common quantum chemistry softwares (so the reader should bear in mind that minus frequencies are actually imaginary). So, the code is something like

return np.sign(eigs) * np.sqrt(np.abs(eigs)) * unit_conversion

Hint 3: “Rotational” Frequencies

For a molecule which is fully geometry optimized, there should be three lowest frequencies exactly to be zero. However, you may find there exist three frequencies near to zero. So the structure used in the computation is not exactly the stationary point on the potential energy surface (PES), although really close to that.

Implementation

Reader should fill all NotImplementedError in the following code:

def harmonic_vib_freq(mole: Molecule) -> np.ndarray or list:
    # Output: Harmonic vibrational frequencies (in unit cm^-1)
    raise NotImplementedError("About 2~15 lines of code")

Molecule.harmonic_vib_freq = harmonic_vib_freq

Solution

sol_mole.harmonic_vib_freq()

array([  -0.0000358,    0.       ,    0.020977 , 1170.1214147, 1202.791828 , 1217.611377 , 1865.752138 , 2359.5107516,
       2492.6021517])

Test Cases

Water

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o_geom.txt")
sol_mole.obtain_hessian("input/h2o_hessian.txt")
sol_mole.print_solution_02()

=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[[ 0.005798   0.         0.        -0.01155    0.         0.        -0.01155    0.         0.       ]
 [ 0.         0.0198215  0.         0.        -0.0394859  0.0199265  0.        -0.0394859 -0.0199265]
 [ 0.         0.         0.0175063  0.         0.00866   -0.0348738  0.        -0.00866   -0.0348738]
 [-0.01155    0.         0.         0.0510613  0.         0.        -0.0050446  0.         0.       ]
 [ 0.        -0.0394859  0.00866    0.         0.1716445 -0.0569462  0.        -0.0143274  0.0224436]
 [ 0.         0.0199265 -0.0348738  0.        -0.0569462  0.1258381  0.        -0.0224436  0.013104 ]
 [-0.01155    0.         0.        -0.0050446  0.         0.         0.0510613  0.         0.       ]
 [ 0.        -0.0394859 -0.00866    0.        -0.0143274 -0.0224436  0.         0.1716445  0.0569462]
 [ 0.        -0.0199265 -0.0348738  0.         0.0224436  0.013104   0.         0.0569462  0.1258381]]
=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[-0.         0.         0.         0.0518147  0.0547485  0.0561059  0.1317344  0.2106858  0.2351242]
=== Harmonic Vibrational Frequencies (in unit cm^-1) ===
[  -0.0000358    0.           0.020977  1170.1214147 1202.791828  1217.611377  1865.752138  2359.5107516 2492.6021517]

Benzene

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/benzene_geom.txt")
sol_mole.obtain_hessian("input/benzene_hessian.txt")
sol_mole.print_solution_02()

=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[[ 0.0696916  0.         0.        ... -0.0234218  0.         0.       ]
 [ 0.         0.0855755  0.        ...  0.        -0.136488   0.       ]
 [ 0.         0.         0.0156863 ...  0.         0.        -0.0153281]
 ...
 [-0.0234218  0.         0.        ...  0.0763606  0.         0.       ]
 [ 0.        -0.136488   0.        ...  0.         0.4833777  0.       ]
 [ 0.         0.        -0.0153281 ...  0.         0.         0.0365666]]
=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[-0.        -0.         0.         0.         0.         0.         0.0086435  0.0086435  0.0185369  0.0185369
  0.0248907  0.0267604  0.0407475  0.0407475  0.0503805  0.0505931  0.051977   0.0536311  0.0536311  0.0558216
  0.0568439  0.0568439  0.0711535  0.0711535  0.0717443  0.0962824  0.1188397  0.1188397  0.1412093  0.1412093
  0.5192258  0.5244019  0.5244019  0.5282622  0.5282622  0.5313338]
=== Harmonic Vibrational Frequencies (in unit cm^-1) ===
[  -0.0164509   -0.0116326    0.0246764    0.6680535    0.6680535    0.6691684  477.9142727  477.9142727  699.8795264
  699.8795269  811.0042569  840.9121042 1037.65959   1037.65959   1153.81349   1156.2446375 1171.9516278 1190.4538092
 1190.4538092 1214.5225828 1225.5931503 1225.5931506 1371.2049457 1371.2049459 1376.8857923 1595.0623276 1772.087214
 1772.0872146 1931.6832235 1931.6832236 3704.0974408 3722.514595  3722.5145955 3736.1908987 3736.1908987 3747.0369861]

3-chloro-1-butene

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/3c1b_geom.txt")
sol_mole.obtain_hessian("input/3c1b_hessian.txt")
sol_mole.print_solution_02()

=== Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[[ 0.0720844 -0.0016129  0.0003104 ... -0.0009051 -0.0017714  0.0003552]
 [-0.0016129  0.0597496  0.0045083 ... -0.0003417 -0.0006569  0.0001635]
 [ 0.0003104  0.0045083  0.0707321 ... -0.0003336 -0.0005921  0.0002486]
 ...
 [-0.0009051 -0.0003417 -0.0003336 ...  0.0073687  0.0026435  0.0017749]
 [-0.0017714 -0.0006569 -0.0005921 ...  0.0026435  0.0022558  0.0007853]
 [ 0.0003552  0.0001635  0.0002486 ...  0.0017749  0.0007853  0.0013844]]
=== Eigenvalue of Mass-Weighted Hessian Matrix (in unit Eh/(amu*a0^2)) ===
[-0.000414   0.         0.         0.         0.         0.0000001  0.0000001  0.0022826  0.0027992  0.0046311
  0.0065116  0.0088534  0.0213522  0.0283228  0.0409907  0.0498683  0.0511909  0.0569821  0.0595559  0.0640405
  0.0690233  0.0835733  0.089757   0.0940067  0.110188   0.1124069  0.1240443  0.1246022  0.1593426  0.4844511
  0.4991677  0.5062146  0.5254742  0.5350153  0.537149   0.5503098]
=== Harmonic Vibrational Frequencies (in unit cm^-1) ===
[-104.5906768    0.0099111    0.0289641    0.0304254    0.3665457    1.5129799    1.7406089  245.5941566  271.9691098
  349.8231157  414.8104445  483.6819181  751.1480607  865.1127581 1040.7507835 1147.9332085 1163.0562807 1227.0817046
 1254.488878  1300.863303  1350.5239405 1486.0657982 1540.0631407 1576.0995807 1706.3629779 1723.4581894 1810.4757465
 1814.5421487 2051.9665893 3577.9088302 3631.8466432 3657.3927489 3726.3186775 3759.995898  3767.4862956 3813.3608088]

References

  • Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations Dover Publication Inc., 1980.

    ISBN-13: 978-0486639413