Harmonic Analysis in Psi4

Molecules vibrate, that we all know. And those vibrational signatures might be of use down the road for future research. In industry, vibration is often measured with the method of infrared spectroscopy. So let’s take our methanol from last time where we only estimated the coordinates in Z-Matrix form:

H11
O11  H11  0.9600
C11  O11  1.4000  H11  108.0000
H12  C11  1.1000  O11  112.0000  H11  -60.0000
H13  C11  1.1000  O11  112.0000  H11   60.0000
H14  C11  1.1000  O11  111.8699  H11 -180.0000

So our plan is to optimize the geometry of the molecule, once we reach a minimal state we want to return that wave function and perform a harmonic analysis on the molecule to give us vibrational modes. Whew sounds like a lot so let’s break it down, first lets get some initial configurations out the way which you can find in my previous blog here.

import psi4

psi4.set_memory('1000mb')
psi4.core.set_num_threads(1)

psi4.set_options({
     'scf_type': 'df',
     'g_convergence': 'gau_tight',
     'freeze_core': 'true',
})

zmatrix = '''\
H11
O11  H11  0.9600
C11  O11  1.4000  H11  108.0000
H12  C11  1.1000  O11  112.0000  H11    0.0000
H13  C11  1.1000  O11  112.0000  H11  120.0000
H14  C11  1.1000  O11  111.8699  H11 -120.0000
0 1
'''

universe = psi4.geometry(zmatrix)
universe.update_geometry()
universe.print_out()

Now we want to run the OPTKING from psi4 and return the wavefunction on the molecule and we and again we are going to optimize at a hartree-fock level. The wave function can be returned with the return_wfn = True

energy, wave_function = psi4.optimize(
        'hf/6-31g*',
        return_wfn =True,
        molecule=universe
)

The frequencies can then be estimated with

frequency = psi4.frequencies(
       'hf/6-31g*',
       ref_gradient=wave_function.gradient(),
       molecule=universe
)

Where we pass in the analytical gradient as the reference gradient for the frequency calculation. The molecule should have the updated geometry. Harmonic analysis will take some time and in my opinion it’s best kept to HF for initial frequency screening and switch into mp2 later on.

An example output (not the methanol, this is benzene) I just happen to have the output can be something like this:

You can see here in the output at least for benzene some of the signature frequencies that might resemble Brian’s Smith IR of benzene in his book. However, you can see some frequencies with an i attached to it meaning imaginary. It might perhaps mean that the molecule’s geometry is not optimized. Would need some further investigating perhaps.

Harmonic analysis is something I might revisit in a couple of years but I’ll leave the method here for anyone that might want to try as well.

Happy Cheminformatics!