The aromaticity of cyanine dyes: ring current and valence bond calculations

Abstract

Cyanine dyes are interesting molecules due to their optical properties, in particular their absorption in the near infrared window. Odd/even effects for these dyes were observed [N. P. Neme, T. L. C. Jansen and R. W. A. Havenith, Phys. Chem. Chem. Phys., 2024, 26, 6235–6241], due to interaction with the polymethine chain. Here, we present a study of the aromatic properties of oxidized cyanine dyes 1(n) using the ipsocentric method to calculate induced ring currents and using valence bond theory. The even members of this series show an induced diatropic ring current, while the odd members show an indued paratropic ring current. Hence, for n even, the ring is aromatic. The aromatic behaviour of the five-membered ring aligns with the possibility to draw ionic resonance structures in the five-membered ring that allows the delocalization of a negative charge within the ring, leading to a six π electron ring. For n odd, the induced ring current is paratropic, while a positive charge delocalizes in the five membered ring (a 4π electron circuit), hence, its aromatic behaviour. This steering of aromatic/antiaromatic properties can have a profound effect on the reactivity of the molecules, and impact their photophysical properties. This may be an extra handle for the design of organic molecules for electronic devices.

---

Introduction

Cyanine dyes belong to the class of fluorescent dyes, that absorb in the near infrared region. Such dyes have been used in various applications, such as molecular switches for surface and electronic applications,¹ sensitizers in organic solar cells,² and they are applied in the field of biological imaging and drug delivery.³ The heptamethine cyanine dyes are frequently used for these applications, as they possess high molar absorption coefficients, and have a tunable structure. This particular dye consists of a polymethine chain and two heterocyclic end-groups, and thus have an extended, linear π system. As a consequence of this extended π system, these dyes have unique photophysical properties; it typically has an intense and narrow absorption band in the near infrared I window (700–1000 nm), combined with a blue shifted shoulder, typically assigned to the 0 → 1 vibronic transition.⁴ A rigid, cyclohexene or cyclopentene, bridge can be introduced into the polymethine chain to prevent photoinduced cis–trans isomerizations, thereby increasing the thermal and photostability of the dye.

Here, we denote the cyanine dyes by the length of the cyanine, and in our notation 1(n), the number of carbon atoms (N_C) in the main chain is equal to N_C = 2n + 1, with n = 1, 2,…. In previous work, the cyanine dye 2(2) (Scheme 1) has been synthesized.^5–7 As a by-product, the oxidized dye 1(2) was found. Interestingly, this dye has, despite having more unsaturated bonds, a blue-shifted main absorption peak compared to its parent. A computational study to explain this unexpected behaviour was conducted and it was shown that this was caused by the interaction of the extra π-bond with the linear conjugated π-system of the polymethine chain.⁸ For pristine cyanine dyes in C_2v symmetry, the symmetries of the HOMO and LUMO alternate, and for even n the HOMO has a₂ symmetry and the LUMO b₁, while for odd n the HOMO has b₁ and the LUMO a₂ symmetry. As the symmetry of the additional bonding π-orbital in the bridge is b₁ and of the anti-bonding orbital a₂, an odd/even effect is expected as the interactions of the bonding/anti-bonding π-orbitals will alternate with the alternating HOMO/LUMO symmetries of the polymethine system. This odd/even affect with changing length of the polymethine chain prompted us to study the properties, and in particular the aromaticity, of the five-membered ring as a function of the length of the polymethine chain.

Scheme 1 The structure of the cyanine molecules 1(n), and 2(2).

Aromaticity, though not always well-defined,⁹ is, according to the magnetic criterion for aromaticity, the ability to sustain a diatropic ring current, induced by an external magnetic field.^9–14 Antiaromaticity is defined as the opposite, the ability to sustain a paratropic induced ring current. The induced ring current can be calculated using the ipsocentric CTOCD-DZ^15,16 approach.^17,18 This approach has many advantages, such as well-converged maps with modest basis sets,¹⁹ interpretable orbital contributions,²⁰ and an even further breakdown of these contributions in terms of virtual transitions.²¹ The diatropic contributions to ring current are caused by virtual translational excitations from occupied to unoccupied orbitals, while paratropic contributions are determined by virtual rotational transitions.^17,18

