Trajectory integration with potential energy discontinuities

https://doi.org/10.1016/j.jcp.2009.11.025Get rights and content

Abstract

Many approximate methods of quantum chemistry yield potential energy surfaces with discontinuities. While clearly unphysical, such features often fall within the typical error bounds of the method, and cannot be easily eliminated. The integration of nuclear trajectories when the potential energy is locally discontinuous is obviously problematic. We propose a method to smooth out the discontinuities that are detected along a trajectory, based on the definition of a continuous function that fits locally the computed potential, and is used to integrate the trajectory across the discontinuity. With this correction, the energy conservation error can be reduced by about one order of magnitude, and a considerable improvement is obtained in the energy distribution among the internal coordinates.

Introduction

Classical trajectories for the nuclear motion are the basis of molecular dynamics treatments of ground and excited state processes. In this approach the potential energy for the nuclear motion is one of the eigenvalues of the electronic hamiltonian, that depends on the nuclear coordinates as parameters. When a direct strategy is used, the potential energy surface (PES) and its gradient are computed at every time step of the trajectory, i.e. for a given molecular geometry by means of quantum chemistry methods [1], [2], [3], [4]. The alternative is to use analytic potential energy functions, either of standard molecular mechanics (MM) type or devised ad hoc for specific systems and processes. The construction of accurate PESs for reaction dynamics and/or multistate processes can be a complex problem, especially in the case of state degeneracies (conical intersections and crossing seams [5]), and its computational cost increases exponentially with the number of internal coordinates. On the other hand, the cost of direct methods depends linearly on the number and duration of the trajectories, which makes them convenient for the simulation of processes in the picosecond time scale.

One of the drawbacks of the direct use of quantum chemistry methods in trajectory simulations is the possible occurrence of discontinuities in the computed PES. Most methods for the solution of the many-electron problem rely on the variational optimization of non-linear parameters, and may therefore admit several solutions. In other words, the energy, as a function of the electron density or wavefunction parameters, may possess one or more local minima in addition to the global one. If the molecular geometry is varied by small steps, as in trajectory calculations or when scanning a potential energy curve, it is not uncommon to switch from one solution to another one, with a “sudden” change in the electronic energy and in all other properties. Even if the change takes place gradually, but in a very small interval of internal coordinates, it is seen as a discontinuity if one samples the PES by finite time or space steps.

The SCF, DFT and MCSCF methods provide good examples of such behaviour [6], [7], [8], [9], [10], [11], [12], [4]. In CASSCF, a variant of the general MCSCF approach, the definition of the active space (number of active electrons and orbitals) suffices to determine the global minimum. However, swapping an active orbital whose occupation number is close to 0 or 2 with a virtual or doubly occupied orbital, may yield another stable solution of slightly higher energy [10], [11], [12]. The particular solution found at a given geometry by the iterative self-consistent procedure depends on many technical details, such as the starting molecular orbitals (MO) and the algorithms that accelerate and/or stabilize the convergence. For instance, the practice of using the MOs of the previous step in a trajectory calculation as the starting point for the SCF or MCSCF procedure may help to maintain a particular solution. However, a solution corresponding to the global minimum for a certain range of molecular geometries, may be just a local minimum at other geometries, and disappear altogether (no minimum) at others yet; in these cases, a sudden switch between two solutions is unavoidable. Configuration Interaction (CI) calculations, of variational or perturbative type, can be based on the SCF or MCSCF MOs in order to improve the treatment of electron correlation. The CI results depend on the MOs, so again they may undergo sudden changes due to switching between different SCF or MCSCF solutions, even if the switch implies a very small change in the SCF or MCSCF energies (a typical case is the conversion from localized to delocalized orbitals). Similar problems are also met with the floating occupation (FO) semiempirical SCF-CI method set up by Granucci et al. [2]. The FO-SCF-CI is an inexpensive alternative to state average (SA) MCSCF for the balanced determination of two or more electronic states, and can be also applied in the ab initio framework [13]. In SA-MCSCF calculations, often used in ab initio multistate molecular dynamics, a further source of discontinuities in the solution is the occurrence of root-switching, i.e. a change in the composition of the subspace of electronic states that are optimized [4]. Notice that root-switching does not affect the FO-SCF-CI results. We end this short overview of possible sources of PES discontinuities, by noting that minor problems are also met in pure MD or in mixed QM/MM treatments, because of cutoffs in the evaluation of small interaction terms or integrals.

