Next Article in Journal
Broadband Solar Absorber and Thermal Emitter Based on Single-Layer Molybdenum Disulfide
Next Article in Special Issue
Ab Initio Neural Network Potential Energy Surface and Quantum Dynamics Calculations on Na(2S) + H2 → NaH + H Reaction
Previous Article in Journal
Cu-ZnO Embedded in a Polydopamine Shell for the Generation of Antibacterial Surgical Face Masks
Previous Article in Special Issue
Mono- and Binuclear Complexes in a Centrifuge-Less Cloud-Point Extraction System for the Spectrophotometric Determination of Zinc(II)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Simulation of Photoreactions Using Restricted Open-Shell Kohn–Sham Theory

Theoretical Chemistry, Leibniz University Hannover, Callinstr. 3A, 30167 Hannover, Germany
*
Author to whom correspondence should be addressed.
Molecules 2024, 29(18), 4509; https://doi.org/10.3390/molecules29184509
Submission received: 29 July 2024 / Revised: 18 September 2024 / Accepted: 20 September 2024 / Published: 23 September 2024

Abstract

:
It is a well-established standard to describe ground-state chemical reactions at an ab initio level of multi-electron theory. Fast reactions can be directly simulated. The most widely used approach is density functional theory for the electronic structure in combination with molecular dynamics for the nuclear motion. This approach is known as ab initio molecular dynamics. In contrast, the simulation of excited-state reactions at this level of theory is significantly more difficult. It turns out that the self-consistent solution of the Kohn–Sham equations is not easily reached in excited-state simulations. The first program that solved this problem was the Car–Parrinello molecular dynamics code, using restricted open-shell Kohn–Sham theory. Meanwhile, there are alternatives, most prominently the Q-Chem code, which widens the range of applications. The present study investigates the suitability of both codes for the molecular dynamics simulation of excited-state motion and presents applications to photoreactions.

1. Introduction

