I want to visualize the vibrational modes of some molecules and present them alongside their predicted spectra. Using Psi4 (1.3.2) and these commands: set normal_modes_write on and molden(wfn, 'struc.molden', density_a=wfn.Da()), I was able to combine the structure and vibrational frequencies into a single Molden input file. But together, they still didn’t give me IR or Raman intensities when read by either Avogadro or Molden. A minimum working example of my h2o.in is here.
Somehow, I don’t know which command did it for me (maybe someone can let me know, so I can learn), I got the IR intensities to be summarized alongside each vibrational mode in the Psi4 output file with the command: e, wfn = freq('hf', return_wfn=True). Example output files after running the above input file are here.
I have cobbled together a Molden input file using my own python script combining the two Molden output files and the overall Psi4 output file, all produced when I run Psi4, into a combined Molden input file that has the IR intensities. But I still am not able to get something that looks like the Raman intensities, or polarizability of each mode.
As an aside, I got Avogadro to read the file produced by my script and use the vibrational information sometimes (but not every time…I don’t know why). Screen capture here. It doesn’t read in the intensities, although I feed them, but only the frequencies. Otherwise, I cannot get Avogadro to read vibrational information at all. Maybe someone could tell me why.
It is highly likely I am going about this all wrong, and I welcome someone helping point me in a better direction. Is there a command I am not finding which gives me these intensities and puts them in the right format for Avogadro or Molden to read?
Thanks for any help, and especially for your patience if this post is wordy and unhelpful.
Before I answer, a crash course in theory. When we compute infrared intensities, we usually use the electronic harmonicity approximation. That tells us that the intensity of an infrared transition is proportional to the first derivative of the electronic dipole along your normal mode. The electronic dipole itself is the first derivative of energy as you turn on an electric field. So if you’ve been counting derivatives, we’re at one derivative with respect to geometry and one derivative with respect to electric field. Now, when we compute Raman intensities, the approximation of choice is that the Raman intensity is proportional to the first derivative of the polarizability along your normal mode. Polarizability is a second derivative of energy as you turn on an electric field. So we’re at one derivative with respect to geometry and two derivatives with respect to electric field.
What does this mean for our ability to compute these things?
For a frequency by finite difference of energies, you have two derivatives with respect to geometry. There isn’t enough information to determine any intensities.
For a frequency by finite difference of gradients, you have one derivative with respect to geometry and one derivative of your wavefunction. A derivative of your wavefunction allows you to get any single-derivative you want. (Assuming you have the right integrals.) You can get an IR intensity.
For a frequency that is perfectly analytic, you have two derivatives of your wavefunction. Again, we can “spend” this to get any two derivatives of our choice, so we can get an IR intensity.
In all these cases, we can’t get a Raman intensity. That fundamentally requires three derivatives, but our frequency calculation only gives two. We’d need to do a special calculation to get Raman intensities, or take a first derivative of polarizabilities, which means doing a finite difference of hessians. I think that Crawford’s group at Virginia Tech is very interested in Raman intensities, but I’m not sure. Either way, you’re not getting those from Psi, out-of-the-box. If you’re interested, I can try to get this topic to the attention of Crawford’s students?
As for infrared intensities, those are automatically handled when you get your frequencies analytically. We have all the information we need to get the IR intensities from frequencies by gradients, too, but those aren’t coded in current Psi. They will be in 1.4, though!
I think that answers everything. Let me know if you still have questions.