The presence of discontinuities in the PES is an unphysical feature, but in many cases cannot be easily eliminated. Two solutions of the electronic problem can be qualitatively and even quantitatively satisfactory, each one in its own range of molecular geometries, and yet not match at intermediate geometries. If the mismatch in the energy of the relevant electronic state does not exceed the typical error of the computational method, it is reasonable to smooth out the discontinuity in order to be able to integrate the trajectories across it. When the PES are prepared beforehand by fitting or interpolation, the discontinuities are automatically smoothed away, and possibly not even noticed. When using the direct approach, one cannot anticipate at which molecular geometry a discontinuity will be met. Even worse, if the convergence of the SCF or MCSCF algorithms is facilitated by using the MOs of the previous trajectory step as a starting point, the switch to another solution is determined by the trajectory itself, the time step, and other technical details. Therefore, we propose a fitting procedure, with a modification of the trajectory integration, to be applied locally when a discontinuity is detected. The mathematical formulation is described in the next section, and the results of tests with a model potential are shown in Section 3.

Section snippets

Method

We shall describe a trajectory by nvar cartesian coordinates X(t), functions of time, and the related velocities V(t). The kinetic and potential energies are also functions of time: T(t)T(V(t)) and U(t)U(X(t)). The PES U(X) comes in two sheets, with a discontinuity at the border between them, that we shall call the discontinuity seam. The seam is a set of points of dimension nvar-1, whose location in the coordinate space is a priori unknown. The total energy E=T+U should be constant, within

Tests with a model potential

We ran several tests with the model potentialU(X)=Ae-BX1+VDH(X1-XD)+r=2nvarCr[Xr-HDH(X1-XD)]2where H(X) is the Heaviside function. The discontinuity is found at X1=XD, and has a “vertical” and a “horizontal” component, VD and HD, respectively. Most of the tests were done with the following values of the parameters: A=0.08,B=0.3,Cr=0.20, XD=6, all in atomic units. The mass was M=20,000a.u. (same for all the coordinates). With this form of the discontinuity the exact solution is easily computed,

Concluding remarks

We have proposed a simple procedure to deal with small discontinuities in the potential energy, in trajectory simulations of the molecular dynamics. Such discontinuities may occur when solving the electronic problem “on the fly”, with quantum chemistry methods, during the integration of the nuclear trajectory. They are often under the “chemical” accuracy, i.e. one or a few kcal/mol, depending on the problem.

Our procedure consists of smoothing out a discontinuity by replacing the potential

Acknowledgments

We are grateful to Professor Yu Zhuang, of Texas Tech University, for very helpful discussions. Most of this work has been carried out in Pisa, during the “Summer Abroad Program 2009” of the NSF-PIRE project “Simulation of Electronic Non-Adiabatic Dynamics for Reactions with Macromolecules, Liquids, and Solids”, of which one of us (P.H.) was one of the selected participants.

References (14)

  • M. Barbatti et al.

    J. Photochem. Photobiol. A

    (2007)
  • A.K. Mazur

    J. Comput. Phys.

    (1997)
  • T. Vreven et al.

    J. Am. Chem. Soc.

    (1997)
  • G. Granucci et al.

    J. Chem. Phys.

    (2001)
  • J.M. Hostettler et al.

    J. Chem. Phys.

    (2009)
  • M. Persico

    Electronic diabatic states: definition, computation and applications

  • J.S. Sears et al.

    J. Chem. Phys.

    (2006)
There are more references available in the full text version of this article.
View full text