The excited-state simulation of molecular systems is significantly more difficult than the ground-state simulation. The excited-state dynamics decides about the outcome of photoreactions. Generally, a system in its singlet ground state S0 is excited by light in the ultraviolet/visible regime (UV/VIS) and reaches a singlet excited state, whereby the photon energy decides which excited state is populated. Such an excitation is assumed to be vertical in an energy diagram, in the sense that it is so fast that the nuclei hardly move during the excitation. At this point, Kasha’s empirical rule [1] applies: Photochemistry occurs from the first excited state only. Higher excitations decay immediately, or, more precisely, within a few femtoseconds to the singlet S1 state. This time span is generally not enough for a photoreaction to occur because the nuclei move too slowly. Like this, upon excitation, the so-called Franck–Condon region of the S1 state is reached within a few femtoseconds. In the Franck-Condon region, the system is in the electronically excited state, but still has the geometry of the ground state. It was the success of Robb and Olivucci [2,3] on the basis of CASSCF (complete active space self-consistent-field) calculations, who showed that the motion to the product state is possible via conical intersections. At these points, the excited state and ground state touch, which facilitates the return to the ground state. It is, meanwhile, clear that conical intersections are indeed omnipresent. However, it turns out, that, in contrast to the claim by Robb and Olivucci, it is rather the motion in the Franck–Condon region, not the motion through the particular conical intersection that decides which photoproduct is formed. An example is methyl chloride: the question if a C-Cl or a C-H bond is broken is decided immediately after the excitation. The system moves either along the C-Cl or along the C-H coordinate. Along both coordinates it will find a conical intersection, which leads to the ground state [4].
Note, however, that the motion in the Franck–Condon region is not the only factor that decides the outcome of a photoreaction. The interplay of a variety of factors ultimately leads to product formation [5].
While it is possible with CASSCF calculations to locate conical intersections, it is not so easy to describe the dynamics through these important points of the potential energy surface. Empirically, these points are passed diabatically. If one uses an adiabatic method like CASSCF, one will want to diabatize the potential energy surfaces. There are several possible approaches. These, however, often do not yield satisfying results in a straight-forward way; hence, the surface hopping approach by Tully [6] is often applied, which, as an empirical approach, represents quite a compromise in an ab initio calculation. Many alternative approaches have been developed over the years, for an overview, see [7,8,9]. Most recently, non-adiabatic dynamics has been combined with machine learning to accelerate the dynamics and to make the simulation of larger systems feasible [10].
For photoreactions that occur in the first excited singlet state, there is an efficient alternative that is fully ab initio: Restricted open-shell Kohn–Sham (ROKS) theory directly yields diabatic surfaces. This is due to the fact that it is a single-configuration method. At the same time, it is a two-determinant method [11]. Two determinants are needed in open-shell calculations to obtain the correct spin symmetry. An exception is the highest-spin case, where a single determinant is enough.
ROKS has a long history. The term was probably first used by Pople, Gill, and Handy [12], who stated that instead of ROKS, an unrestricted formulation should be chosen, without giving a clear reason for this statement. In 1988, we published the relevant equations for the first time and directly implemented them in the CPMD code, in order to simulate photoreactions [13]. This was successful for the simple cases we started from, namely n → π situations. We called the method low spin excitation (LSE) or restricted open-shell density functional theory (RODFT) at the time. Shortly thereafter, Filatov and Shaik published the same method and called it ROKS [14]. They already addressed the problem of orbital rotations during excited-state SCF simulations (see below). A very similar development is the ROSS approach (restricted open-shell singlet state) [15]. In this method, the homogeneous electron gas approximation is applied in the end of the derivation. Due to error compensation, ROSS in combination with GGA functionals provides slightly better results than ROKS/GGA. However, ROKS can be more easily generalized to arbitrary open-shell situations and, hence, is the preferable method [11]. The general picture resulting from the sum rule for up to five unpaired electrons is depicted in Figure 1.
It turns out that every configuration can be composed from two Slater determinants. The general formula for N unpaired electrons [11] reads as follows:
E j C = N + 1 + M 2 M E j S D N + 1 M 2 M E j 1 S D
Hereby, N is the number of unpaired electrons, M is the multiplicity, and j denotes the jth lowest configuration or determinant. For this energy expression, ROKS operators can be obtained by functional variation [11]. In the high-spin case, the second term vanishes.
ROKS has several drawbacks, one of them being the orbital rotation problem. Naive minimization of the low-spin state results in convergence to the high-spin state, with the HOMO and LUMO forming a linear combination that does not reflect the molecular symmetry anymore. This linear combination is possible if HOMO and LUMO have the same symmetry, e.g., in π - π transitions, not in n - π transitions. In 1993, we presented a solution for the orbital rotation problem that works for many cases [16], based on the work by Hirao and Nakatsuji [17]. Alternatively, the algorithm can be viewed as a modification of the procedure proposed by Goedecker and Umrigar [18]. We did not succeed in combining this approach with Car–Parrinello molecular dynamics, because the orbital orthogonality problem became very difficult to solve. Hence, Born–Oppenheimer molecular dynamics must be used with ROKS, unless the resulting error is small. A small error is observed for a small exchange interaction between HOMO and LUMO.
There are other approaches in the literature addressing the orbital rotation problem [19,20]. To find out which are the most general and stable ones, one must compare practically existing codes. Ideally, an excited-state calculation must be self-consistent and should be a black-box calculation as far as possible.
It should be mentioned that the orbital rotation problem is not unique to ROKS. From what we know, the same or similar problems arise in all excited-state SCF calculations. In CASSCF calculations, the problem is solved by using state-average CASSCF. By mixing with the ground state, the excited-state orbitals cannot escape to an unphysical solution. Despite these problems, it is highly desirable to achieve an excited-state SCF when describing photoreactions. Not only is it closer to the experimental picture than a repeated excitation from the ground state, an excited-state SCF has also the advantage that the Hellmann–Feynman theorem can be used for the calculation of the gradients. Hence, one has immediately analytic gradients available that can be used in excited-state geometry optimizations and, in particular, in excited-state molecular dynamics [13].
ROKS must not be confused with Δ -SCF, which deals with single determinants, making it far easier to implement, but less accurate. As has been shown by Hait and Head-Gordon, recently, that Δ -SCF and ROKS can be summarized under the more general expression orbital-optimized DFT (OO-DFT) [20].
For ROKS/BLYP a significant red-shift is observed, rendering the Franck–Condon region too low in energy and reducing the kinetic energy that is set free during a photoreaction. Nevertheless, it was possible to simulate quite a few interesting photoreactions, for example the rhodopsin isomerization [21] and Feringa’s nanorotor [22].
Despite these successes, ROKS was always viewed as a cheap and dirty alternative to the in- principle exact TDDFT approach. However, more recently, ROKS has attracted renewed attention, based on the statement that ROKS is superior to the alternative method of time-dependent density functional theory (TDDFT) in many respects [23]. It should be emphasized that this is only true for the treatment of photoreactions. For the computation of electronic spectra, ROKS is unsuited and TDDFT is the method of choice. TDDFT in its original meaning [24] allows, in principle, for computing the time development of the electronic structure. In practice, it is dominantly used in connection with time-dependent perturbation theory in order to compute electronic spectra [25]. The result may be called linear-response TDDFT (LR-TDDFT). However, as almost all TDDFT calculations performed nowadays are LR-TDDFT calculations, the abbreviated term TDDFT is normally used for the computation of spectra. In contrast with ROKS, TDDFT yields the complete spectrum at once, not only the lowest excited states. It has, however, certain drawbacks. TDDFT fails to describe charge-transfer situations, while these are correctly described by ROKS, where the occupation of LUMO leads to the correct result. More important is the fact that time-dependent perturbation theory fails when two states have the same energy, because the energy difference is in the denominator of the TDDFT energy expression. Hence, exactly at the conical intersections, TDDFT tends to fail seriously. Extensions of TDDFT perform better [26,27,28], but their usefulness in molecular dynamics simulations of photochemistry, where several thousand points must be computed self-consistently, must still be demonstrated.
At this point, we will discuss the shortcomings of ROKS. Besides the optimization problem, there is the red-shift of typically 0.6 eV, which slows down reactions. As the Hartree–Fock variant ROHF is normally too high in energy by about 1.0 eV, a mixture of GGA exchange and exact exchange like in hybrid functionals should yield better results. Unfortunately, it is difficult at best to conduct ROKS/B3LYP calculations with the CPMD or the GAUSSIAN codes. The Q-Chem code might represent the solution for at least some of the problems discussed above. It is the aim of the present study to test the Q-Chem implementation in comparison with the CPMD code in order to explore its suitability for the simulation of photoreactions. Unlike previous large benchmark studies on other excited-state methods [29], where only vertical transitions were investigated, we focus not only on vertical transitions but also on the shape of the potential energy surfaces by exploring the geometry changes upon excitation, as reflected by the adiabatic transitions. In addition, we studied the excited-state dynamics that cannot be fully automated and that are expensive in terms of CPU time, and are thus performed for selected systems only.
For simplicity, we did not address internal conversion or intersystem crossing in the present study, even if such processes matter in photoreactions [30]. We simply investigated the dynamics in the first few picoseconds after the system has reached the first excited state.