A complementary view to aromaticity is given in terms of resonating Lewis structures. The wavefunction of a system can be written in terms of resonating Lewis structures using valence bond theory.^22,23 In this way, the contribution of each resonance structure can be calculated, together with a resonance energy, using different orbital models.^24,25 In this work, we will study the aromaticity of the bridging five-membered ring in the cyanine dyes 1(n), with n = 2–11. We present plots of the induced current density, together with a detailed orbital analysis, and we interpret the electronic structure of the molecules in terms of contributing resonance structures. We show that an odd/even effect exists, and that the five-membered ring is alternating aromatic/antiaromatic, shown by a diatropic/paratropic ring current. This behaviour can be explained in terms of the orbital contributions and virtual transitions governing the ring current. Furthermore, the valence bond analysis shows an odd/even effect, due to contributions of ionic structures and that the number of π electrons in these ionic structures in the five-membered ring alternates between four and six.

Computational information

The geometries for the molecules 1(n), with n = 2–11 were optimized within the C_2v point group, using the PBE0 functional^26–29 combined with the 6-311G** basis set.³⁰ All stationary points were characterized as genuine minima by hessian calculations (no imaginary frequencies). For each molecule, UPBE0 and UHF calculations were performed in order to find broken spin symmetry solutions with a lower energy than the closed shell solution. All calculations were performed with the GAMESS-UK³¹ package, except the TDDFT^32,33 calculations, which were performed using Gaussian16³⁴ (the PBE0 and 6-311G** basis set were used here as well).

The ring current calculations were performed using the CTOCD-DZ method^{16–18,35,36} as implemented in GAMESS-UK^31,37 and SYSMO³⁸ using the PBE0 functional and 6-311G** basis set. The ring current is plotted in a plane 1 a₀ above one of the aromatic rings. For the ring current calculation of the broken spin symmetry species 1-BS(9), UHF was used (see text), following established protocols published in the literature.^39,40 In all the plots, diatropic (paratropic) current is anticlockwise (clockwise).

The valence bond calculations on 1(2) and 1(3) were performed with TURTLE,^41,42 as implemented in GAMESS-UK, using the 6-311G** basis set. The σ orbitals were obtained from a preceding RHF calculation, and were kept frozen. The p-orbitals to describe the π system were kept strictly atomic, and were optimized using the VBSCF procedure.^43,44 The structures depicted in Scheme 2 were included in the wavefunction, and their weights were determined using the Gallup–Norbeck (GN) scheme.⁴⁵

Scheme 2 The considered valence bond structures for 1(2) and 1(3). The number of symmetry equivalent structures is indicated.

Results and discussion

Geometries

The optimized C–C bond lengths for 1(2) and 1(3) are depicted in Scheme 3, and the C¹–C² and C²–C³ bond lengths for the series 1(n), n = 2–11 are summarized in Table 1. The main differences between the molecules are found for the C¹–C² and C²–C³ bond lengths, the other C–C and C–N bond lengths are very similar in the series. A clear difference is discernible between the molecules with even n and those with odd n. The C¹–C² bond length for the molecules with odd n is consistently shorter than that for the even n molecules in the series, while the C²–C³ bond length is consistently longer. The larger bond length alternation in the n odd series hints towards a less degree of aromaticity of the five-membered ring in the odd series than in the even series. The difference between the odd and even series in bond length alternation diminishes when the molecule is elongated, suggesting a decrease in the difference in aromaticity in the five-membered ring.

Scheme 3 PBE0/6-311G** optimized geometries for 1(2) and 1(3). Indicated bond lengths are in Å.

Table 1 The PBE0/6-311G** optimized C¹–C² and C²–C³ bond lengths (in Å) for 1(n), n = 2–11

Molecule	C^1–C^2	C^2–C^3	Molecule	C^1–C^2	C^2–C^3
1(2)	1.355	1.453	1(3)	1.337	1.495
1(4)	1.352	1.459	1(5)	1.339	1.487
1(6)	1.351	1.461	1(7)	1.340	1.484
1(8)	1.351	1.462	1(9)	1.341	1.483
1(10)	1.351	1.461	1(11)	1.342	1.480

