Infrared intensities of furan, pyrrole and thiophene: beyond the double harmonic approximation

Infrared intensities of the fundamental, overtone and combination transitions in furan, pyrrole and thiophene have been calculated using the variational normal coordinate code MULTIMODE.

We use pure vibrational wavefunctions, and quartic force fields and cubic dipole moment vector surfaces, generated by density functional theory.

The results are compared graphically with second-order perturbation calculations and with relative intensities from experiment for furan and pyrrole.


This paper is a continuation of our previous studies of vibrational transitions in five-membered aromatic ring compounds.1

We have found that both perturbative and variational approaches, based on the fourth degree Taylor expansion of the potential, yield energies of vibrational levels for fundamental and overtone transitions in excellent agreement with experimental data.

The next logical step is to expand the dipole moment vector as a Taylor series in normal coordinates and calculate anharmonic contributions to the infrared intensities.

This approach was tested recently for water and formaldehyde.2

While for small molecules there are many avenues along which to pursue infrared intensities of higher quality, for larger molecules (with more than five atoms) expansion in normal coordinates appears to be the only practical approach.

A number of scientists are now recognising the importance of comparing theoretical spectra with observation.

For example Hirano et al3. have recently used the Morbid program4 to predict and compare spectra for MgNC and MgCN.

The paper is organised as follows.

In Section 2 we discuss the potential and dipole fields, variational and perturbation approaches used to evaluate the transition dipole moments and calculation of infrared intensities from these moments.

In Section 3 we present and discuss the results for furan, pyrrole and thiophene.

In Section 4 we summarise the findings and outline future directions.


To study molecular vibrations it is necessary to solve the Schrödinger equation for nuclear motion.

While there are many different coordinate systems in which the kinetic and potential energy operators can be expressed, for larger molecules the only expansion of the kinetic energy operator, which is not prohibitively complicated, is in terms of normal coordinates.5

Therefore it is advantageous to express the potential energy also in terms of the normal coordinates.

If analytical second derivatives φjk are calculated at the equilibrium geometry and at the geometries displaced along the normal coordinates, cubic and semi-diagonal quartic force constants, φijk and φiijk, can be obtained by numerical differentiation: If infrared intensities are desired, also dipole moment field needs to be evaluated.

Since the dipole moment is a derivative of the energy with respect to the electric field, the equivalent of the quartic expansion of the potential is a cubic expansion of the dipole surfaces.

The second and third derivatives of dipole moment, dαij and dαiij (where α is a Cartesian component of the dipole moment vector dαj referred to Eckart axes), are calculated in a manner analogous to the cubic and quartic force constants: Potential and dipole moment surfaces were calculated by the Density Functional Theory (DFT) approach using the B97-1 functional6 and a TZ2P basis set.

Calculations were performed using the Cadpac program package.7

The step size δQi was chosen such that it would correspond to an energy step of 1 mEh on the harmonic potential surface.

A fuller discussion of the finite difference calculation, the effects of the step size, as well as values for the anharmonic constants may be found in .ref. 1

The variational calculations were performed using the program MULTIMODE.8–10

This program performs vibrational self-consistent force field (VSCF) calculations, followed by the truncated Configuration Interaction (CI) expansion in virtual VSCF functions.

This has an advantage of explicitly treating all resonances among states included in the CI expansion.

Terms coupling up to three normal coordinates were used in the VSCF Hamiltonian.

The CI expansion included all single, double, triple and quadruple VCI excited states with the sum of vibrational quantum numbers less than or equal to four.

Finally, the matrix elements Rαij were evaluated using the dipole moment surface Dα(Q) and the CI vectors Ψ(Q): Rαij = ∫Ψi(Q)Dα(Q)Ψj(Q)dQSecond-order perturbation calculations were carried out using the program SPECTRO.11

While this program is capable of treating low-level resonances by removing the divergent terms and explicitly solving 2 × 2 CI problems (or by providing matrix elements for larger CI problems), second-order perturbation theory does not handle resonances involving states with the sum of vibrational quantum numbers being larger than 2.

While there are many resonances present in the systems studied, none of the resonances are very strong.

Therefore, we decided not to include any resonances in the present perturbation calculations with SPECTRO.

For more information about the programs MULTIMODE and SPECTRO and their use, see for example refs. 1 and 2.

The infrared intensities of a vi → vj transitions were evaluated using the usual expression: where NA is Avogadro's constant, νa is the wavenumber of the transition, and Nv is the fraction of molecules in the state v.

If we assume an equilibrium distribution, then for the transitions originating from the ground vibrational state the term (Nv – Nv) tends to 1.

Evaluating the fraction and converting to practical units leads to a working formula with νa in cm–1 and Rαij in debye.