2. Results and Discussions

2.1. Excitation Energies and Orbitals

The vertical excitation energies of CPMD and Q-Chem are in good agreement (Table 1 and Figure 2), with Q-Chem being a bit higher in energy (0.2 eV on average). This is attributed to the different basis sets. Q-Chem with its smaller Gaussian basis sets produces higher excitation energies; that is, we observe error compensation. As expected, B3LYP cures the red shift of BLYP, but not as strongly as expected (0.2 eV on average). TDDFT is again blue-shifted by 0.2 eV compared with ROKS, and is, hence, in better agreement with experiment for these vertical excitation energies. A clear deviation of roughly 1.5 eV between some CPMD and Q-Chem results is observed for the adiabatic excitations. This deviation is attributed to a failure of Q-Chem to converge to the right SCF solution, and will be investigated more deeply in the next chapter.
The orbitals plotted with Q-Chem (Figure 3) demonstrate that the HOMOs and LUMOs are properly computed.

2.2. Molecular Dynamics

2.2.1. Ethene

Molecular dynamics is an important and sensitive check for the stability of a wavefunction. Only a stable self-consistent field (SCF) calculation will lead to stable dynamics; that is, to conservation of the total energy during a reaction. To understand the differences between CPMD and Q-Chem, we performed excited-state dynamics simulations for ethene. The result obtained with the CPMD code is shown in Figure 4. Born–Oppenheimer molecular dynamics and Car–Parrinello molecular dynamics as computed with the CPMD code yield similar dynamics. In contrast, Q-Chem, where we can use Born–Oppenheimer dynamics only, fails to converge after a few steps. This is also true if non-default optimization algorithms are employed, like level-shifting. It indicates that Q-Chem has a problem with rotating the singly occupied orbitals properly in non-planar π systems, leading to too low adiabatic energies. The modified Goedecker algorithm [16,18] as implemented in the CPMD code performs better and permits the computation of several thousand excited-state SCF calculations, as afforded in Born–Oppenheimer molecular dynamics runs.
Movie files of the ethene isomerization are contained in the Supplementary Information together with input files. The movies demonstrate that BOMD as implemented in the CPMD code is the method of choice. The development of the orbitals with time is most convincing for BOMD. With Q-Chem, it is presently not possible to plot the orbitals during a molecular dynamics run.

2.2.2. Toluol and Chlorine

Students in Germany are taught the elementary rule “Sonne, Siedehitze, Seitenkette” (“sun, heat, side chain”), meaning that the photoreaction of aromates with halides occurs preferentially at elevated temperatures and at the side chain. An example is the reaction of toluol with chlorine (Figure 5). In the ROKS simulation, we observe the immediate photodissociation of chlorine, both with CPMD and Q-Chem (Figure 6). The dissociation leads diabatically to the open-shell ground state. At this point, it is possible to switch to the unrestricted Kohn–Sham (UKS) formalism to continue the calculation in the ground state. With Q-Chem, the resulting radicals and the toluol molecules dissociate from each other, and the periodic boundary conditions of the CPMD code guarantee that the fragments encounter each other again. This makes the CPMD code clearly superior for the simulation of reactions involving several molecules as well as reactions involving a solvent.
The next step, namely the attack of chlorine radicals at the side chain, is neither observed in the excited state nor after return to the ground state. The potential is shallow, but the reaction entropy prevents a reaction on the picosecond timescale. The final step, namely the radical recombination that leads to product formation, is observed spontaneously again in unrestricted ground-state CPMD simulations.
Movies of the chlorine dissociation are contained in the Supplementary Information. They reveal a disadvantage of Born–Oppenheimer molecular dynamics simulations—the calculation may oscillate between several SCF solutions. In a movie, one observes a flickering. Fortunately, such flickering is observed only rarely. Otherwise, CPMD and BOMD yield a similar picture.