---

When checking for triplet instabilities in the wavefunction, by trying to converge on a spin symmetry broken unrestricted solution with lower energy than the closed shell solution, it was found that using the PBE0 functional all closed shell solutions were stable. At this level of theory, no broken symmetry solutions were found for the molecules in the series. However, when unrestricted Hartree–Fock (UHF) calculations were performed, broken symmetry solutions were found for 1(9) and longer. The UHF solutions were found to be lower in energy than the closed shell RHF solutions for those molecules (Table 2). The energy lowering with respect to the closed shell RHF solution obtained using the broken symmetry approach becomes larger when the molecule becomes longer. The degree of spin contamination (Table 2) also increases upon chain elongation.

Table 2 ΔE (= E_BS − E_CS, the energy difference between the broken symmetry solution and the closed shell RHF energy, in eV) and the expectation value of the S² operator for 1(n), n = 9–11

Molecule	ΔE	〈S^2〉
1(9)	−0.49	2.48
1(10)	−0.62	2.59
1(11)	−1.15	3.08

---

Induced current density

The aromaticity of the five-membered ring was studied using the ring current criterion. Plots of the induced π current density for 1(2), 1(3), and 2(2) are shown in Fig. 1. A clear diatropic current is discernible for 1(2) in the five-membered ring, indicating that this ring has aromatic character. This ring current is dominated by the contributions from the HOMO and HOMO−1 (Fig. 2), and can thus be seen as a classical 4π electron current, as has been observed in prototypical aromatic molecules.^17,18 The contributions of the other π orbitals to the current in the five-membered ring are negligible.

Fig. 1 Plots of the induced π current density for 1(2), 1(3), and 2(2).

Fig. 2 Orbital contributions to the π current density for 1(2) and 1(3). Upper panels, the sum of the contributions of the HOMO and HOMO−1, lower panels the sum of the contributions of the remaining π orbitals.

The ring current pattern for 1(3) is distinctly different: here, a paratropic π ring current is found (Fig. 1). This ring current is also dominated by the (paratropic) contributions from the HOMO and HOMO−1 (Fig. 2). Again, the remaining π orbitals do not contribute to the π current density significantly. Hence, in 1(3), the five-membered ring shows the signature of anti-aromatic character.

As a reference, the ring current has also been calculated for 2(2) (Fig. 1); this molecule has a saturated C₂H₄ bridge (Scheme 1). No ring current is visible for 2(2), only localized currents around the nitrogen atoms and on the carbon backbone. This shows that the ring current in the five-membered ring of 1(2) is not a mere consequence of the presence of a ring, but that the unsaturated C¹–C² bond actually plays an important role in the aromaticity/anti-aromaticity of the five-membered ring.

Now, a pattern emerges: the members of the series 1(n), with n being even, have an aromatic five-membered ring, whereas those with n being odd, have an anti-aromatic five membered ring. This is corroborated by plots of the π ring currents for the other members of the series (Fig. S1). For the compounds 1(n), with n = 4, 6, 8, and 10, a diatropic ring current is found, while the compounds 1(n), with n = 5, 7, 9, and 11, possess a paratropic ring current. This odd/even behaviour can be explained using the ipsocentric rules for induced current density,^17,18,20 in terms of orbital contributions and virtual translational (rotational) transitions that contribute to the diatropic (paratropic) ring current. For this analysis, we consider the π orbitals of the HC=CH moiety to interact with the π orbitals of the methine chain (Fig. 3). In C_2v, the HOMO and LUMO of the methine chain transform according to the a₂ and b₁ irreducible representations for even n, and according to the b₁ and a₂ representations for odd n. The HOMO and LUMO of the HC=CH moiety transform according to the b₁ and a₂ representations, respectively. Consequently, for even n, the HOMO of 1(n) is stabilized by the interaction of the methine chain (a₂ symmetry) with the anti-bonding, a₂, π orbital of the HC=CH moiety. The LUMO of b₁ symmetry of the methine chain interacts with the occupied b₁ π orbital of the HC=CH moiety, resulting in a destabilization of the LUMO of 1(n). The result of these interactions is a widening of the HOMO–LUMO gap in the series 1(n), with n even (Fig. 3).

