Rosa E. Bulo

Recent progress in adaptive multiscale molecular dynamics simulations of soft matter

Understanding mesoscopic phenomena in terms of the fundamental motions of atoms and electrons poses a severe challenge for molecular simulation. This challenge is being met by multiscale modeling techniques that aim to bridge between the microscopic and mesoscopic time and length scales. In such techniques different levels of theory are combined to describe a system at a number of scales or resolutions. Here we review recent advancements in adaptive hybrid simulations, in which the different levels are used in separate spatial domains and matter can diffuse from one region to another with an accompanying resolution change. We discuss what it means to simulate such a system, and how to enact the resolution changes. We show how to construct efficient adaptive hybrid quantum mechanics/molecular mechanics (QM/MM) and atomistic/coarse grain (AA/CG) molecular dynamics methods that use an intermediate healing region to smoothly couple the regions together. This coupling is formulated to use only the native forces inherent to each region. The total energy is conserved through the use of auxiliary bookkeeping terms. Error control, and the choice of time step and healing region width, is obtained by careful analysis of the energy flow between the different representations. We emphasize the CG → AA reverse mapping problem and show how this problem is resolved through the use of rigid atomistic fragments located within each CG particle whose orientation is preconditioned for a possible resolution change through a rotational dynamics scheme. These advancements are shown to enable the adaptive hybrid multiscale molecular dynamics simulation of macromolecular soft matter systems.

Toward a Practical Method for Adaptive QM/MM Simulations

We present an accurate adaptive multiscale molecular dynamics method that will enable the detailed study of large molecular systems that mimic experiment. The method treats the reactive regions at the quantum mechanical level and the inactive environment regions at lower levels of accuracy, while at the same time molecules are allowed to flow across the border between active and environment regions. Among many other things, this scheme affords accurate investigation of chemical reactions in solution. A scheme like this ideally fulfills the key criteria applicable to all molecular dynamics simulations: energy conservation and computational efficiency. Approaches that fulfill both criteria can, however, result in complicated potential energy surfaces, creating rapid energy changes when the border between regions is crossed. With the difference-based adaptive solvation potential, a simple approach is introduced that meets the above requirements and reduces fast fluctuations in the potential to a minimum. In cases where none of the current adaptive QM/MM potentials are able to properly describe the system under investigation, we use a continuous force scheme instead, which, while no longer energy conserving, still retains a related conserved quantity along the trajectory. We show that this scheme does not introduce a significant temperature drift on time scales feasible for QM/MM simulations.

PyADF — A scripting framework for multiscale quantum chemistry

Applications of quantum chemistry have evolved from single or a few calculations to more complicated workflows, in which a series of interrelated computational tasks is performed. In particular multiscale simulations, which combine different levels of accuracy, typically require a large number of individual calculations that depend on each other. Consequently, there is a need to automate such workflows. For this purpose we have developed PYADF, a scripting framework for quantum chemistry. PYADF handles all steps necessary in a typical workflow in quantum chemistry and is easily extensible due to its object‐oriented implementation in the Python programming language. We give an overview of the capabilities of PYADF and illustrate its usefulness in quantum‐chemical multiscale simulations with a number of examples taken from recent applications.

NMR solvent shifts of acetonitrile from frozen density embedding calculations.

We present a density functional theory (DFT) study of solvent effects on nuclear magnetic shielding parameters. As a test example we have focused on the sensitive nitrogen shift of acetonitrile immersed in a selected set of solvents, namely water, chloroform, and cyclohexane. To include the effect of the solvent environment in an accurate and efficient manner, we employed the frozen-density embedding (FDE) scheme. We have included up to 500 solvent molecules in the NMR computations and obtained the cluster geometries from a large set of conformations generated with molecular dynamics. For small solute-solvent clusters comparison of the FDE results with conventional supermolecular DFT calculations shows close agreement. For the large solute-solvent clusters the solvent shift values are compared with experimental data and with values obtained using continuum solvent models. For the water --> cyclohexane shift the obtained value is in very good agreement with experiments. For the water --> chloroform NMR solvent shift the classical force field used in the molecular dynamics simulations is found to introduce an error. This error can be largely avoided by using geometries taken from Car-Parrinello molecular dynamics simulations.

Multiscale Modeling of Chemistry in Water: Are We There Yet?