2.2.3. Silver Citrate

Silver citrate is used as the basis for the formation of nanostructured silver coatings [31,32]. Upon exposure to light, it displays a fast photoreaction, leading to the formation of three CO2 molecules (Figure 7). The CPMD and Q-Chem simulations show that the reaction is not immediate in contrast with the chlorine dissociation (Figure 8). The bond breaking is caused by heat rather than by a specific anti-binding interaction. As a consequence, bond breaking occurs after some time in the excited state and may also be observed after return to the ground state. In this sense, the reaction is autocatalytic—the high amount of energy that is set free may cause further C–C bond breaking reactions, leading to CO2 formation.
Interestingly, C–C bond breaking is often connected with a hydrogen transfer reaction (Figure 9 and Table 2).
We find a tendency to a higher reactivity with the B3LYP functional—here, two C–C bonds break, while we observe the breaking of just one C–C bond with the BLYP functional. This is in line with the higher excitation energy obtained with B3LYP.
Movies of the silver citrate photoreaction are contained in the Supplementary Information. Again, CPMD and BOMD performed with the CPMD code yield a similar picture.

3. Methods

Geometry optimizations and single point calculations were performed using the CPMD code, version 4.3 [33,34,35], with the Becke–Lee–Yang–Parr (BLYP) functional [36,37]. In the toluol and citrate simulations, the Grimme dispersion correction was applied [38]. The time step was chosen as 2 a.u. (0.0484 fs) for both CPMD and BOMD simulations, in order to minimize the energy drift in the BOMD simulations. The fictitious electron mass used by the CP dynamics was chosen as 200 a.u. in the CPMD simulations. Troullier–Martins pseudopotentials, as optimized for the BLYP functional, were employed for describing the core electrons [39,40]. The plane-wave cutoff, which determines the size of the basis set, was set to 70.0 Rydberg. Depending on the size of the molecules, orthorhombic simulation cells were used with extensions ranging from 15 × 15 × 15 a . u . 3 (7.94 × 7.94 × 7.94 Å3) up to 20 × 15 × 15 a . u . 3 (10.58 × 10.58 × 10.58 Å3) and 18 × 18 × 18 a . u . 3 (9.53 × 9.53 × 9.53 Å3). The system sizes employed for the toluol/chlorine and silver citrate systems are given in Table 2. The systems were equilibrated in the ground state at a certain temperature, then let go freely in the NVE ensemble. After simulation in the ground state, the systems were vertically placed in the excited state; that is, by keeping the ground-state coordinates and velocities.
Table 2 lists the parameters used in the single simulation runs and timings.
Q-Chem calculations using the BLYP and B3LYP [41] functionals were performed with version 6.0 [42]. 6-311G** and LANL2DZ basis sets were used. A time step of 10 a.u. was chosen in the Born–Oppenheimer molecular dynamics calculations.

4. Conclusions

The CPMD code was clearly ahead of its time for several decades. This mainly refers to the ability to perform both Car–Parrinello molecular dynamics and Born–Oppenheimer molecular dynamics. Due to full parallelization and the use of FORTRAN libraries like BLAS and LAPACK, the CPMD code makes optimal use of the CPU time provided by a linux cluster, and is hence favoured by computation centers. Meanwhile, other codes like Q-Chem become competitive. This is particularly true for the availability of higher-level functionals like B3LYP. For the examples we investigated, implementation in the CPMD code was proven to be more stable for molecular dynamics simulations. The CPMD code was explicitly developed for performing molecular dynamics right from the beginning, in contrast with the Q-Chem code. This also resulted in a more convenient format of the relevant output data in CPMD and better restart options.
With both codes, it was possible to simulate the photoreaction of silver citrate, which, to the best of our knowledge, has not been simulated before at this level of theory. The simulation of the reaction of toluol with chlorine is more straight-forward with the CPMD code because it is a bimolecular reaction. A viable repair of the problem, that molecules fly away from each other unless periodic boundary conditions are used, can be achieved by restraints or repulsive potentials that are implemented in several program codes. An example is the nanoreactor by Martínez et al. [43].
The desired result of reactive simulations of photoreactions, namely movies with orbitals, is presently obtained with the CPMD code only. Making these movies is important not only for illustration, but also for judging if there are any unphysical jumps.
It should be emphasized at this point that excited-state SCF simulations are never perfect black-box calculations, and one must not expect that every run will be successful. When doing BOMD simulations with the CPMD code, the ODIIS 2 or PCG MINIMIZE options greatly improve the SCF convergence.
Future work will also include comparison to the CP2k code, which is in several respects similar to Q-Chem (Gaussian basis functions, BOMD only).
As Q-Chem continues to be developed, it will soon outperform the CPMD code for several applications. In particular, it can be expected that machine learning will lead to optimized functionals that will not be available in the CPMD code. For the simulation of chemical reactions in the condensed phase, however, the CPMD code might stay the best program code for some time still. At present, artificial intelligence plays a minor role in the simulation of rare events like chemical reactions. On a longer time scale, it might be used to analyse reactive trajectories generated with the traditional program codes and to recognize reaction patterns. This becomes important with growing CPU time and growing complexity of the systems under investigation. Nevertheless, the traditional program codes will stay at the core of any simulation of chemical reactions and it is important to maintain them.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules29184509/s1, suppl.zip: sample CPMD version 4.3 input files ethene_CPMD_in1 (equilibration), ethene_CPMD_in2 (production, ground state), ethene_ CPMD_in3 (production, excited state), sample Q-Chem input file ethene_qchem_in1, movies for the ethene, toluol/chlorine, and silver citrate photoreactions (Q-Chem/BOMD, CPMD/BOMD and CPMD/CPMD).