Fig. 3 Orbital interaction diagram for n even and odd.

For the molecules of the series with odd n, the situation is reversed (Fig. 3): the methine chain HOMO of b₁ symmetry interacts with the occupied π orbital of the HC=CH moiety, leading to destabilization of the HOMO of 1(n), with n odd. The LUMO of a₂ symmetry interacts with the empty, a₂, π orbital of HC=CH, resulting in stabilization of the LUMO of 1(n). Hence, for n odd, the HOMO–LUMO is smaller, and consequently, the paratropic contribution to ring current arising from the virtual allowed rotational transition from the b₁ occupied orbitals to the empty a₂ LUMO is increased. This is further corroborated by a spectral decomposition of the ring current (Fig. 4a). If only the rotational transitions from the HOMO and HOMO−1 of 1(3) are taken into consideration, a strong, paratropic ring current is obtained in the five-membered ring. The absolute value of the mixing coefficient in the perturbation expansion for the paratropic HOMO → LUMO transition was calculated to be 7.4. The next largest mixing coefficient is that for the HOMO−1 → LUMO transition of 4.7, with other coefficients smaller than 3.5. This indicates that the dominant virtual transition is indeed the rotationally allowed HOMO → LUMO transition, together with the HOMO−1 → LUMO transition, but other virtual transitions contribute to the current as well, as shown in Fig. 4a. If all these other transitions (translational and rotational) are taken into consideration, except the two rotational HOMO → LUMO and HOMO−1 → LUMO transitions, a strong diatropic current is obtained, similar to the current in the 1(n) series with even n.

Fig. 4 (a) The contributions from the rotationally allowed HOMO → LUMO and HOMO−1 → LUMO transitions, and the contributions from all other translationally and rotationally allowed transitions to the π ring current of 1(3), and (b) the contributions of the HOMO → LUMO, HOMO−1 → LUMO, HOMO → LUMO+7, and HOMO−1 → LUMO+8 transitions, and the contributions from all other translationally and rotationally allowed transitions to the π ring current of 1(2).

It must be noted that the diatropic ring current in all 1(n) molecules is not determined by only a few translational transitions, but more transitions have to be taken into account to obtain the diatropic ring current. The spectral decomposition (Fig. 4b) shows that the ring current is dominated by the translational and rotational virtual transitions from the HOMO and HOMO−1 to the LUMO, LUMO+1, LUMO+7, and LUMO+8. The remaining virtual transitions do not significantly contribute to the diatropic ring current.

Further support for this alternation of the HOMO–LUMO gap is obtained from TDDFT calculations. The excitation energy to the first ¹B₂ state has been calculated for the series (Table 3). A decreasing trend in excitation energy is observed in the series 1(n), with n even, and, separately, a decreasing trend is observed for 1(n), with n odd. However, the excitation energies for the n even series are higher than those for n odd, in agreement with the model proposed. The excitation energy for 1(2) is higher than for 2(2), while the excitation energy for 1(3) is lower than that of 2(3), substantiating the hypothesis that the orbital interactions in the n even cases would increase the gap, whereas in the n odd cases it would decrease the gap.

Table 3 The excitation energies (E_exc, in eV) to the first ¹B₂ state, and the character (only weights higher than 10%) of the transition for 1(n), n = 2–11, and 2(n), n = 2, 3

Molecule	Eexc	Character	Molecule	Eexc	Character
1(2)	4.04	76.4% H−1 → L	1(3)	2.02	79.6% H → L
		23.9% H → L+1			19.3% H−1 → L
1(4)	3.02	91.5% H → L	1(5)	1.77	74.6% H → L
					23.8% H−1 → L
1(6)	2.44	96.5% H → L	1(7)	1.63	78.5% H → L
					19.0% H−2 → L
1(8)	2.07	98.1% H → L	1(9)	1.50	85.8% H → L
					10.7% H−2 → L