This paper critically evaluates the state of the art in combined quantum mechanical/molecular mechanical (QM/MM) approaches to the computational description of chemistry in water and supplies guidelines for the setup of customized multiscale simulations of aqueous processes. We differentiate between structural and dynamic performance, since some tasks, e.g., the reproduction of NMR or UV-vis spectra, require only structural accuracy, while others, i.e., reaction mechanisms, require accurate dynamic data as well. As a model system for aqueous solutions in general, the approaches were tested on a QM water cluster in an environment of MM water molecules. The key difficulty is the description of the possible diffusion of QM molecules into the MM region and vice versa. The flexible inner region ensemble separator (FIRES) approach constrains QM solvent molecules within an active (QM) region. Sorted adaptive partitioning (SAP), difference-based adaptive solvation (DAS), and buffered-force (BF) are all adaptive approaches that use a buffer zone in which solvent molecules gradually adapt from QM to MM (or vice versa). The costs of SAP and DAS are relatively high, while BF is fast but sacrifices conservation of both energy and momentum. Simulations in the limit of an infinitely small buffer zone, where DAS and SAP become equivalent, are discussed as well and referred to as ABRUPT. The best structural accuracy is obtained with DAS, BF, and ABRUPT, all three of similar quality. FIRES performs very well for dynamic properties localized deep within the QM region. By means of elimination DAS emerges as the best overall compromise between structural and dynamic performance. Eliminating the buffer zone (ABRUPT) improves efficiency and still leads to surprisingly good results. While none of the many new flavors are perfect, all together this new field already allows accurate description of a wide range of structural and dynamic properties of aqueous solutions.

Proton Transfer in Aqueous Solution: Exploring the Boundaries of Adaptive QM/MM

In this chapter, we review the current state-of-the-art in quantum mechanical/molecular mechanical (QM/MM) simulations of reactions in aqueous solutions, and we discuss how proton transfer poses new challenges for its successful application. In the QM/MM description of an aqueous reaction, solvent molecules in the QM region are diffusive and need to be either constrained within the region, or their description (QM versus MM) needs to be updated as they diffuse away. The latter approach is known as adaptive QM/MM. We review several constrained and adaptive QM/MM methods, and classify them in a consistent manner. Most of the adaptive methods employ a transition region, where every solvent molecule can continuously change character (from QM to MM, and vice versa), temporarily becoming partially QM and partially MM. Where a conventional QM/MM scheme partitions a system into a set of QM and a set of MM atoms, an adaptive method employs multiple QM/MM partitions, to describe the fractional QM character. We distinguish two classes of adaptive methods: Discontinuous and continuous. The former methods use at most two QM/MM partitions, and cannot completely avoid discontinuities in the energy and the forces. The more recent continuous adaptive methods employ a larger number of QM/MM partitions for a given configuration. Comparing the performance of the methods for the description of solution chemistry, we find that in certain cases the low-cost constrained methods are sufficiently accurate. For more demanding purposes, the continuous adaptive schemes provide a good balance between dynamical and structural accuracy. Finally, we challenge the adaptive approach by applying it to the difficult topic of proton transfer and diffusion. We present new results, using a well-behaved continuous adaptive method (DAS) to describe an alkaline aqueous solution of methanol. Comparison with fully QM and fully MM simulations shows that the main discrepancies are rooted in the presence of a QM/MM boundary, and not in the adaptive scheme. An anomalous confinement of the hydroxide ion to the QM part of the system stems from the mismatch between QM and MM potentials, which affects the free diffusion of the ion. We also observe an increased water density inside the QM region, which originates from the different chemical potentials of the QM and MM water molecules. The high density results in locally enhanced proton transfer rates.

Energy extrapolation schemes for adaptive multi-scale molecular dynamics simulations

This paper evaluates simple schemes to extrapolate potential energy values using the set of energies and forces extracted from a molecular dynamics trajectory. In general, such a scheme affords the maximum amount of information about a molecular system at minimal computational cost. More specifically, schemes like this are very important in the field of adaptive multi-scale molecular dynamics simulations. In this field, often the computation of potential energy values at certain trajectory points is not required for the simulation itself, but solely for the a posteriori analysis of the simulation data. Extrapolating the values at these points from the available data can save considerable computational time. A set of extrapolation schemes are employed based on Taylor series and central finite difference approximations. The schemes are first tested on the trajectories of molecular systems of varying sizes, obtained at MM and QM level using velocity-Verlet integration with standard simulation time steps. Remarkably good accuracy was obtained with some of the approximations, while the failure of others can be explained in terms of the distinct features of a molecular dynamics trajectory. We have found that, for a Taylor expansion of the potential energy, both a first and a second order truncation exhibit errors that grow with system size. In contrast, the second order central finite difference approximation displays an accuracy that is independent of the size of the system, while giving a very good estimate of the energy, and costing as little as a first order truncation of the Taylor series. A fourth order central finite difference approximation requires more input data, which is not always available in adaptive multi-scale simulations. Furthermore, this approximation gives errors of similar magnitude or larger than its second order counterpart, at standard simulation time steps. This leads to the conclusion that a second order central finite difference approximation is the optimal choice for energy extrapolation from molecular dynamics trajectories. This finding is confirmed in a final application to the analysis of an adaptive multi-scale simulation.

