Question about Basis Sets to IR Frequencies, Plotting, and Mode Selection.

[Sulstice]

Hello,

Question 1

I was wondering if you could give me the scaling factors for mp2 level of theory and basis sets: aug-cc-pvdz and cc-pvqz for harmonic analysis.

Question 2

I read a lot about IR intensities in our previous threads about it being released in psi4 1.5, is this feature released yet?

Question 3

There is a keyword called "mode" to specify the frequency of interest and I saw it in the documentation but didn't know what parameters I could specify. Is there a way to isolate only anharmonic if that makes sense.

Thank you, -Sulstice

[loriab]

1. Psi doesn't tabulate these, and I've mostly seen them for DFT. Scaling factors handy for getting oriented with a vib spectrum, but they're not uniquely determined for a method/basis.
2. Yes, IR intensities are available as of psi4 1.6 in late May provided analytic gradients are available (as opposed to Hessian by finite difference of energies).
3. I'm not sure what "mode" you were looking at. I can clarify if you post a link. Once upon a time, there was a frequency(..., mode=(continuous|sow|reap)) but that had to do with farming out finite difference jobs and is now defunct anyways. There's irrep https://github.com/psi4/psi4/blob/master/psi4/driver/driver.py#L1511-L1517 that allows a partial frequency job for only a symmetry subset. Example here https://github.com/psi4/psi4/blob/master/tests/fd-freq-gradient/input.dat#L75 . There's no singling out of particular modes.

[Sulstice]

1.) Thank you, I read it from a graduate student thesis I remember summarizing it with different levels of basis sets and theories but definitely needed some validation on that front.

2.) Awesome I will give it a shot and check it out, let you know if I run into any problems.

3.) mode=(continuous|sow|reap)) ah I think I got confused on the parameter meant. I thought it was a way to single out individual modes.

This is great and moves me a long! Thank you!

[Sulstice]

Would you mind teaching me if I have this right because I don't actually know and this is how far I got, I still can't seem to get intensities out.

import psi4
import textwrap
import numpy as np

psi4.core.set_num_threads(8)
psi4.set_memory('30000mb')

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

psi4.core.set_output_file('water.out', False)

def run_calculation(molecule):

    universe = psi4.geometry(molecule)
    universe.update_geometry()
    mass = np.asarray([16.01, 1.0, 1.0])
    geometry = np.asarray(universe.geometry())
    irrep_labels = universe.irrep_labels()
    dipole_derivatives = None
    project_translation = True
    project_rotation = True
    symbols = [universe.symbol(at) for at in range(universe.natom())]
    theory = 'mp2/aug-cc-pvdz'

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

    hessian, wave_function_2 = psi4.hessian(
       theory,
       ref_gradient=wave_function.gradient(),
       return_wfn= True
    )

    basisset = wave_function_2.basisset()

    wave_function_2.hessian().print_out()

    vibinfo, vibtext  = psi4.driver.qcdb.vib.harmonic_analysis(
        np.array(hessian),
        geometry,
        mass,
        basisset,
        irrep_labels,
        dipole_derivatives,
        project_translation,
        project_rotation
    )

    print(vibtext)
    print(psi4.driver.qcdb.vib.print_vibs(vibinfo, shortlong=True, normco='q', atom_lbl=symbols))

if __name__ == '__main__':

    water_zmatrix = '''\
        O
        H 1 0.9894093
        H 1 0.9894093 2 100.02688
    '''
    run_calculation(textwrap.dedent(water_zmatrix))

The output from the script for water is:

  Vibration                       7                   8                   9           
  Freq [cm^-1]                1563.0797           4068.4404           4208.3554       
  Irrep                                                                               
  Reduced mass [u]              1.0740              1.0365              1.0684        
  Force const [mDyne/A]         1.5461             10.1080             11.1487        
  Turning point v=0 [a0]        0.2678              0.1690              0.1636        
  RMS dev v=0 [a0 u^1/2]        0.1962              0.1216              0.1196        
  Char temp [K]               2248.9237           5853.5800           6054.8864       
  ----------------------------------------------------------------------------------
      1   O                0.00  0.00 -0.27    0.00  0.00  0.19   -0.00 -0.26  0.00   
      2   H               -0.00  0.41  0.54    0.00  0.58 -0.39   -0.00  0.52 -0.44   
      3   H               -0.00 -0.41  0.54   -0.00 -0.58 -0.39    0.00  0.52  0.44 

Everything else but that, what am I missing so far?

[loriab]

At first read-through, that looks reasonable. What are your goals, though? I ask because all the hessian ... print_vibs section can accomplished by replacing the hessian() call with a frequency() call. IR intensities are absent because dipole_derivatives=None. frequency() will pass the data around behind the scenes. Also, there's good physical reasons to do the optimization and frequency at exactly the same level of theory (mtd + basis). Hessian is more expensive than opt, so if anything, the freq is the cheaper level of theory.

[Sulstice]

I'm looking for anharmonic modes of vibration of complexes. Good point on the optimization and frequency, for testing purposes I bring the level of theory down to hartree-fock just to make it go faster. On production code runs, the level of theory and basis set are consistent.

I have been playing around with both frequency() and hessian(). I chose the hessian last night because of this:

https://psicode.org/psi4manual/master/api/psi4.driver.qcdb.vib.harmonic_analysis.html#psi4.driver.qcdb.vib.harmonic_analysis

https://psicode.org/psi4manual/master/freq.html

Where the harmonic analysis is documented showing the IR intensities. The first parameter was a hessian matrix so I went back to go look at how to produce that. Maybe I went down a different rabbit hole.

https://github.com/psi4/psi4/blob/821134f62396ba27f9bcb8fbfa93ea2c370b7616/tests/pytests/test_vibanalysis.py

Line 17-40 I kind of copied your guys test to get the code running. I was actually confused with dipole derivatives and how to produce them and pass them in appropriately.

Ah okay! the frequency has it built in and I can see it. I got confused on the docs.

theory = 'hf/6-31g*'

    energy, wave_function = psi4.optimize(
        theory,
        return_wfn = 'yes',
        molecule=universe
    )

    frequencies = psi4.frequencies(
       theory,
       ref_gradient=wave_function.gradient(),
       molecule=universe

    )

Ouput:

Freq [cm^-1]                1557.5017           4053.2831           4197.9898       
  Irrep                           A1                  A1                  B2          
  Reduced mass [u]              1.0830              1.0449              1.0829        
  Force const [mDyne/A]         1.5478             10.1143             11.2435        
  Turning point v=0 [a0]        0.2672              0.1686              0.1627        
  RMS dev v=0 [a0 u^1/2]        0.1966              0.1219              0.1198        
  IR activ [km/mol]            92.6794             13.8599             85.5429        
  Char temp [K]               2240.8982           5831.7720           6039.9727       
  ----------------------------------------------------------------------------------
      1   O                0.00 -0.00 -0.07    0.00 -0.00  0.05    0.00 -0.07 -0.00   
      2   H               -0.00  0.43  0.56    0.00  0.59 -0.39    0.00  0.56 -0.43   
      3   H                0.00 -0.43  0.56    0.00 -0.59 -0.39    0.00  0.56  0.43   

So I want to make sure I understand this correctly, the IR active means it is my epsilon in beer-lambert's law and all I would need to do to get absorbance is times it by the path length of my cell, and the concentration of my sample. And then calculate transmittance from how much was absorbed to how much light was emitted?