1(10)	1.81	98.4% H → L	1(11)	1.37	90.6% H → L
2(2)	3.76	98.0% H → L	2(3)	3.50	98.8% H →L

---

As has been noted earlier, at the UHF level, the closed shell solution becomes unstable, and a broken symmetry solution with lower energy is found. This symmetry breaking may have a consequence for the induced π ring current. The smallest molecule in the series for which a broken symmetry solution is found is 1(9). For the closed-shell solution, the ring current follows the general pattern observed for the molecules in the series with n odd. A paratropic ring current is thus found (Fig. 5). However, the broken symmetry solution is of b₂ symmetry, meaning that a b₁ and an a₂ orbital are singly occupied. In the spin-up (α) manifold, the highest occupied orbital has b₁ symmetry (and the lowest unoccupied a₂), whereas in the spin-down (β) manifold, the highest occupied orbital transforms according to the a₂ representation, with the lowest unoccupied β orbital transforming according to the b₁ irreducible representation. This means that the induced current density in the spin-up (α) channel is similar to that of closed-shell 1(9), which is a paratropic ring current. However, the consequence of the switch in occupation in the spin-down (β) channel, is that the ring current in this spin manifold follows the n even rules and results in diatropic contribution (Fig. 5). The sum of the α and β current densities lead to localized currents around the double bonds, and the five-membered ring becomes non aromatic (Fig. 5).

Fig. 5 Plots of the induced π current density for 1(9), 1-BS(9), the induced π current density for α electrons (1(9)-α), and the induced π current density for β electrons (1(9)-β).

Valence bond calculations

A remaining question that needs to be answered is whether this alternating aromaticity/anti-aromaticity with n even/odd can be explained on the basis of resonance structures. For this purpose, we performed valence bond calculations. In both cases, structures can be drawn with the positive charge delocalized over the methine chain (Scheme 2, structures A, B, and E). For n odd, additional structures C can be drawn, with the positive charge on the methine backbone, in the corner of the five-membered ring, and additional structures D, with the positive charge delocalized in the five-membered ring (note that these structures cannot be drawn for the n even series). Additionally, when drawing double ionic structures (the structures with two positive charges and one negative charge), for both series, the structure H can be drawn, with the negative charge on the top carbon of the five-membered ring. Only for the n even series, structures can be drawn with the negative charge delocalized in the five-membered ring (structures F and G). Hence, based on drawing resonance structures, it is possible to delocalize six π electrons in the five-membered ring of the n even series, but in the n odd series, there are four π electrons in the five-membered ring.

The calculated GN weights for the valence bond wavefunctions are summarized in Table 4. The non-aromatic structures A and B have the highest weight for both the molecules 1(2) and 1(3). For 1(2), structure E is also non-aromatic, whereas the structures F, G, and H, each with six π electrons in the five-membered ring, have a cumulative weight of 0.114. Hence, also based on this valence bond analysis, it can be concluded that the five-membered ring in 1(2), and in general in 1(n) with n even, has aromatic character. However, for n odd, the situation is different. In this case, structure H is non-aromatic, but structures C, D, and E have four π electrons in the five-membered ring, and contribute therefore to the anti-aromatic nature for this ring in the 1(n) series, with odd n. The cumulative weight of these anti-aromatic structures is 0.179, which is non-negligible, and these contributions give rise to the presence of the paratropic ring current in the n odd members of the series.

Table 4 The GN weights of the valence bond structures (depicted in Scheme 2) in the wavefunctions for 1(2) and 1(3)

Structure	1(2)	1(3)
2 × A	0.451	0.422
2 × B	0.334	0.344
2 × C	—	0.172
2 × D	—	0.001
1 × E	0.100	0.006
2 × F	0.092	0.047
2 × G	0.016	—
1 × H	0.006	0.007

---

Although the weights of the ionic aromatic and anti-aromatic structures is modest, it is important to notice that the magnetic response is dominated by mixing in excited states into the ground state via the angular (paratropic contribution) or the linear (diatropic contribution) operator.^17,18 These excited states consists of ionic structures,⁴⁶ and the accessibility of these excited states determines the magnetic response.

