Skip to main content
European Commission logo print header

On-the-fly nonadiabatic quantum dynamics suitable for large biomolecules:<br/>Developing the DD-vMCG method

Final Report Summary - DQDPROT (On-the-fly nonadiabatic quantum dynamics suitable for large biomolecules:<br/>Developing the DD-vMCG method)

The rapid development of molecular spectroscopy techniques have made it possible to study femto- and, recently, atto-second events, lifting the veil of the nuclear and electronic dynamics, respectively. Accordingly, the amount of experimental data on molecular and chemical reaction mechanisms is constantly increasing. In particular, experiments on non-adiabatic effects, governing dynamics at the intersection of electronic states of the system upon interaction with light are growing in importance and demand accurate theoretical methods able to account for quantum effects. Although it is now relatively straightforward to solve numerically the time-dependent Schroedinger equation for a small system consisting of a few nuclei and corresponding electrons, such calculations become quickly unfeasible with increasing system dimensionality. At the same time, classical and semi-classical methods are usually incapable of capturing explicit quantum effects such as proton tunnelling and dynamics in the regions where nuclear and electronic motion becomes strongly coupled, e.g. at a conical intersection between potential surfaces.

There are plethora of large systems of biological and medicinal importance, in which such effects play an important role, with the green fluorescent protein (GFP) being a prominent example. An invaluable tool in cell biology in the last twenty years, it was the first discovered member of the family of fluorescent proteins used in bioimaging techniques and promising other applications such as photodynamic cancer treatment and optical highlighting. GFP’s mechanism of action involves initial absorption of blue light, promoting an excited state proton transfer (ESPT) via the pre-organised hydrogen-bond network, leading to an anionic chromophore emitting green light. This is the only known ESPT taking part in biologically active molecules, and although thoroughly studied experimentally, its mechanism is still a matter of debate. To study this process theoretically, a quantum dynamical treatment of the mechanism is required as ESPT may depend on both proton tunnelling and non-adiabatic effects. At the same time, even a reduced-dimensional model that may in principle capture the physics of the process, will consist of a least about 50 atoms, making it impossible to be treated with conventional quantum dynamical methods.

A method of the on-the-fly quantum dynamics, called direct dynamics variational multi-configuration Gaussian (DD-vMCG), has been developed. It overcomes some of the crucial difficulties of the conventional quantum nuclear dynamics and is in principle applicable to a molecular system of arbitrary size. In this method, localized Gaussian basis functions are used to express the nuclear wavepacket, and direct dynamics becomes possible by means of the local harmonic approximation of the potential energy surface around the Gaussian functions centres. It opened a unique opportunity to study sensitive quantum effects in larger molecules in a nearly black-box fashion with minimal approximations adopted in the theory of the method. The main aim of the current project was to make the direct quantum dynamical calculations more efficient and possible for large non-adiabatic systems of biological importance via further development of the DD-vMCG program, and applying it to investigate properties of GFP and, possibly, PYP (photoactive yellow protein).

The project as originally envisioned was split into three main objectives: (a) update the DD-vMCG code and develop a new on-the-fly diabatization method for efficient quantum dynamics of large non-adiabatic molecular systems, (b) develop a new quantum mechanical/molecular mechanical (QM/MM) interface to the DD-vMCG program and (c) perform full-dimensional calculations of GFP and PYP photoactivation.

The early stages of the work were planned to be devoted to making light improvements over the existing code and showing convergence of the direct dynamics to the conventional quantum dynamics methods using a ground state proton transfer reaction occurring in the medium-sized salicylaldimine molecule. Conventionally, direct dynamics is performed with the use of frozen-width Gaussian basis functions due to a high stability of such approach. During this study, however, we came to an understanding of the importance of choosing the right Gaussian basis function width parameter, or even of using the non-frozen-width basis functions, for highly anharmonic modes. The work in this direction was started, but due to a high instability of the dynamics with the use of variable-width Gaussians, was postponed, and carefully adjusted width parameters were used to study salicylaldimine proton-transfer reaction. Recently, a paper on these calculations was accepted in the Journal of Chemical Physics.

