Mass-zero constrained dynamics for simulations based on orbital-free density functional theory.

A new algorithm for efficient and fully time-reversible integration of first-principles molecular dynamics based on orbital-free density functional theory (OFDFT) is presented. The algorithm adapts to this nontrivial case, the recently introduced Mass-Zero (MaZe) constrained dynamics. The formalism ensures that full adiabatic separation is enforced between nuclear and electronic degrees of freedom and, consequently, that the exact Born-Oppenheimer probability for the nuclei is sampled. Numerical integration of the MaZe dynamics combines standard molecular dynamics algorithms, e.g., Verlet or velocity Verlet, with the SHAKE method to impose the minimum conditions on the electronic degrees of freedom as a set of constraints. The developments presented in this work, which include a bespoke adaptation of the standard SHAKE algorithm, ensure that the quasilinear scaling of OFDFT is preserved by the new method for a broad range of kinetic and exchange-correlation functionals, including nonlocal ones. The efficiency and accuracy of the approach are demonstrated via calculations of static and dynamic properties of liquid sodium in the constant energy and constant temperature ensembles.

Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.

Classical molecular dynamics simulations have recently become a standard tool for the study of electrochemical systems. State-of-the-art approaches represent the electrodes as perfect conductors, modeling their responses to the charge distribution of electrolytes via the so-called fluctuating charge model. These fluctuating charges are additional degrees of freedom that, in a Born-Oppenheimer spirit, adapt instantaneously to changes in the environment to keep each electrode at a constant potential. Here, we show that this model can be treated in the framework of constrained molecular dynamics, leading to a symplectic and time-reversible algorithm for the evolution of all the degrees of freedom of the system. The computational cost and the accuracy of the new method are similar to current alternative implementations of the model. The advantage lies in the accuracy and long term stability guaranteed by the formal properties of the algorithm and in the possibility to systematically introduce additional kinematic conditions of arbitrary number and form. We illustrate the performance of the constrained dynamics approach by enforcing the electroneutrality of the electrodes in a simple capacitor consisting of two graphite electrodes separated by a slab of liquid water.

Characterization of the Photochemical Properties of 5-Benzyluracil via Time-Dependent Density Functional Theory.

We present a detailed study of the excited state properties of 5-benzyluracil (5BU) in the gas phase and in implicit solvent using different electronic structure approaches ranging from time-dependent density functional theory in the linear response regime (LR-TDDFT) to a set of different wave-function-based methods for excited states, namely perturbed coupled cluster (CC2), algebraic diagrammatic construction method to second order (ADC(2)), and perturbed configuration interaction (CIS(D)). 5BU has been used to investigate DNA base-amino acid interactions. In particular, it served as a model of protein-DNA photoinduced cross-linking. While LR-TDDFT is computationally the most efficient first-principles approach for static and dynamic simulations of this bichromophoric system, its accuracy is difficult to assess due to the presence of excited states with charge transfer character. In this work, the performance of different exchange correlation functionals is compared against accurate benchmarks obtained either from high level wave-function-based methods or directly from experimental absorption spectra. Our investigation shows that accurate results for the excitation energies can be obtained using the hybrid meta-GGA functional M06. In view of dynamical studies of the relaxation of 5BU after photoexcitation, we also show that the PBE functional, while failing in the Franck-Condon region, provides qualitatively good results for the characterisation of a possible photocyclization path.

Computing thermal Wigner densities with the phase integration method.

We discuss how the Phase Integration Method (PIM), recently developed to compute symmetrized time correlation functions [M. Monteferrante, S. Bonella, and G. Ciccotti, Mol. Phys. 109, 3015 (2011)], can be adapted to sampling/generating the thermal Wigner density, a key ingredient, for example, in many approximate schemes for simulating quantum time dependent properties. PIM combines a path integral representation of the density with a cumulant expansion to represent the Wigner function in a form calculable via existing Monte Carlo algorithms for sampling noisy probability densities. The method is able to capture highly non-classical effects such as correlation among the momenta and coordinates parts of the density, or correlations among the momenta themselves. By using alternatives to cumulants, it can also indicate the presence of negative parts of the Wigner density. Both properties are demonstrated by comparing PIM results to those of reference quantum calculations on a set of model problems.