Author Contributions

Conceptualization, I.F.; methodology I.F.; software, I.F.; validation, R.B. and I.F.; formal analysis, R.B., J.G., D.M. and I.F.; investigation, R.B., J.G., D.M. and I.F.; resources, I.F.; data curation, R.B. and I.F.; writing—original draft preparation, I.F.; writing—review and editing, R.B., L.Á. and I.F.; visualization, R.B., L.Á. and I.F.; supervision, L.Á. and I.F.; project administration, I.F.; funding acquisition, I.F. All authors have read and agreed to the published version of the manuscript.

Funding

Part of the calculations were performed on the local cluster of the Leibniz University of Hannover at the LUIS (Leibniz University IT Services), project number 424969120, and on the Höchstleistungsrechner Nord, HLRN, maintained by the North German Supercomputing Alliance, project nic00086.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article and Supplementary Materials.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
AIMDAb-Initio Molecular Dynamics
BLYPBecke-Lee-Yang-Parr density functional
BOMDBorn-Oppenheimer Molecular Dynamics
CASSCFComplete Active Space Self-Consistent Field
CPMDCar-Parrinello Molecular Dynamics
GGAGeneralized Gradient Approximation
SCFSelf-Consistent Field

References

  1. Kasha, M. Characterization of electronic transitions in complex molecules. Discuss. Farad. Soc. 1950, 9, 14. [Google Scholar] [CrossRef]
  2. Bernardi, F.; De, S.; Olivucci, M.; Robb, M.A. The Mechanism of Ground-State-Forbidden Photochemical Pericyclic Reactions: Evidence for Real Conical Intersections. J. Am. Chem. Soc. 1990, 112, 1737. [Google Scholar] [CrossRef]
  3. Bernardi, F.; Olivucci, M.; Robb, M.A. Potential energy surface crossings in organic photochemistry. Chem. Soc. Rev. 1996, 25, 321. [Google Scholar] [CrossRef]
  4. Frank, I. Classical motion of the nuclei in a molecule: A concept without alternatives. Chem. Sel. 2020, 5, 1872. [Google Scholar] [CrossRef]
  5. Hammes-Schiffer, S.; Stuchebrukhov, A.A. Theory of Coupled Electron and Proton Transfer Reactions. Chem. Rev. 2010, 110, 6939. [Google Scholar] [CrossRef] [PubMed]
  6. Tully, J.C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061. [Google Scholar] [CrossRef]
  7. Domcke, W.; Yarkony, D.R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325. [Google Scholar] [CrossRef]
  8. Curchod, B.F.E.; Martínez, T.J. Ab Initio Nonadiabatic Quantum Molecular Dynamics. Chem. Rev. 2018, 118, 3305. [Google Scholar] [CrossRef]
  9. Nelson, T.R.; White, A.J.; Bjorgaard, J.A.; Sifain, A.E.; Zhang, Y.; Nebgen, B.; Fernandez-Alberti, S.; Mozyrsky, D.; Roitberg, A.E.; Tretiak, S. Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials. Chem. Rev. 2020, 120, 2215. [Google Scholar] [CrossRef]
  10. Westermayr, J.; Gastegger, M.; Vörös, D.; Panzenboeck, L.; Joerg, F.; González, L.; Marquetand, P. Deep learning study of tyrosine reveals that roaming can lead to photodamage. Nat. Chem. 2022, 14, 914. [Google Scholar] [CrossRef]
  11. Schulte, M.; Frank, I. Restricted open-shell Kohn-Sham theory: N unpaired electrons. Chem. Phys. 2010, 373, 283. [Google Scholar] [CrossRef]
  12. Pople, J.A.; Gill, P.M.W.; Handy, N.C. Spin-Unrestricted Character of Kohn-Sham Orbitals for Open-Shell Systems. Int. J. Quantum Chem. 1995, 56, 303. [Google Scholar] [CrossRef]
  13. Frank, I.; Hutter, J.; Marx, D.; Parrinello, M. Molecular dynamics in low-spin excited states. J. Chem. Phys. 1998, 108, 4060. [Google Scholar] [CrossRef]
  14. Filatov, M.; Shaik, S. Spin-restricted density functional approach to the open-shell problem. Chem. Phys. Lett. 1998, 288, 689. [Google Scholar] [CrossRef]
  15. Gräfenstein, J.; Cremer, D.; Kraka, E. Density functional theory for open-shell singlet biradicals. Chem. Phys. Lett. 1998, 288, 593. [Google Scholar] [CrossRef]
  16. Grimm, S.; Nonnenberg, C.; Frank, I. Restricted open-shell Kohn-Sham theory for π-π* transitions. I. Polyenes, cyanines, and protonated imines. J. Chem. Phys. 2003, 119, 11574. [Google Scholar] [CrossRef]
  17. Hirao, K.; Nakatsuji, H. General SCF operator satisfying correct variational condition. J. Chem. Phys. 1973, 59, 1457. [Google Scholar] [CrossRef]
  18. Goedecker, S.; Umrigar, C.J. Critical assessment of the self-interaction-corrected-local-density-functional method and its algorithmic implementation. Phys. Rev. A 1997, 55, 1765. [Google Scholar] [CrossRef]
  19. Filatov, M. Spin-restricted ensemble-referenced Kohn–Sham method: Basic principles and application to strongly correlated ground and excited states of molecules. WIREs Comput. Mol. Sci. 2014, 5, 146. [Google Scholar] [CrossRef]
  20. Hait, D.; Head-Gordon, M. Orbital Optimized Density Functional Theory for Electronic Excited States. J. Phys. Chem. Lett. 2021, 12, 4517. [Google Scholar] [CrossRef]
  21. Röhrig, U.; Guidoni, L.; Laio, A.; Frank, I.; Röthlisberger, U. A molecular spring for vision. J. Am. Chem. Soc. 2004, 126, 15328. [Google Scholar] [CrossRef] [PubMed]
  22. Grimm, S.; Bräuchle, C.; Frank, I. Light-driven unidirectional rotation in a molecule: ROKS simulation. Chem. Phys. Chem. 2005, 6, 1943. [Google Scholar] [CrossRef] [PubMed]
  23. Kunze, L.; Hansen, A.; Grimme, S.; Mewes, J.M. PCM-ROKS for the description of charge-transfer states in solution: Singlet-triplet gaps with chemical accuracy from open-shell Kohn-Sham reaction-field calculations. J. Phys. Chem. Lett. 2021, 12, 8470. [Google Scholar] [CrossRef] [PubMed]
  24. Runge, E.; Gross, E.K.U. Density-functional theory for time-dependent systems. Phys. Rev. 1984, 52, 997. [Google Scholar] [CrossRef]
  25. Casida, M.E.; Salahub, D.R. Asymptotic correction approach to improving approximate exchange-correlation potentials: Time-dependent density-functional theory calculations of molecular excitation spectra. J. Chem. Phys. 2000, 113, 8918. [Google Scholar] [CrossRef]
  26. Shao, Y.; Head-Gordon, M.; Krylov, A.I. The spin–flip approach within time-dependent density functional theory: Theory and applications to diradicals. J. Chem. Phys. 2003, 118, 4807. [Google Scholar] [CrossRef]
  27. Wang, F.; Ziegler, T. Time-dependent density functional theory based on a noncollinear formulation of the exchange-correlation potential. J. Chem. Phys. 2004, 121, 12191. [Google Scholar] [CrossRef]
  28. Lee, S.; Filatov, M.; Lee, S.; Choi, C.H. Eliminating spin-contamination of spin-flip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix. J. Chem. Phys. 2018, 149, 104101. [Google Scholar] [CrossRef]
  29. Jacquemin, D.; Wathelet, V.; Perpete, E.A.; Adamo, C. Extensive TD-DFT Benchmark: Singlet-Excited States of Organic Molecules. J. Chem. Theory Comput. 2009, 5, 2420. [Google Scholar] [CrossRef]
  30. Marian, C. Spin–orbit coupling and intersystem crossing in molecules. WIREs Comput. Mol. Sci. 2012, 2, 187. [Google Scholar] [CrossRef]
  31. Zhang, Q.; Li, N.; Goebl, J.; Lu, Z.; Yin, Y. A Systematic Study of the Synthesis of Silver Nanoplates: Is Citrate a ”Magic” Reagent? J. Am. Chem. Soc. 2011, 133, 18931. [Google Scholar] [CrossRef] [PubMed]
  32. Bastus, N.G.; Merkoci, F.; Piella, J.; Puntes, V. Synthesis of Highly Monodisperse Citrate-Stabilized Silver Nanoparticles of up to 200 nm: Kinetic Control and Catalytic Properties. Chem. Mater. 2014, 26, 2836. [Google Scholar] [CrossRef]
  33. Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471–2474. [Google Scholar] [CrossRef] [PubMed]
  34. Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: Cambridge, UK, 2009. [Google Scholar]
  35. Hutter, J.; Alavi, A.; Deutsch, T.; Bernasconi, M.; Goedecker, S.; Marx, D.; Tuckerman, M.; Parrinello, M. CPMD Version 4.3. Copyright IBM Corp 1990–2015, Copyright MPI für Festkörperforschung Stuttgart 1997–2001. Available online: https://github.com/CPMD-code/CPMD/releases/tag/4.3 (accessed on 18 September 2024).
  36. Becke, A. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098. [Google Scholar] [CrossRef] [PubMed]
  37. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle–Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. [Google Scholar] [CrossRef]
  38. Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463. [Google Scholar] [CrossRef]
  39. Troullier, N.; Martins, J.L. Efficient Pseudopotentials for Plane-Wave Calculations. Phys. Rev. B 1991, 43, 1993. [Google Scholar] [CrossRef]
  40. Boero, M.; Parrinello, M.; Terakura, K.; Weiss, H. Car-Parrinello study of Ziegler-Natta heterogeneous catalysis: Stability and destabilization problems of the active site models. Mol. Phys. 2002, 100, 2935–2940. [Google Scholar] [CrossRef]
  41. Becke, A. A New Mixing of Hartree-Fock and Local-Density Functional Theories. J. Chem. Phys. 1993, 98, 1372. [Google Scholar] [CrossRef]
  42. Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package. J. Chem. Phys. 2021, 155, 084801. [CrossRef]
  43. Wang, L.P.; Titov, A.; McGibbon, R.; Liu, F.; Pande, V.S.; Martínez, T.J. Discovering chemistry with an ab initio nanoreactor. Nat. Chem. 2014, 6, 1044. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Energy diagram for up to five unpaired electrons [11]. Microstates as described by Slater determinants (SD) are linearly combined to states described by spin-adapted configurations (C).
