The Use Of Non-Unitary Initial Sampling In The Multi-Configurational Ehrenfest Method
Introduction
Solving the Schrödinger equation for a Hamiltonian with many degrees of freedom, though academically understood is, in practise, impossible. The computation cost of modelling the problem grows exponentially with each extra degree of freedom. Thus, it has been necessary to develop other methods to model quantum systems using semiclassical behaviour or ideas to lessen the computational cost. Many of these models can be traced to the work of Miller and co-workers, in the early 1970s, which introduced the Initial Value Representation [IVR] (1, 2). IVR was initially presented as a more general expression of the S matrix in terms of classical quantities which could be used to treat non-classical transitions more accurately (1).
However, the method was not practical: requiring a root search for the classical trajectories which proved too computationally expensive, even for simple chemical reactions (3). Consequently, Heller proposed a system using ‘frozen gaussian’ [FG] wave packets (4-6). The FG method works by breaking the initial wave function into a superposition of Gaussian wave packets, called coherent states. These are then propagated along classical trajectories, the wave packet is then reunited at various time intervals (5). But FG method ‘lacked a rigorous semiclassical basis, in the sense of passing over to the [van Vleck] limit of the quantum propagator’ (3). But in 1986 Herman and Klux derived the FG approximation using the van Vleck propagator (7-10). This revealed a new pre-factor that had been missing from Heller’s method. Later work showed that Herman and Klux’s [HK] method could be written in terms of the IVR thus linking the Frozen Gaussian approximation back to the Initial value Representation method (11, 12). Despite the clear success and modern use of HK methods (13-16) they are still semiclassical and so do not fully describe the quantum wave function. Hence, there was a clear necessity to build a model to accurately describe the quantum mechanics of a system one of which was the Coupled Coherent State [CCS] method. The Coupled Coherent States method was first proposed in 2000 by Shalashilin and Child (17).
The quantum Schrödinger equation is approximated as a finite sum thus giving rise to the propagation of the wave function in a continuous basis of coherent states in phase space. This is notably different as previous models used coordinate space. Hence, CCS permitted the calculation of quantum amplitudes, unlike in the semiclassical approaches, to be, in theory, exact (18). Insert some stuff about it being good!The CCS method was then generalised to a allowing multiple potential energy surfaces (19). The Multi-Configurational Ehrenfest [MCE] method was first developed in 2009 by Shalashilin (19) and a revised version (20) was published in 2010. The first version [MCEv1] was shown to have good agreement with the MCTDH benchmark result for the spin-boson model (19). The revised method [MCEv2] also showed good agreement with the MCTDH benchmark calculations for the simulation 24-dimensional model of pyrazine (20).
Further MCEv2 has been used to develop a method for ab initio nonadiabatic molecular dynamics (21). It is with this method that it has been possible to run simulations of real molecules undergoing photoinduced dynamics (22-25). However, a lack of agreement between MCEv1 and MCEv2 was discovered (26) while looking at model systems. Thus, there clearly is a need to improve the sampling methods employed in MCEv2. This report will give an outline of the theory behind both versions of the MCE method and the mathematical equations that underpin them. Then the lack of coherence between both versions of MCE when looking at model systems (26) will be briefly discussed. Subsequently a possible solution, the aim of this project, will be proposed as a way of improving the sampling method to improve the overall agreement of the simulations with benchmark calculations.
The Multi-Configurational Ehrenfest method
The self-consistent mean-field formalism originally developed by Billing (27) is the starting point for the MCE approach. The Ehrenfest approximation which assumes that a system can be split into quantum and classical parts. Thus, its Hamiltonian can be written as a sum of the quantum system, classical bath and the interactions between the two, where are the system’s quantum operators and are the coordinates and momenta classical vectors for the bath. Hence the trajectory of the classical bath is given by the Hamilton’s equations, resulting in the Hamiltonian. This is averaged over the wave function. The evolution of the time dependent coefficients of the Ehrenfest equations, giving the time evolution of the quantum system are found by,MCE uses the Ehrenfest dynamics equations as starting off point. However, rather than treating the bath classically the MCE method treats all degrees of freedom [DOF] at a fully quantum level.
All the following equations are given in atomic units, i. e. The equations are given here for the interaction of two potential energy surfaces [PESs], but they could be generalised to more states or to one giving the equations for CCS.
Equations for MCEv1
The ansatz of the wave function for MCEv1 is given by, The quantum system states consist of the orthonormal number states. The kets and are the electronic states and describes the nuclear wave function. The bath is a product of M one-dimensional CS of individual bath states,. Individual configurations are coupled via the amplitudes ; the time dependence of these coefficients and the trajectory of give the evolution of the wave function. It is show in Ref. (19) that the time dependence of is found by applying the variational principle to the quantum Lagrangian. This has the Hamiltonian, whose off-diagonal elements give the coupling between the two PES. The index ‘ord’ shows that the creation operators precede the annihilation operators, this is normal ordering. The Lagrangian is expanded and the variational principle, applied which gives the equations of motion for and. Variation of eventually yields where the Ehrenfest Hamiltonian for a system with two electronic states has the form. A full derivation of these results can be found in Ref. (19). The oscillating exponent of the classical action is where, is the classical coherent state action. Using Eq. and variation of the amplitudes it is possible to write and as a product of a smooth preexponential factor. This then makes the coupled equations.
Equations for MCEv2
The second form of MCE was developed because in the original method the coefficients of the system basis states were coupled, concurrently within the same configuration (intra-configurational) and between configurations (inter-configurational).
The consequential effect being the coherent states in the bath became coupled. To alleviate this problem the coupling of the configurations is separated from the coupling of the PES by the presence of an extra amplitude This gives a new ansatz of. Propagation of the wave function is now completed by a concurrent process propagating and. The evolution of are found using equations and the coefficients of can be found using the Schrödinger equation. In Ref. (20) it is shown how Eq. can be substituted into the Schrödinger equation to give a set of linear equations, where, is the overlap matrix of the Ehrenfest basis wave functions. The elements of the matrix are. Further, in Eq. and within the matrix is. [20]Differences in result between versions of MCE As it has been shown, mathematically, both versions of the MCE method have a very similar. MCEv1 is capable of simulating model systems (19) and MCEv2 has been successfully applied to real molecules (21-25).
However, when given the same model system there is not the agreement that would be expected (26). The Spin Boson model was used to make a comparison between the two formulations of the MCE method. The difference between the results is immediately clear: MCEv2 has overly emphasised oscillations for both symmetric and asymmetric wells whereas MCEv1 agrees, almost exactly, with the benchmark MCTDH result. Further, increasing ‘the size of the basis set or the number of repetitions did not improve the agreement for the MCEv2 results’(26). Comparisons of the two versions of the MCE method using the Spin Boson model with symmetric [left] and asymmetric [right] wells. All results are compared to the numerically exact MCTDH results. Copied from Ref. (26)
Therefore, it was concluded there was an issue with the way the basis was set up. Either there is a lack of coupling; the basis set is acting as an ensemble of independent noninteracting basis functions or the basis is overly coupled; the basis is being guided by a set of overly similar trajectories resulting in it acting as a larger single basis function. It was shown to be the later by looking at the overlap between the coherent states which was much higher for a longer period of time in MCEv2 than MCEv1. Consequently, the basis set did not spread out as much in phase space. This is shown in Figure 2, both formulations had the same starting point but while MCEv1 spreads with time MCEv2 maintains a high density at the centre. Thus, the coherent states cannot describe a large enough area of phase space sufficiently to then simulate the full quantum nature of the system (26).
A comparison of how the density of basis functions in phase space change over time for both versions of the MCE method. The Spin Boson model with an asymmetric well has been used. Copied from Ref. (26)
A Possible Solution: Non-Unitary Initial Sampling
Despite the progresses trains and cloning have made to the MCE method there is still room for improvement (26). This could be provided by the use of non-unitary sampling. The Spin-boson model will be used as a simple example as it has only two electronic states, this could be applied to systems with a greater number of states. Currently, the initial Ehrenfest configuration places an amplitude of 1 for the populated electronic state and an amplitude of 0 for the unoccupied state,, which gives the initial propagating wave function. So, for all configurations the Ehrenfest amplitudes are the same, 1 and 0. Only the values of are selected randomly. Thus, if the distribution of is dense all the trajectories will be very similar. The possible solution is to not only randomly generate but also the amplitudes, and. This results in the following equation for the initial propagating wave function. The new amplitude functions take the form, where and are two random numbers between and. The values of will be generated as they normally as described in Refs. (19, 20). The values of are generated in a similar way to described in Ref. For simplicity the matrix that represents the overlap of the Ehrenfest configurations Eq. will be expressed as a single entity. The identity operator is expressed via the elements of the inverse overlap matrix.
It can now be seen that the values of can be found by using Eq. where. As described above and in Ref. the inverse is not found but instead a set of linear equations are solved to find the vector of amplitudes. Though perhaps counterintuitive, for a quantum system, this new sampling makes logical sense. By setting the amplitudes for the Ehrenfest states to be 1 and 0, for the populated and unpopulated states respectively, when finding the amplitudes of models the classical behaviour of the system. But by allowing the Ehrenfest amplitudes to vary for the populated and unpopulated states the model is now allowing for some classically disallowed trajectories. In a quantum system, even initially, there is the possibility of trajectories existing in the unoccupied state. Find paper to back this up!
