Multistate metadynamics for automatic exploration of conical intersections
Abstract
We introduce multistate metadynamics for automatic exploration of conical intersection seams between adiabatic BornOppenheimer potential energy surfaces in molecular systems. By choosing the energy gap between the electronic states as a collective variable the metadynamics drives the system from an arbitrary groundstate configuration toward the intersection seam. Upon reaching the seam, the multistate electronic Hamiltonian is extended by introducing biasing potentials into the offdiagonal elements, and the molecular dynamics is continued on a modified potential energy surface obtained by diagonalization of the latter. The offdiagonal bias serves to locally open the energy gap and push the system to the next intersection point. In this way, the conical intersection energy landscape can be explored, identifying minimum energy crossing points and the barriers separating them. We illustrate the method on the example of furan, a prototype organic molecule exhibiting rich photophysics. The multistate metadynamics reveals plateaus on the conical intersection energy landscape from which the minimum energy crossing points with characteristic geometries can be extracted. The method can be combined with the broad spectrum of electronic structure methods and represents a generally applicable tool for the exploration of photophysics and photochemistry in complex molecules and materials.
I Introduction
The BornOppenheimer (BO) approximation, the cornerstone of molecular physics, breaks down when molecules get electronically excited. This is most dramatically reflected in the presence of conical intersections (CI) Yarkony (1996); Domcke et al. (2011) between BO potential energy surfaces (PES). A wide range of lightinduced processes including vision in biology BonačićKoutecký et al. (1987); Polli et al. (2010), photostability of DNA Schultz et al. (2004), a plethora of organic photochemical reactions Bernardi et al. (1996), and charge separation in photovoltaic devices Musser et al. (2015) are governed by the efficient nonradiative transitions mediated by a strong coupling between electronic and nuclear dynamics at CIs. In this sense, they play a role in photophysics and photochemistry analogous to the transition states in the groundstate chemistry. However, unlike the transition states, CIs are not isolated points on the PES but form a multidimensional seam. A significant effort has been invested both experimentally and theoretically to identify and characterize CIs Sorgues et al. (2003); Mitrić et al. (2006); Polli et al. (2010); Liebel et al. (2014); Kowalewski et al. (2015, 2016). A number of algorithms has been developed with the aim to search for the minimum energy configurations within the CI seam, since their knowledge allows us to predict the fate of a photoexcited molecule Manaa and Yarkony (1993); Bearpark et al. (1994); Ciminelli et al. (2004). In addition to these local optimization algorithms, methods have been also proposed for location of minimum energy crossing points (MECP), which are closest in terms of a massweighted distance to a predefined reference structure Levine et al. (2008). However, all these methods require a reasonable initial guess for the CI geometry. Alternatively, more elaborate approaches for the automated search of CI geometries based on the anharmonic downward distortion following (ADDF) and artificial forceinduced reaction (AFIR) methods Maeda et al. (2009, 2011); Harabuchi et al. (2013); Maeda et al. (2015), as well as the nudged elastic band method Mori and Martínez (2013), have also been introduced and applied for the exploration of conical intersection seams. Recent studies have emphasized the importance of conical intersections in coherent vibronic dynamics that involves nonadiabatic driving forces associated with the vibrational motion along the tuning modes Duan et al. (2016). While such effects can only be properly investigated in the frame of quantum dynamics simulation, systematic and unbiased sampling of the conical intersection seam might help in construction of reliable (diabatic) Hamiltonians that are needed for such simulations. Since CIs are usually found in the range of several eV above the groundstate equilibrium structure, it would be desirable to develop an accelerated molecular dynamics technique that is able to reach and sample the CI seam starting directly from the groundstate minimum. The metadynamics, introduced by Parrinello and coworkers Laio and Parrinello (2002); Barducci et al. (2008); Raiteri et al. (2006); Bonomi and Parrinello (2010); Tribello et al. (2010); Branduardi et al. (2012), represents an ingeniously simple sampling method that uses collective variables (CVs) to drive transitions between different barrierseparated basins on the PES, allowing for systematic sampling of the PES as well as determination of free energies. Depending on the problem to be solved, a wide variety of CVs has been used, ranging from simple geometrical quantities to complex variables constructed by machine learning and dimensionality reduction techniques Tribello et al. (2010); Ceriotti et al. (2011). In this article, we introduce a multistate metadynamics method for automatic exploration of conical intersection seams and identification of MECP.
Ii Method
The central idea of our multistate metadynamics is to use the energy gap between the ground and excited electronic states as a collective variable to drive the system toward the CI starting from any groundstate structure. For this purpose, the system is propagated using the Newtonian equations of motion, which for the th particle read
(1) 
where represents the groundstate PES and is a historydependent bias potential:
(2)  
which is dependent on the modified energy gap , and and represent the fixed height and width, respectively. The bias potential is updated at regular time steps but only if the value of the gap is larger than a numerical threshold , which is enforced by the presence of the Heaviside theta function in Eq. (2). Starting from a minimum on the groundstate PES, the bias potential will drive the system toward the CI seam by systematically lowering the energy gap [see Fig. 1(a) for illustration]. After the intersection point is reached for the first time, the molecular electronic Hamiltonian, which is initially diagonal in the BO approximation, is extended by introducing a further biasing potential into the offdiagonal elements according to
(3) 
Subsequently, the metadynamics is continued on a locally modified potential energy surface obtained by diagonalization of the above Hamiltonian which gives rise to a modified PES with the effective energy gap
(4) 
The offdiagonal bias serves to locally open the gap between the eigenstates of the modified Hamiltonian. The metadynamics then drives the system to the next intersection point. In order to prevent the return to the previously sampled regions of the CI seam, the is made dependent on a collective variable , which should be capable of distinguishing different molecular configurations and is updated only if the energy gap is below according to
(5)  
One possible choice of is to use the distance matrix of a molecule, which is unique but inconvenient because of its high dimensionality and symmetry ambiguities. A better choice is to use its scalar invariants such as its lowest eigenvalue, or various topological indices Mihalić et al. (1992). As a CV for the offdiagonal bias in our simulations, we choose the 3DWiener number Bogdanov et al. (1989) defined as
(6) 
where is the number of atoms and are interatomic nonhydrogen distances. Its correlation with molecular shape made it a reliable topographical descriptor in numerous studies where structural unambiguity is important. Since the conical intersection seam is of dimension , where is the number of degrees of freedom, and is a scalar function of the distances, there is no restriction to any specific region on the seam. Therefore, the most possible structural variability is achieved and can be further extended by running several trajectories with different initial conditions.
We wish to point out that as long as no Gaussians have been added to the offdiagonal bias , exactly equals the BO gap . This can be derived from Eq. (4) and is due to the construction of that reduces to in case of the nonexisting offdiagonal coupling. Since is updated only in the vicinity of the CI seam, this holds true for the complete pathway until the intersection is reached for the first time [see Fig. 1(a)]. At this point, the metadynamics potential is already overcompensating the minimum on the groundstate PES where the simulation had been started, resulting in a persistent force toward structures with small . If now a Gaussian bias is added to at the current value of the collective variable , will increase, thus opening again the gap between the two states [see Fig. 1(b)]. However, the bias potential disfavors nonzero values of the energy gap, which in turn forces the molecule to change the value of the collective variable such that the intersection seam is reached again [see Fig. 1(c)]. In this way, the whole CI seam can be “unzipped”, allowing for automatic exploration of its energy landscape. In order to calculate the force in Eq. (1), the gradient of the modified energy gap is needed:
(7) 
which requires the calculation of the gradients of the BO energies that can be provided by any suitable electronic structure method, as well as additional differentiation of :
(8)  
Like the potential itself, it is given by a sum of Gaussians, but each of them is multiplied by an additional factor dependent on the gradient of the collective variable. For the 3DWiener number, is defined by
(9) 
The algorithm can be summarized with the following steps:

Initialize metadynamics trajectory with effective energy gap as CV and choose an appropriate CV for the offdiagonal bias potential.

Integrate the Newtonian equations of motion using a conventional molecular dynamics algorithm. If the actual time step is a multiple of , proceed with one of the following steps:

If , add new Gaussian to the diagonal metadynamics potential .

If , add new Gaussian to the offdiagonal coupling potential and diagonalize the electronic Hamiltonian to obtain modified PES.


Stop the simulation after a defined number of steps.

Identify the plateau regions in the CI energy landscape and select structures for MECP optimization.
The use of a threshold in Eqs. (2) and (5) for the switch between on and offdiagonal metadynamics enables the addition of Gaussians to only if the trajectory is close to the CI seam. At the same time, it prevents biasing small values of that would lead back to larger energy gaps again. A reasonable choice for is to use the same energy that is employed as Gaussian width for the . The Gaussian width and height for the offdiagonal bias should be chosen sufficiently small in order to bias the previously visited configurations as local as possible. Thus, it can be ensured that the system stays in close proximity to the intersection seam most of the simulation time. Unlike conventional metadynamics, the aim of this algorithm is not to converge the metadynamics potential to the free energy surface and the results are therefore less sensitive to the choice of parameters.
Periods of the trajectory with low represent structures close to the intersection seam. It is convenient to apply a local CI optimization method on snapshots from these regions in order to group the resulting MECPs by energy.
Iii Results and Discussion
In order to illustrate the ability of the multistate metadynamics to sample conical intersections, we choose the furan molecule, which exhibits versatile photochemistry that is still not completely understood and has been a subject of a number of recent experimental and theoretical studies Stenrup and Larson (2011); Gavrilov et al. (2008); Fuji et al. (2010); Spesyvtsev et al. (2015); Oesterling et al. (2017). Here, the electronic structure of furan is described in the frame of the ab initio complete active space selfconsistent field (CASSCF) method as implemented in the MOLPRO quantum chemical package Werner et al. (2011) with an active space of 10 electrons distributed over 9 orbitals together with the 631G* atomic basis set, as in previous quantumchemical studies Oesterling et al. (2017). The Newtonian equations of motion in the metadynamics simulations were integrated using the velocity Verlet algorithm Verlet (1967) with a step size of and the temperature was kept constant at utilizing the Berendsen thermostat Berendsen et al. (1984). Multistate metadynamics simulations were run using initial conditions sampled from a 2pslong MD trajectory propagated in the electronic ground state. Gaussian bias potentials with a width of and a height of were added every 100 steps () to the diagonal bias. For the offdiagonal coupling , we apply the 3DWiener number (excluding hydrogen atoms) as CV with Gaussians of width and height. The update threshold was set to . As can be seen from Fig. 2(a), the energy gap decreases during the multistate metadynamics simulation from an initial FranckCondon value of to zero within the first 1.4 ps. At this point, the CI is reached for the first time and the offdiagonal coupling shown in Fig. 2(b) is automatically switched on. Together with the previously builtup bias potential , this forces the system to continue exploring the CI energy landscape. The time evolution of the ground and excitedstate PES for three runs with different initial conditions sampled from the ground state trajectory is shown in Fig. 3. In each case, the system is initially driven from a local minimum toward the conical intersection, leading to the closure of the energy gap. Subsequently, the dynamics explores the CI landscape and the energy gap remains zero within the predefined numerical threshold. As can be seen from Fig. 3, CI energy landscape exhibits several plateaus and depending on the CI character, some of them are separated by relatively large barriers. Low barriers are in general present for transitions between various ringpuckering CIs, as is seen in the trajectory depicted in Fig. 3(a). On the contrary, the trajectory shown in Fig. 3(b) leads to the CI structures exhibiting ring opening as well as fragmentation. Molecular fragmentations are typically found after longer simulation times since they are characterized by large Wiener numbers. It is interesting to note that both molecular geometries that have previously been held responsible for the ultrafast deactivation of furan Fuji et al. (2010); Oesterling et al. (2017), a ring puckering of the oxygen atom and the lowenergetic ring opening, have been identified in successive order in the third trajectory shown in Fig. 3(c). The found geometries have been subsequently fully optimized using the BearparkRobb local optimization algorithm Bearpark et al. (1994) and their branching plane vectors have been determined (cf. Fig. 3).
Iv Conclusion
In summary, we have developed a multistate extension of the metadynamics with the aim to automatically explore conical intersection seams between BO PES. Biasing metadynamics potentials are introduced as offdiagonal elements into the multistate electronic Hamiltonian and MD simulations are run on modified potential energy surface obtained by diagonalization of the latter, using the energy gap as a collective variable. The method can be easily implemented in the frame of any electronic structure method capable of providing energy gradients and excitation energies. It can be applied to explore photochemical reaction pathways, nonradiative relaxation channels, and photophysics of complex molecular systems. As an illustration, we have performed simulations starting from the groundstate structure of the furan and have demonstrated that the conical intersection landscape can be efficiently mapped, allowing us to systematically identify a large number of minimum energy crossing points that can mediate nonradiative relaxation and photochemical reactivity.
Acknowledgement
Funding by the European Research Council (ERC) Consolidator Grant DYNAMO (Grant No. 646737) is gratefully acknowledged.
References
 Yarkony (1996) D. R. Yarkony, Rev. Mod. Phys. 68, 985 (1996).
 Domcke et al. (2011) W. Domcke, D. R. Yarkony, and H. Köppel, eds., Conical Intersections I and II (World Scientific, Singapore, 2004, 2011).
 BonačićKoutecký et al. (1987) V. BonačićKoutecký, J. Koutecký, and J. Michl, Angew. Chem. Int. Ed. 26, 170 (1987).
 Polli et al. (2010) D. Polli, P. Altoè, O. Weingart, K. M. Spillane, C. Manzoni, D. Brida, G. Tomasello, G. Orlandi, P. Kukura, R. A. Mathies, M. Garavelli, and G. Cerullo, Nature 467, 440 (2010).
 Schultz et al. (2004) T. Schultz, E. Samoylova, W. Radloff, I. V. Hertel, A. L. Sobolewski, and W. Domcke, Science 306, 1765 (2004).
 Bernardi et al. (1996) F. Bernardi, M. Olivucci, and M. A. Robb, Chem. Soc. Rev. 25, 321 (1996).
 Musser et al. (2015) A. J. Musser, M. Liebel, C. Schnedermann, T. Wende, T. B. Kehoe, A. Rao, and P. Kukura, Nat. Phys. 11, 352 (2015).
 Sorgues et al. (2003) S. Sorgues, J. M. Mestdagh, J. P. Visticot, and B. Soep, Phys. Rev. Lett. 91, 103001 (2003).
 Mitrić et al. (2006) R. Mitrić, V. BonačićKoutecký, J. Pittner, and H. Lischka, J. Chem. Phys. 125, 024303 (2006).
 Liebel et al. (2014) M. Liebel, C. Schnedermann, and P. Kukura, Phys. Rev. Lett. 112, 198302 (2014).
 Kowalewski et al. (2015) M. Kowalewski, K. Bennett, K. E. Dorfman, and S. Mukamel, Phys. Rev. Lett. 115, 193003 (2015).
 Kowalewski et al. (2016) M. Kowalewski, K. Bennett, J. R. Rouxel, and S. Mukamel, Phys. Rev. Lett. 117, 043201 (2016).
 Manaa and Yarkony (1993) M. R. Manaa and D. R. Yarkony, J. Chem. Phys. 99, 5251 (1993).
 Bearpark et al. (1994) M. J. Bearpark, M. A. Robb, and H. B. Schlegel, Chem. Phys. Lett. 223, 269 (1994).
 Ciminelli et al. (2004) C. Ciminelli, G. Granucci, and M. Persico, Chem. Eur. J. 10, 2327 (2004).
 Levine et al. (2008) B. G. Levine, J. D. Coe, and T. J. Martínez, J. Phys. Chem. B 112, 405 (2008).
 Maeda et al. (2009) S. Maeda, K. Ohno, and K. Morokuma, J. Phys. Chem. A 113, 1704 (2009).
 Maeda et al. (2011) S. Maeda, R. Saito, and K. Morokuma, J. Phys. Chem. Lett. 2, 852 (2011).
 Harabuchi et al. (2013) Y. Harabuchi, S. Maeda, T. Taketsugu, N. Minezawa, and K. Morokuma, J. Chem. Theory Comput. 9, 4116 (2013).
 Maeda et al. (2015) S. Maeda, T. Taketsugu, K. Ohno, and K. Morokuma, J. Am. Chem. Soc. 137, 3433 (2015).
 Mori and Martínez (2013) T. Mori and T. J. Martínez, J. Chem. Theory Comput. 9, 1155 (2013).
 Duan et al. (2016) H.G. Duan, R. J. D. Miller, and M. Thorwart, J. Phys. Chem. Lett. 7, 3491 (2016).
 Laio and Parrinello (2002) A. Laio and M. Parrinello, Proc. Nat. Acad. Sci. USA 99, 12562 (2002).
 Barducci et al. (2008) A. Barducci, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 100, 020603 (2008).
 Raiteri et al. (2006) P. Raiteri, A. Laio, F. L. Gervasio, C. Micheletti, and M. Parrinello, J. Phys. Chem. B 110, 3533 (2006).
 Bonomi and Parrinello (2010) M. Bonomi and M. Parrinello, Phys. Rev. Lett. 104, 190601 (2010).
 Tribello et al. (2010) G. A. Tribello, M. Ceriotti, and M. Parrinello, Proc. Nat. Acad. Sci. USA 107, 17509 (2010).
 Branduardi et al. (2012) D. Branduardi, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 8, 2247 (2012).
 Ceriotti et al. (2011) M. Ceriotti, G. A. Tribello, and M. Parrinello, Proc. Nat. Acad. Sci. USA 108, 13023 (2011).
 Mihalić et al. (1992) Z. Mihalić, S. Nikolić, and N. Trinajstić, J. Chem. Inf. Model. 32, 28 (1992).
 Bogdanov et al. (1989) B. Bogdanov, S. Nikolić, and N. Trinajstić, J. Math. Chem. 3, 299 (1989).
 Stenrup and Larson (2011) M. Stenrup and Å. Larson, Chem. Phys. 379, 6 (2011).
 Gavrilov et al. (2008) N. Gavrilov, S. Salzmann, and C. M. Marian, Chem. Phys. 349, 269 (2008).
 Fuji et al. (2010) T. Fuji, Y.I. Suzuki, T. Horio, T. Suzuki, R. Mitrić, U. Werner, and V. BonačićKoutecký, J. Chem. Phys. 133, 234303 (2010).
 Spesyvtsev et al. (2015) R. Spesyvtsev, T. Horio, Y.I. Suzuki, and T. Suzuki, J. Chem. Phys. 143, 014302 (2015).
 Oesterling et al. (2017) S. Oesterling, O. Schalk, T. Geng, R. D. Thomas, T. Hansson, and R. de VivieRiedle, Phys. Chem. Chem. Phys. 19, 2025 (2017).
 Werner et al. (2011) H.J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2, 242 (2011).
 Verlet (1967) L. Verlet, Phys. Rev. 159, 98 (1967).
 Berendsen et al. (1984) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984).