Additionally, other developments of the DD-vMCG program were found to be necessary for a more efficient and stable performance of the code and a lot more effort was spent on its upgrades than expected originally. Much work was required on the integration scheme used for propagation of the basis functions due to numerical instabilities that became apparent in the studied system when approaching convergence. Specifically, a new method was implemented to regularize the so-called С-matrix, inversion of which is required at every integration step, and the size of which grows proportionally to the size of the system times the number of basis functions used in propagation. Regularisation is necessary due to the possibility of the C-matrix becoming singular, and quickly becomes a computational bottleneck with the growth of the system size. The new method showed good performance and allowed for usage of matrix decomposition techniques, leading to a significant improvement of the dynamics scaling. The handling of the database of energies, gradients and Hessians was also significantly improved and its performance optimized. Also, we implemented a Hessian update procedure, allowing for a further speed-up of direct dynamics, and being especially important for expensive excited-state calculations. These developments, together with an in-depth review of the current state of the DD-vMCG program were published in the paper in the International Reviews in Physical Chemistry.

In parallel, we started ab initio calculations of the GFP ESPT reaction. We decided to stick initially to the cluster model of the active site, containing 48 atoms. We found stationary points along the reaction path at both the ground and first excited states, and performed benchmark calculations with the different quantum chemical methods. As a reasonable trade-off between accuracy and computational effort, we have chosen TD-DFT (wB97X-D functional) method, taking into account necessity to perform hundreds or thousands of energy, gradient and optionally Hessian calculations during the wavepacket propagation. Initially, we aimed at fitting a reduced-dimensional reaction potential energy surface in terms of the small subset of normal modes obtained at the ground state transition state geometry, in order to perform conventional quantum dynamics benchmark calculations, but it was observed that in order to be able to even qualitatively describe the minimum-energy reaction path profile and relative energies of reaction minima, one needs at least 30 normal modes. Therefore, necessity of performing direct quantum dynamical simulations became obvious.

With much development already performed, DD-vMCG method was still not scaling well enough to efficiently propagate 30- and higher-dimensional wavepackets due to the rapid growth of the C-matrix. Therefore, we decided to implement the possibility to express the wavepacket as a product of reduced-dimensional Gaussians as opposed to the full-dimensional Gaussians used in the vMCG. That would make the C-matrix block-diagonal and remove the computational bottleneck related to its inversion. In order to do so, significant work must have been done at the code, and we have only recently successfully finished it. Therefore, the direct dynamics of the GFP ESPT became possible at reasonable computational costs only by the end of the funded period but is a great success by itself. We note, that a proposal was written to EPSRC with the aim to continue investigation of GFP ESPT reaction both in a cluster model and including MM environment.

Therefore, we have slightly deviated from the proposed schedule, spending more effort and time at the development of the DD-vMCG code. Execution of objective 2 was partially postponed to a future work, as we found it more important to make the direct quantum dynamics viable for larger molecular systems than to write a new QM/MM interface (a QM/MM interface to the program QoMMa is already present in the code). However, as a part of the code development project, we have written an interface to the Q-Chem software, which itself has a good QM/MM interface with CHARMM program, allowing us to perform hybrid calculations in future. Also, our cluster model of GFP showed that ESPT is actually found to take place all on a single surface and so the diabatisation procedure is not required. However, as proposed in the objective 1, we have developed and implemented a new on-the-fly diabatisation procedure and a paper is in preparation. Also, we plan to write a paper on the ab initio and preliminary quantum dynamics results of the GFP ESPT model.

Overall, although not all the tasks could be completed within the project lifespan, we have made significant progress within each objective, and expect to meet the project goals over the longer term. In terms of scientific progress the project has been a success, leading to a significant improvement of the DD-vMCG code and methodology of quantum nuclear dynamics that will have a considerable effect on the field of theoretical reaction dynamics.