This formula does not sum over hot and stimulated emission bands, which can all be treated explicitly.

It also allows us to treat all kinds of transitions on an equal footing, which is especially important for variational calculations which mix transitions of various characters (e.g. fundamental, overtone, combination).

Results and discussion

The stick spectra representing vibrational transitions from the ground vibrational state, obtained by variational and perturbational calculations for furan and pyrrole, are compared with scaled experimental results of Mellouki et al12. in Figs. 1 and 2, respectively.

Since, to our knowledge, there is no consistent set of experimental infrared intensities for thiophene available, Fig. 3 compares only calculated results.

The peaks corresponding to fundamental transitions νi are labeled by i.

Transitions for which experimental intensities were not available are denoted with a star.

It can readily be seen that the agreement between the stick spectra calculated by variational and second-order perturbation approaches is, in general, very good, especially for furan and thiophene.

For pyrrole there are some disagreements, most notably for the states involving ν16.

As observed previously1 this mode shows large negative anharmonicity and is probably not adequately described by a quartic Taylor expansion.

Also the density of states is higher for pyrrole, which leads to an increase in the number and strength of resonances.

It seems that inclusion of resonances would be necessary for a proper perturbation description of intensities in pyrrole; however, this would eliminate much of the simplicity of the second-order perturbative approach, and this is one of its most attractive features.

(A referee has observed that it is difficult to see from Figs. 2(a) and (b) that the variational treatment represents an improvement over the perturbational approach.

We agree, but we anticipate that as the molecules become larger, as a force field becomes more complicated with more complicated resonances, the difference between the variational and perturbational spectra will become more marked.

We underline a fact that only the variational approach can include all resonances, and therefore it must be recommended when it is possible.)

In the region above 2000 cm–1 there are more states shown in the variational calculation panels than in the corresponding second-order perturbation calculation and experiment panels.

This is due to the fact that MULTIMODE yields all states (as long as they are included in the CI expansion), while SPECTRO produces only states with the sum of vibrational quantum numbers less than or equal to two.

Higher excited states can be produced, but the quality of the results is not the same.

Also, only states which were properly assigned are included in the experimental panel.

One of the consequences of the higher number of states and explicit treatment of all couplings in MULTIMODE are apparent in the lower intensities for high-frequency C–H and N–H stretches obtained from MULTIMODE, as compared with those obtained from SPECTRO.

In fact, all of these modes are strongly coupled to other modes, resulting in several modes each with lower intensity.

However, if a stick spectrum is replaced by gaussians with finite width, the summed intensities are comparable.

Comparison with experiment is complicated by the lack of firm experimental data.

Mellouki et al12. published relative intensities for furan and pyrrole.

However, the authors caution that the data are not accurate.

For comparison with our stick spectra we have scaled the relative intensities of ref. 12 by factors of 22.2 for furan and 163 for pyrrole, respectively.

Despite the above mentioned uncertainties in the experimental spectrum, it can be seen that the correspondence between theoretical and experimental results is good for these two molecules.

It can be observed that overtone and combination bands obtained from MULTIMODE are shifted towards higher wavenumbers with respect to those obtained from SPECTRO.

This is due to the relatively small size of the CI.

To calculate the wavenumbers of overtones it is necessary to include states with the sum of vibrational quantum numbers larger than four.

A proper description of combination bands further requires the inclusion of quintuple excitations in the CI.

Nevertheless, the intensities of overtone and combination bands are in reasonable agreement with those obtained from SPECTRO and with the scaled experimental results.


In this paper we have demonstrated the success of the variational code MULTIMODE to predict integrated band intensities for molecules.

All one needs are quartic potential energy surfaces and cubic dipole surfaces, expanded in normal coordinates.

An ab initio or DFT quantum chemistry code which calculates analytic second derivatives is also required.

The quality of our results is immediately apparent by comparing Fig. 1(a) (MULTIMODE) with Fig. 1(c) (experiment), for furan, and likewise Figs. 2(a) and (c) for pyrrole.

We have every reason to expect that the MULTIMODE predicted spectra (Fig. 3(a)) for thiophene would be close to observation.

The new feature of our science is its automatic capability of working beyond the double harmonic approximation, thus including anharmonicity, overtones and combination bands in our calculations.

Of course we have only looked at the wider detail.

Our earlier studies show that frequencies determined by DFT are in error by 1%, and the calculated intensities are also only very approximate, and without fine structure.

Even so, we believe that these qualitative spectra will be useful in molecular identification, at least.

A caveat to these results is that ideally we should perform as large a CI as possible within MULTIMODE in order to ensure that the overtones and combination bands are fully converged.

The larger the molecule, the more mixing is present.

Finally, we are in the process of calculating ro-vibrational wavefunctions, from which band spectra can be obtained (for low J).