Issue with GBS for Molecular Vibronic Spectra: Formic Acid

Hello!

To test Strawberry Fields, I’ve decided to compute the vibronic spectra of formic acid using GBS based on the electronic structure obtained from Firefly (or GAMESS). However, I would get two issues:

  • Compared with the data provided in the strawberryfields.apps.data package, the number of modes are not the same;
  • The vibronic spectra are very different.

More explicitly, on my side, Firefly gives a total of 15 modes, so 9 vibronic modes, compared with 7 for data in strawberryfields.apps.data. As for the spectrum, I get the following:

As with my other post, I wish I could attach the data but I can’t since I’m a new user. Let me know how I should share the Firefly output files.

Thank you!

Andre

Hello @andre,

Thank you for your question!

Regarding your first point:

Compared with the data provided in the strawberryfields.apps.data package, the number of modes are not the same

The data for formic acid provided in strawberryfields.apps.data was taken directly from this paper by Joonsuk Huh et al. Please, see Sec. 2 of the supporting information.

You are right, formic acid has 9 vibrational modes. However, due to the symmetry of the molecule, the Duschinsky matrix for this vibronic transition is block diagonal with two blocks, one of dimension 7x7 (totally symmetric block) and the other with dimension 2x2 (non-totally symmetric block). If the overall shape of the spectrum is dominated by the totally symmetric block you can neglect the other two modes.

Second point:

  • The vibronic spectra are very different.

There may be several aspects related to this.

1- The vibronic spectrum you want to simulate corresponds to the photoelectron emission of formic acid. So, you need to compute the vibrational modes of the ground states of the neutral and positively charged molecule in order to create the input data for the GBS simulations.

2- Which type of electronic structure calculations are you performing to optimize the molecular geometries and compute the vibrational frequencies? Could you please provide us with the frequencies you have obtained? If they are different with respect to the ones reported in the above mentioned paper, the spectrum will be different.

3- In principle, you don’t need to account for the symmetry to obtain these values. I would expect that using the B3LYP functional with a Pople-type basis set e.g., 6-311++G** would be enough to reproduce the frequencies in the data set. GAMESS and Firefly should give you similar results.

4- Assuming correctness of points 1-3, how many samples are you taking to reconstruct the vibronic spectrum?

Hope this helps.

Thank you!

Hi @Alain_Delgado_Gran,

Thank you very much for the detailed answer.

You are right, formic acid has 9 vibrational modes. However, due to the symmetry of the molecule, the Duschinsky matrix for this vibronic transition is block diagonal with two blocks, one of dimension 7x7 (totally symmetric block) and the other with dimension 2x2 (non-totally symmetric block). If the overall shape of the spectrum is dominated by the totally symmetric block you can neglect the other two modes.

Makes sense. For the data I use, however, it doesn’t look like there is a 7x7 totally symmetric block (data is included below). Would a non-symmetric Duschinsky matrix cause issues in the back-end?

1- The vibronic spectrum you want to simulate corresponds to the photoelectron emission of formic acid. So, you need to compute the vibrational modes of the ground states of the neutral and positively charged molecule in order to create the input data for the GBS simulations.

That is what we do. However, we are relaxing the geometry of the excited state to avoid imaginary frequencies. Is this something you do too?

As for points 2-3, we do a Hartree-Fock-type calculation with the basis set being TZV or similar to 6-311+G**. However, we don’t expect the spectrum to be radically different; peaks may shift but the overall shape should be similar.

4- Assuming correctness of points 1-3, how many samples are you taking to reconstruct the vibronic spectrum?

In the spectrum shown in the last post, there are exactly 4675 samples. Referring to the vibronic spectra tutorial, it seems that with 10 samples you can already find the main peaks of the spectrum, pointing towards the fact that 4675 samples is more than enough to recover the overall shape.

Moreover, using a conventional method, we were able to use the same data and compute a spectrum that is approximately correct, compared to experimental data.

This may indicate that I’m simply not using Strawberry Fields correctly. Thus, I’m including the code used to run the simulation along with the data files.