This observed odd/even effect for aromaticity may have further implications for the applicability of this class of molecules in for example electronic circuits. The conductance of molecules is dependent on the aromaticity of the five-membered ring,^47,48 and an odd/even effect is, hence, also expected for quantum interference effects and should thus be considered when designing single molecule circuits.

Conclusions

The five-membered ring in the cyanine dyes 1(n), with n = 2–11 show alternating aromatic and anti-aromatic character for n is even or odd, respectively. This odd/even effect has been explained using the ipsocentric symmetry selection rules. For the odd members of the series, the occupied b₁ π orbital of the HC=CH moiety interacts with the b₁ HOMO of the methine chain, while the unoccupied a₂ π orbital of the HC=CH moiety interacts with the LUMO of the methine chain. These interactions lead to a destabilization of the HOMO and a stabilization of the LUMO, thereby reducing the HOMO–LUMO gap. The paratropic contribution to the ring current is increased due to this lower energy gap, resulting in a final paratropic ring current, and hence, anti-aromatic character for the n odd members. A spectral decomposition analysis of the induced current density and valence bond calculations further corroborate this interpretation. The valence bond calculations show that aromatic resonance structures, with six π electrons in the five-membered ring, contribute to the wavefunction for the n even members of the series, whereas for the n odd members, anti-aromatic structures, with four π electrons in the five-membered ring, contribute to the wavefunction. This odd/even effect may have consequences for the optical properties and reactivity of the cyanine dyes.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: plots of the induced π current density for 1(n), n = 4–11 (Fig. S1). Cartesian coordinates of the optimized geometries for 1(n), n = 2–11 (Tables S1–S10). See DOI: https://doi.org/10.1039/d6cp00085a.

Acknowledgements

RWAH and NPM thank Kristin Becker (University of Groningen) for fruitful discussions. NPM acknowledges the Dieptestrategie Programme of the Zernike Institute for Advanced Materials (University of Groningen, the Netherlands) for financial support. We thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Hábrók high performance computing cluster. This work also made use of the Dutch national e-infrastructure with the support of the SURF cooperative using grants no. EINF-3508 and NWO-2022.004. This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy (DOE). The U.S. government retains, and the publisher, by accepting the article for publication, acknowledges that the U.S. government retains, a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/doe-public-access-plan).

References