Protonated thiophene-based oligomers as formed within zeolites: understanding their electron delocalization and aromaticity†

In an earlier work, protonated thiophene-based oligomers were identified inside ZSM-5 zeolites. The novel compounds exhibited π-π* absorption wavelengths deep within the visible region, earmarking them for possible use as chromophores in a variety of applications. In this computational study, we determine the factors that cause such low-energy transitions, and describe the electronic structure of these remarkable compounds. DFT calculations of conjugated thiophene-based oligomers with up to five monomer units reveal that the main absorption band of each protonated oligomer is strongly red-shifted compared to the unprotonated form. This effect is counterintuitive, since protonation is expected to diminish aromaticity, and thereby increase the HOMO-LUMO gap. We find that upon protonation the π-electrons remain delocalized over the entire π-conjugated molecule, but the positive charge is localized predominantly on the protonated side of the molecule. A possible explanation for this ground-state charge localization is the participation of the C-H bond in the π-system of the protonated ring, locally providing aromatic stabilization for the positive charge. The addition of the proton stabilizes all electronic orbitals, but due to the ground state π-electron distribution away from the added nucleus, the HOMO is stabilized less than the LUMO. The main absorption peak upon protonation corresponds to the charge transfer excitation involving the frontier orbitals, and the small band gap explains the observed red shift. Analogue calculations on thiophene within a ZSM-5 zeolite cluster model confirm the same trends upon protonation as observed in the non-interacting compounds. Understanding the electronic structure of these compounds is very relevant to correlate UV-Vis bands with acidic strength and possibly environment in zeolites and to improve their performance in catalytic and energy related applications.

Explicit Solvation Matters: Performance of QM/MM Solvation Models in Nucleophilic Addition

Nucleophilic addition onto a carbonyl moiety is strongly affected by solvent, and correctly simulating this solvent effect is often beyond the capability of single-scale quantum mechanical (QM) models. This work explores multiscale approaches for the description of the reversible and highly solvent-sensitive nucleophilic N|···C=O bond formation in an Me2N–(CH2)3–CH=O molecule. In the first stage of this work, we rigorously compare and test four recent quantum mechanical/molecular mechanical (QM/MM) explicit solvation models, employing a QM description of water molecules in spherical regions around both the oxygen and the nitrogen atom of the solute. The accuracy of the models is benchmarked against a reference QM simulation, focusing on properties of the solvated Me2N–(CH2)3–CH=O molecule in its ring-closed form. In the second stage, we select one of the models (continuous adaptive QM/MM) and use it to obtain a reliable free energy profile for the N|···C bond formation reaction. We find that the dual-sphere approach allows the model to accurately account for solvent reorganization along the entire reaction path. In contrast, a simple microsolvation model cannot adapt to the changing conditions and provides an incorrect description of the reaction process.

Accurate Quantum Mechanics/Molecular Mechanics Simulation of Aqueous Solutions with Tailored Molecular Mechanics Models

In recent years, quantum mechanical/molecular mechanical (QM/MM) methods have emerged that are designed specifically for chemical reactions in water. Despite the many advances, a remaining problem is that the patchwork of QM and MM descriptions changes the solvent structure. In a solvent as intricately connected as water, such structural changes can alter a chemical process even across large distances. Examples of structural artifacts in QM/MM water include density accumulation at the QM/MM boundary, decreased order, and density differences between regions. These issues are mostly apparent if the difference between the QM and the MM model is very large, which is often the case with water models. Here, we assess the QM/MM performance of simple MM models that are specifically parametrized to match selected data from a QM simulation of bulk water. To this end, we introduce a novel MM model (PM6-(DH+)-EFF) that reproduces PM6-DH+ water properties. We also assess a recent PBE-DFT-based MM model (PBE-EFF) that reproduces structural properties of bulk water simulated with PBE-DFT. Both models consist solely of tabulated potential energy terms for interactions between atom pairs. We compare the matched QM/MM results (PBE-DFT/PBE-EFF and PM6-DH+/PM6(-DH+)-EFF) with those from mismatched QM/MM simulations (PM6-DH+/PBE-EFF). The mismatched simulations reflect issues similar to those reported for other mismatched QM/MM pairs. The matched simulations yield very good results with water structures that barely deviate from the QM reference. In view of these findings, we strongly recommend adoption of specifically parametrized MM models in the QM/MM simulation of chemical processes in water.