I’d very much appreciate if you could help me fix this issue.

Thank you!

test.py (952 Bytes) formicAcidHF_Hess.txt (122.4 KB) formicAcidCationHF_Hess.txt (125.0 KB)

Hi @andre,

Thank you for your message and the attached files.

Just to be clear. In order to compute the vibrational frequencies of the initial and final states you have to optimize the geometries of the neutral and the cation molecules. I noted that your electronic structure calculations were performed at the level of Hartree-Fock. That’s why your frequencies are blueshifted w.r.t. the values used in the Vibronic spectra tutorial. If you want to refine these numbers you can run density functional theory (DFT) simulations using the B3LYP functional. This way you would include electronic correlation effects and refine your frequencies. I assume that in Firefly you can do this by adding DFTTYP=b3lyp to the $CONTRL group in your input file.

Having said that, I would also expect that this refinement should not change drastically the relative intensities of the vibronic spectrum.

We will take a look at the python script and get back to you soon!

Thank you.

Hi @andre. Please also note that the normal modes used to compute the Duschinsky with the duschinsky function need to be mass-weighted. I am not sure if Firefly prints out mass-weighted data but you can use:

Li = Li[6:].T * m[:, np.newaxis]**0.5
Lf = Lf[6:].T * m[:, np.newaxis]**0.5

and then call the duschinsky function. For more information please check the documentation here.

I also highly recommend using DFT or a post-HF method with a basis set with proper size to get the input vibrational data. Please let us know if you face any issues.

Hi @sjahangiri,

Thanks for the tip for mass-weighting the normal modes, if made explicitly in the documentation I’m sure it would help others. However, I still believe the resulting spectrum (shown below, with 1000 samples) is far from being in agreement with the experimental results (and with the results obtained using the conventional method as shown in my second post):

If you believe the discrepancies could be solved by simply using DFT, please let me know. Thanks again!

Hi @andre. Please see a detailed explanation of the input format, including the normal modes format, for the duschinsky function here. I recommend using DFT, for example PBE0 or B3LYP functionals, to compute the input data rather than HF. That, very likely, will give you the correct spectrum. Please let us know about the results. Thanks.

Hi @sjahangiri,

Concerning the lack of proper documentation, I meant that it would be beneficial to have the process of using GAMESS data documented (i.e., skip over the first 5-6 modes and mass-weight the normal modes).

As for the vibronic spectrum calculation of formic acid, we’ve tried with DFT (B3LYP) and although the spectrum does look better (spectrum and data attached below), I am not entirely satisfied and would like some pointers on how to obtain better results.

Also, if possible, could you please explain the procedure of removing the two independent modes to perform the computation with 7 modes (as described above by @Alain_Delgado_Gran)?

Thank you very much!

formicAcidDFTHESS.TXT (190.4 KB) formicAcidCationDFTHESS.TXT (194.0 KB)

Hi @andre.

Thanks for the suggestion on documenting the process of using GAMESS data, we will consider that in the next refactor.

The fact that the DFT spectrum has a better agreement with the experimental one, compared to the one you obtained with HF, suggests that the shape of the spectrum is very sensitive to the vibrational input data obtained from electronic structure calculations. So, my recommendation is to further improve the electronic structure calculations by using a functional that gives more accurate vibrational frequencies, using a better basis set or using post-HF ab initio methods.

Regarding your last question, if you look at the vibrational modes of formic acid, you can see two groups of vibrations: one group of in-plane vibrations that contains 7 modes and a group of out-of-plane vibrations with two modes. In this paper, the authors have observed that the total vibronic spectrum is dominated by the modes of the first group and the second one can be safely neglected because of the lack of structural displacement. The vibrational frequencies for these two out-of-plane modes are 687.27, 1056.82 cm^-1, which are obtained from a DFT calculation using PBE0 functional. For more information, you can have a look at page 8 of the reference paper.

Thank you very much @sjahangiri and @Alain_Delgado_Gran, I’ve marked the thread as solved.