1. G. T. Dempsey, M. Bates, W. E. Kowtoniuk, D. R. Liu, R. Y. Tsien and X. Zhuang, J. Am. Chem. Soc., 2009, 131, 18192–18193.
2. Y. Ooyama and Y. Harima, Eur. J. Org. Chem., 2009, 2903–2934.
3. G. Patonay, J. Salon, J. Sowell and L. Strekowski, Molecules, 2004, 9, 40–49.
4. A. A. Ishchenko, Russ. Chem. Rev., 1991, 60, 865–884.
5. K. Becker, N. P. Neme, T. L. C. Jansen, M. S. Pschenichnikov, J. C. Hummelen, R. W. A. Havenith and M. D. Witte, unpublished work.
6. K. Becker, Exploring the Terra Incognita of Dye-Sensitised Upconversion of NIR-II Light, PhD thesis University of Groningen, Groningen, The Netherlands, 2025.
7. N. P. Neme, Quantum chemistry and cyanine dye-sensitization: tuning photophysical properties for advanced applications, PhD thesis University of Groningen, Groningen, The Netherlands, 2025.
8. N. P. Neme, T. L. C. Jansen and R. W. A. Havenith, Phys. Chem. Chem. Phys., 2024, 26, 6235–6241.
9. P. von, R. Schleyer and H. Jiao, Pure Appl. Chem., 1996, 68, 209–218.
10. F. London, J. Phys. Radium, 1937, 8, 397–409.
11. L. Pauling, J. Chem. Phys., 1936, 4, 673–677.
12. J. A. Pople, J. Chem. Phys., 1956, 24, 1111.
13. J. A. Elvidge and L. M. Jackman, J. Chem. Soc., 1961, 859–866.
14. J. A. N. F. Gomes and R. B. Mallion, Chem. Rev., 2001, 101, 1349–1384.
15. S. Coriani, P. Lazzeretti, M. Malagoli and R. Zanasi, Theor. Chim. Acta, 1994, 89, 181–192.
16. T. Keith and R. F. W. Bader, Chem. Phys. Lett., 1993, 210, 223–231.
17. E. Steiner and P. W. Fowler, Chem. Commun., 2001, 2220–2221.
18. E. Steiner and P. W. Fowler, J. Phys. Chem. A, 2001, 105, 9553–9562.
19. E. Steiner and P. W. Fowler, Int. J. Quantum Chem., 1996, 60, 609–616.
20. E. Steiner and P. W. Fowler, Phys. Chem. Chem. Phys., 2004, 6, 261–272.
21. E. Steiner, A. Soncini and P. W. Fowler, J. Phys. Chem. A, 2006, 110, 12882–12886.
22. S. Shaik and P. C. Hiberty, A Chemist's guide to Valence Bond Theory, John Wiley & Sons, Inc., Hoboken, New Jersey, 2008.
23. D. L. Cooper, J. Gerratt and M. Raimondi, in Advances in Chemical Physics: Ab initio Methods in Quantum Chemistry - II, ed. K. P. Lawley, John Wiley & Sons, London, 1987, vol. LXIX, pp. 319–397.
24. M. Zielinski, R. W. A. Havenith, L. W. Jenneskens and J. H. van Lenthe, Theor. Chem. Acc., 2010, 127, 19–25.
25. Z. Rashid, J. H. van Lenthe and R. W. A. Havenith, J. Phys. Chem. A, 2012, 116, 4778–4788.
26. C. Adama and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
27. M. Ernzerhof and G. E. Scuseria, J. Chem. Phys., 1999, 110, 5029–5036.
28. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
29. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 77, 1396–1396.
30. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
31. M. F. Guest, I. J. Bush, H. J. J. van Dam, P. Sherwood, J. M. H. Thomas, J. H. van Lenthe, R. W. A. Havenith and J. Kendrick, Mol. Phys., 2005, 103, 719–747.
32. R. E. Stratmann, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 1998, 109, 8218–8224.
33. C. A. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications, Oxford University Press, New York, 5th edn, 2012.
34. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Rev. C.01.
35. P. Lazzeretti, M. Malagoli and R. Zanasi, Chem. Phys. Lett., 1994, 220, 299–304.
36. R. Zanasi, J. Chem. Phys., 1996, 105, 1460–1469.
37. R. W. A. Havenith and P. W. Fowler, Chem. Phys. Lett., 2007, 449, 347–353.
38. P. Lazzeretti and R. Zanasi, SYSMO package (University of Modena), 1980. Additional routines by P. W. Fowler, E. Steiner, R. W. A. Havenith, A. Soncini.
39. A. Soncini, J. Chem. Theory Comput., 2007, 3, 2243–2257.
40. A. Soncini and P. W. Fowler, Chem. Phys. Lett., 2008, 450, 431–436.
41. J. H. van Lenthe, F. Dijkstra and R. W. A. Havenith, in Valence Bond Theory, ed. D. L. Cooper, Elsevier, Amsterdam, 2002, vol. 10, pp. 79–112.
42. J. Verbeek, J. H. Langenberg, C. P. Byrman, F. Dijkstra, R. W. A. Havenith, J. J. Engelberts, M. Zielinski, Z. Rashid and J. H. van Lenthe, TURTLE, an ab initio VB/VBSCF program, Utrecht, The Netherlands, 2026.
43. J. H. van Lenthe and G. G. Balint-Kurti, Chem. Phys. Lett., 1980, 76, 138–142.
44. J. H. van Lenthe and G. G. Balint-Kurti, J. Chem. Phys., 1983, 78, 5699–5713.
45. G. A. Gallup and J. M. Norbeck, Chem. Phys. Lett., 1973, 21, 495–500.
46. R. W. A. Havenith, J. Org. Chem., 2006, 71, 3559–3563.
47. N. P. Arasu and H. Vázquez, Chem. Phys. Chem., 2021, 22, 864–869.
48. A. Borges and G. C. Solomon, J. Phys. Chem. C, 2017, 121, 8272–8279.

---