Figure 1. Energy diagram for up to five unpaired electrons [11]. Microstates as described by Slater determinants (SD) are linearly combined to states described by spin-adapted configurations (C).
Molecules 29 04509 g001
Figure 2. Values of Table 1 plotted against each other.
Figure 2. Values of Table 1 plotted against each other.
Molecules 29 04509 g002
Figure 3. Orbitals of acetone as computed with Q-Chem (from top to bottom: ground state, HOMO; vertically excited state, HOMO; ground state, LUMO, vertically excited state, LUMO).
Figure 3. Orbitals of acetone as computed with Q-Chem (from top to bottom: ground state, HOMO; vertically excited state, HOMO; ground state, LUMO, vertically excited state, LUMO).
Molecules 29 04509 g003
Figure 4. Dihedral angle defining the ethene photoreaction, namely the rotation about the central double bond. Black: BOMD, ground state, red: BOMD, excited state, green: CPMD, ground state, blue: CPMD, excited state. Both BOMD and CPMD calculations were performed with the CPMD code.
Figure 4. Dihedral angle defining the ethene photoreaction, namely the rotation about the central double bond. Black: BOMD, ground state, red: BOMD, excited state, green: CPMD, ground state, blue: CPMD, excited state. Both BOMD and CPMD calculations were performed with the CPMD code.
Molecules 29 04509 g004
Figure 5. The reaction of toluol and chlorine. In the excited state, chlorine dissociates. This is followed by ground-state reactions.
Figure 5. The reaction of toluol and chlorine. In the excited state, chlorine dissociates. This is followed by ground-state reactions.
Molecules 29 04509 g005
Figure 6. Distances of the chlorine/toluol photoreaction. Black: Cl–Cl distance, red and green: Cl-C distances. Upper panel: CPMD, BLYP; Lower panel: Q-Chem, B3LYP. We compare to the B3LYP simulation performed with Q-Chem, because the BLYP simulation with Q-Chem did not converge already in the ground state. In the Q-Chem simulation the components of the system simply fly away from each other. In the CPMD simulation, using periodic boundary conditions, they encounter again and may react on a longer time scale.
Figure 6. Distances of the chlorine/toluol photoreaction. Black: Cl–Cl distance, red and green: Cl-C distances. Upper panel: CPMD, BLYP; Lower panel: Q-Chem, B3LYP. We compare to the B3LYP simulation performed with Q-Chem, because the BLYP simulation with Q-Chem did not converge already in the ground state. In the Q-Chem simulation the components of the system simply fly away from each other. In the CPMD simulation, using periodic boundary conditions, they encounter again and may react on a longer time scale.
Molecules 29 04509 g006
Figure 7. The silver citrate reaction. Only the first two steps are observed in the simulations.
Figure 7. The silver citrate reaction. Only the first two steps are observed in the simulations.
Molecules 29 04509 g007
Figure 8. Distances of the silver citrate photoreaction. Plotted are the three C–C distances that may be broken in order to form CO2. Upper panel: CPMD, BLYP, Lower panel: Q-Chem, B3LYP. Even if a higher temperature was employed in the BLYP calculation, the reactivity is a bit higher for the B3LYP functional. In both cases, the reaction is not immediate.
Figure 8. Distances of the silver citrate photoreaction. Plotted are the three C–C distances that may be broken in order to form CO2. Upper panel: CPMD, BLYP, Lower panel: Q-Chem, B3LYP. Even if a higher temperature was employed in the BLYP calculation, the reactivity is a bit higher for the B3LYP functional. In both cases, the reaction is not immediate.
Molecules 29 04509 g008
Figure 9. The intramolecular hydrogen transfer often connected with CO2 evolution in the silver citrate reaction. Color code: gold: silver; red: oxygen; brown: carbon; white: hydrogen.
Figure 9. The intramolecular hydrogen transfer often connected with CO2 evolution in the silver citrate reaction. Color code: gold: silver; red: oxygen; brown: carbon; white: hydrogen.
Molecules 29 04509 g009
Table 1. Vertical and adiabatic ROKS excitation energies computed with the CPMD and Q-Chem codes in eV.
Table 1. Vertical and adiabatic ROKS excitation energies computed with the CPMD and Q-Chem codes in eV.
CompoundCPMD
BLYP
Q-Chem
BLYP
Q-Chem
B3LYP
TDDFT/B3LYP
[29]
Exp.
[13]
vertical excitation energies
Formaldehyde3.5013.7613.8283.874.07
Acetaldehyde3.8264.0694.155-4.28
Acetone3.9164.1564.2744.344.48
Formamide5.1865.5875.6225.555.65
Acetamide5.1065.5165.6085.57-
Formaldimine4.5764.8734.970-5.0–5.4
Acetaldimine4.8465.1275.240-5.0–5.4
Ethene6.5696.6687.2747.717.11
Propene6.2706.3837.021-6.58
1-Butene5.9416.0506.646-6.61
2-Butene6.1036.1456.804-6.13
Butadiene4.4884.4775.1265.755.93
Hexatriene3.5023.3874.0364.695.15
Propenal3.2443.3723.535--
Pentadienal3.1063.1613.355--
adiabatic excitation energies
Formaldehyde3.1513.3623.352-3.50
Acetaldehyde3.4743.5543.552-3.69
Acetone3.5673.6133.627-4.73
Formamide4.6014.0534.220--
Acetamide4.5823.9514.143--
Formaldimine3.7912.0103.003--
Acetaldimine4.0343.1603.127--
Ethene4.5422.9422.887--
Propene2.9232.9412.901--
1-Butene4.3682.8312.788--
2-Butene4.3922.8982.871--
Butadiene3.9144.0042.508--
Hexatriene3.1482.3082.306--
Propenal2.9232.6593.088--
Pentadienal2.7682.3242.369--
Table 2. Overview of the single simulation runs.
Table 2. Overview of the single simulation runs.
Sim.
Number
AIMD
Method
Equil.
Temp.
[K]
CPU Time
Exc. State
[CPU min per ps]
Cell
Size
3]
Reaction
CPMD code
Chlorine + toluene
1BOMD/BLYP300193,6671000dissociation
2CPMD/BLYP30032671000dissociation
3CPMD/BLYP60035831000dissociation
4CPMD/BLYP90034831000dissociation
Trisilver citrate
5BOMD/BLYP600350,0001728dissociation + proton transfer
6BOMD/BLYP600342,7372744dissociation (2 times)
7BOMD/BLYP600450,7882197dissociation + proton transfer
8CPMD/BLYP60011,9251728-
9CPMD/BLYP60015,5692744dissociation + proton transfer
10CPMD/BLYP60013,2503375-
Trisilver citrate + 8 water molecules
11BOMD/BLYP600424,0162744-
Q-Chem code
Chlorine + toluene
12BOMD/B3LYP30010,131-dissociation
Trisilver citrate
13BOMD/BLYP30038,637-dissociation + proton transfer
14BOMD/B3LYP30052,482-dissociaton (2 times)
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Büchel, R.; Álvarez, L.; Grage, J.; Maniscalco, D.; Frank, I. On the Simulation of Photoreactions Using Restricted Open-Shell Kohn–Sham Theory. Molecules 2024, 29, 4509. https://doi.org/10.3390/molecules29184509

AMA Style

Büchel R, Álvarez L, Grage J, Maniscalco D, Frank I. On the Simulation of Photoreactions Using Restricted Open-Shell Kohn–Sham Theory. Molecules. 2024; 29(18):4509. https://doi.org/10.3390/molecules29184509

Chicago/Turabian Style

Büchel, Ralf, Luis Álvarez, Jan Grage, Dominykas Maniscalco, and Irmgard Frank. 2024. "On the Simulation of Photoreactions Using Restricted Open-Shell Kohn–Sham Theory" Molecules 29, no. 18: 4509. https://doi.org/10.3390/molecules29184509

APA Style

Büchel, R., Álvarez, L., Grage, J., Maniscalco, D., & Frank, I. (2024). On the Simulation of Photoreactions Using Restricted Open-Shell Kohn–Sham Theory. Molecules, 29(18), 4509. https://doi.org/10.3390/molecules29184509

Article Metrics

Back to TopTop