Path integral based calculations of symmetrized time correlation functions. I.

In this paper, we examine how and when quantum evolution can be approximated in terms of (generalized) classical dynamics in calculations of correlation functions, with a focus on the symmetrized time correlation function introduced by Schofield. To that end, this function is expressed as a path integral in complex time and written in terms of sum and difference path variables. Taylor series expansion of the path integral's exponent to first and second order in the difference variables leads to two original developments. The first order expansion is used to obtain a simple, path integral based, derivation of the so-called Schofield's quantum correction factor. The second order result is employed to show how quantum mechanical delocalization manifests itself in the approximation of the correlation function and hinders, even in the semiclassical limit, the interpretation of the propagators in terms of sets of guiding classical trajectories dressed with appropriate weights.

Path integral based calculations of symmetrized time correlation functions. II.

Schofield's form of quantum time correlation functions is used as the starting point to derive a computable expression for these quantities. The time composition property of the propagators in complex time is exploited to approximate Schofield's function in terms of a sequence of short time classical propagations interspersed with path integrals that, combined, represent the thermal density of the system. The approximation amounts to linearization of the real time propagators and it becomes exact with increasing number of propagation legs. Within this scheme, the correlation function is interpreted as an expectation value over a probability density defined on the thermal and real path space and calculated by a Monte Carlo algorithm. The performance of the algorithm is tested on a set of benchmark problems. Although the numerical effort required is considerable, we show that the algorithm converges systematically to the exact answer with increasing number of iterations and that it is stable for times longer than those accessible via a brute force, path integral based, calculation of the correlation function. Scaling of the algorithm with dimensionality is also examined and, when the method is combined with commonly used filtering schemes, found to be comparable to that of alternative semiclassical methods.

Iterative linearized approach to nonadiabatic dynamics.

This paper presents a new approach to propagating the density matrix based on a time stepping procedure arising from a Trotter factorization and combining the forward and backward incremental propagators. The sums over intermediate states of the discrete quantum subsystem are implemented by a Monte Carlo surface hopping-like procedure, while the integrals over the continuous variables are performed using a linearization in the difference between the forward and backward paths of these variables leading to classical-like equations of motion with forces determined by the quantum subsystem states. The approach is tested on several models and numerical convergence is explored.

Mapping the hydropathy of amino acids based on their local solvation structure.

In spite of its relevant biological role, no general consensus exists on the quantitative characterization of amino acid's hydropathy. In particular, many hydrophobicity scales exist, often producing quite different rankings for the amino acids. To make progress toward a systematic classification, we analyze amino acids' hydropathy based on the orientation of water molecules at a given distance from them as computed from molecular dynamics simulations. In contrast with what is usually done, we argue that assigning a single number is not enough to characterize the properties of an amino acid, in particular when both hydrophobic and hydrophilic regions are present in a residue. Instead we show that appropriately defined conditional probability densities can be used to map the hydrophilic and hydrophobic groups on the amino acids with greater detail than possible with other available methods. Three indicators are then defined based on the features of these probabilities to quantify the specific hydrophobicity and hydrophilicity of each amino acid. The characterization that we propose can be used to understand some of the ambiguities in the ranking of amino acids in the current scales. The quantitative indicators can also be used in combination with standard bioinformatics tools to predict the location of transmembrane regions of proteins. The method is sensitive to the specific environment of the amino acids and can be applied to unnatural and modified amino acids, as well as to other small organic molecules.

LAND-map, a linearized approach to nonadiabatic dynamics using the mapping formalism.

We present a new approach for calculating quantum time correlation functions for systems whose dynamics exhibits relevant nonadiabatic effects. The method involves partial linearization of the full quantum path-integral expression for the time correlation function written in the nonadiabatic mapping Hamiltonian formalism. Our analysis gives an algorithm which is both numerically efficient and accurate as we demonstrate in test calculations on the spin-boson model where we find results in good agreement with exact calculations. The accuracy of our new approach is comparable to that of calculations performed using other approximate methods over a relatively broad range of model parameters. However, our method converges relatively quickly when compared with most alternative schemes. These findings are very encouraging in view of the application of the new method for studying realistic nonadiabatic model problems in the condensed phase.