Linear scaling computation of forces for the domain-decomposition linear Poisson-Boltzmann method.

The Linearized Poisson-Boltzmann (LPB) equation is a popular and widely accepted model for accounting solvent effects in computational (bio-) chemistry. In the present article, we derive the analytical forces using the domain-decomposition-based LPB-method with a van-der Waals or solvent-accessible surface. We present an efficient strategy to compute the forces and its implementation, allowing linear scaling of the method with respect to the number of atoms using the fast multipole method. Numerical tests illustrate the accuracy of the computation of the analytical forces and compare the efficiency with other available methods.

Multi-center decomposition of molecular densities: A mathematical perspective.

The aim of this article is to analyze from a mathematical perspective some existing schemes to partition a molecular density into several atomic contributions with a specific focus on Iterative Stockholder Atom (ISA) methods. We provide a unified mathematical framework to describe the latter family of methods and propose a new scheme, named L-ISA (for linear approximation of ISA), which generalizes the so-called additive variational Hirshfeld method. We prove several important mathematical properties of the ISA and L-ISA minimization problems and show that the so-called ISA algorithms can be viewed as alternating minimization schemes, which, in turn, enables us to obtain new convergence results for these numerical methods. Specific mathematical properties of the ISA decomposition for diatomic systems are also presented. Numerical results on diatomic systems illustrate the proven mathematical properties.

Subdiffusive-Brownian crossover in membrane proteins: a generalized Langevin equation-based approach.

In this work, we propose a generalized Langevin equation-based model to describe the lateral diffusion of a protein in a lipid bilayer. The memory kernel is represented in terms of a viscous (instantaneous) and an elastic (noninstantaneous) component modeled through a Dirac δ function and a three-parameter Mittag-Leffler type function, respectively. By imposing a specific relationship between the parameters of the three-parameter Mittag-Leffler function, the different dynamical regimes-namely ballistic, subdiffusive, and Brownian, as well as the crossover from one regime to another-are retrieved. Within this approach, the transition time from the ballistic to the subdiffusive regime and the spectrum of relaxation times underlying the transition from the subdiffusive to the Brownian regime are given. The reliability of the model is tested by comparing the mean-square displacement derived in the framework of this model and the mean-square displacement of a protein diffusing in a membrane calculated through molecular dynamics simulations.

Theoretical analysis of screened many-body electrostatic interactions between charged polarizable particles.

This paper builds on two previous studies [Lindgren et al., J. Comput. Phys. 371, 712 (2018) and Quan et al., "A domain decomposition method for the Poisson-Boltzmann solvation models," SIAM J. Sci. Comput. (to be published); e-print arXiv:1807.05384] to devise a new method to solve the problem of calculating electrostatic interactions in a system composed by many dielectric particles, embedded in a homogeneous dielectric medium, which in turn can also be permeated by charge carriers. The system is defined by the charge, size, position, and dielectric constant of each particle, as well as the dielectric constant and the Debye length of the medium. The effects of taking into account the dielectric nature of the particles are explored in selected scenarios where the presence of electrolytes in the medium can significantly influence the total undergoing interactions. The description of the mutual interactions between all particles in the system as being truly of many-body nature reveals how such effects can effectively influence the magnitudes and even directions of the resulting forces, especially those acting on particles that have a null net charge. Particular attention is given to a situation that can be related to colloidal particles in an electrolyte solution, where it is shown that polarization effects alone can substantially raise or lower-depending on the dielectric contrast between the particles and the medium-the energy barrier that divides particle coagulation and flocculation regions, when an interplay between electrostatic and additional van der Waals forces is considered. Overall, the results suggest that for an accurate description of the type of system in question, it is essential to consider particle polarization if the separation between the interacting particles are comparable to or smaller than the Debye length of the medium.

Dynamic simulations of many-body electrostatic self-assembly.

Two experimental studies relating to electrostatic self-assembly have been the subject of dynamic computer simulations, where the consequences of changing the charge and the dielectric constant of the materials concerned have been explored. One series of calculations relates to experiments on the assembly of polymer particles that have been subjected to tribocharging and the simulations successfully reproduce many of the observed patterns of behaviour. A second study explores events observed following collisions between single particles and small clusters composed of charged particles derived from a metal oxide composite. As before, observations recorded during the course of the experiments are reproduced by the calculations. One study in particular reveals how particle polarizability can influence the assembly process.This article is part of the theme issue 'Modern theoretical chemistry'.

A coherent derivation of the Ewald summation for arbitrary orders of multipoles: The self-terms.

In this work, we provide the mathematical elements we think essential for a proper understanding of the calculus of the electrostatic energy of point-multipoles of arbitrary order under periodic boundary conditions. The emphasis is put on the expressions of the so-called self-parts of the Ewald summation where different expressions can be found in the literature. Indeed, such expressions are of prime importance in the context of new generation polarizable force field where the self-field appears in the polarization equations. We provide a general framework, where the idea of the Ewald splitting is applied to the electric potential and, subsequently, all other quantities such as the electric field, the energy, and the forces are derived consistently thereof. Mathematical well-posedness is shown for all these contributions for any order of multipolar distribution.

Computation of forces arising from the polarizable continuum model within the domain-decomposition paradigm.

The domain-decomposition (dd) paradigm, originally introduced for the conductor-like screening model, has been recently extended to the dielectric Polarizable Continuum Model (PCM), resulting in the ddPCM method. We present here a complete derivation of the analytical derivatives of the ddPCM energy with respect to the positions of the solute's atoms and discuss their efficient implementation. As it is the case for the energy, we observe a quadratic scaling, which is discussed and demonstrated with numerical tests.

A new discretization for the polarizable continuum model within the domain decomposition paradigm.

We present a new algorithm to solve the polarizable continuum model equation in a framework compatible with the strategy previously developed by us for the conductor-like screening model based on Schwarz's domain decomposition method (ddCOSMO). The new discretization is systematically improvable and is fully consistent with ddCOSMO so that it reproduces ddCOSMO results for large dielectric constants.

Quantum, classical, and hybrid QM/MM calculations in solution: general implementation of the ddCOSMO linear scaling strategy.

We present the general theory and implementation of the Conductor-like Screening Model according to the recently developed ddCOSMO paradigm. The various quantities needed to apply ddCOSMO at different levels of theory, including quantum mechanical descriptions, are discussed in detail, with a particular focus on how to compute the integrals needed to evaluate the ddCOSMO solvation energy and its derivatives. The overall computational cost of a ddCOSMO computation is then analyzed and decomposed in the various steps: the different relative weights of such contributions are then discussed for both ddCOSMO and the fastest available alternative discretization to the COSMO equations. Finally, the scaling of the cost of the various steps with respect to the size of the solute is analyzed and discussed, showing how ddCOSMO opens significantly new possibilities when cheap or hybrid molecular mechanics/quantum mechanics methods are used to describe the solute.

The significance of multipole interactions for the stability of regular structures composed from charged particles.

Identifying the forces responsible for stabilising binary particle lattices is key to the controlled fabrication of many new materials. Experiments have shown that the presence of charge can be integral to the formation of ordered arrays; however, a complete analysis of the forces responsible has not included many of the significant lattice types that may form during fabrication. A theory of many-body electrostatic interactions has been applied to six lattice stoichiometries, AB, AB2, AB3, AB4, AB5 and AB6, to show that induced multipole interactions can make a very significant (>80 %) contribution to the total lattice energy of arrays of charged particles. Particle radii ratios which favour global minima in electrostatic energy are found to be the same or a close match to those observed by experiment. Although certain lattice types exhibit local energy minima, the calculations show that many-body rather than two-body interactions are ultimately responsible for the structures observed by experiment. For a lattice isostructural with CFe4, a particle size ratio not previously observed is found to be particularly stable due to many-body effects.

Domain decomposition for implicit solvation models.

This article is the first of a series of papers dealing with domain decomposition algorithms for implicit solvent models. We show that, in the framework of the COSMO model, with van der Waals molecular cavities and classical charge distributions, the electrostatic energy contribution to the solvation energy, usually computed by solving an integral equation on the whole surface of the molecular cavity, can be computed more efficiently by using an integral equation formulation of Schwarz's domain decomposition method for boundary value problems. In addition, the so-obtained potential energy surface is smooth, which is a critical property to perform geometry optimization and molecular dynamics simulations. The purpose of this first article is to detail the methodology, set up the theoretical foundations of the approach, and study the accuracies and convergence rates of the resulting algorithms. The full efficiency of the method and its applicability to large molecular systems of biological interest is demonstrated elsewhere.

Reduced basis surrogates for quantum spin systems based on tensor networks.

Within the reduced basis methods approach, an effective low-dimensional subspace of a quantum many-body Hilbert space is constructed in order to investigate, e.g., the ground-state phase diagram. The basis of this subspace is built from solutions of snapshots, i.e., ground states corresponding to particular and well-chosen parameter values. Here, we show how a greedy strategy to assemble the reduced basis and thus to select the parameter points can be implemented based on matrix-product-state calculations. Once the reduced basis has been obtained, observables required for the computation of phase diagrams can be computed with a computational complexity independent of the underlying Hilbert space for any parameter value. We illustrate the efficiency and accuracy of this approach for different one-dimensional quantum spin-1 models, including anisotropic as well as biquadratic exchange interactions, leading to rich quantum phase diagrams.

A Quasi Time-Reversible Scheme Based on Density Matrix Extrapolation on the Grassmann Manifold for Born-Oppenheimer Molecular Dynamics.

This Letter introduces the so-called Quasi Time-Reversible scheme based on Grassmann extrapolation (QTR G-Ext) of density matrices for an accurate calculation of initial guesses in Born-Oppenheimer Molecular Dynamics (BOMD) simulations. The method shows excellent results on four large molecular systems that are representative of real-life production applications, ranging from 21 to 94 atoms simulated with Kohn-Sham (KS) density functional theory surrounded with a classical environment with 6k to 16k atoms. Namely, it clearly reduces the number of self-consistent field iterations while at the same time achieving energy-conserving simulations, resulting in a considerable speed-up of BOMD simulations even when tight convergence of the KS equations is required.

Surrogate models for quantum spin systems based on reduced-order modeling.

We present a methodology to investigate phase diagrams of quantum models based on the principle of the reduced basis method (RBM). The RBM is built from a few ground-state snapshots, i.e., lowest eigenvectors of the full system Hamiltonian computed at well-chosen points in the parameter space of interest. We put forward a greedy strategy to assemble such a small-dimensional basis, i.e., to select where to spend the numerical effort needed for the snapshots. Once the RBM is assembled, physical observables required for mapping out the phase diagram (e.g., structure factors) can be computed for any parameter value with a modest computational complexity, considerably lower than the one associated to the underlying Hilbert space dimension. We benchmark the method in two test cases, a chain of excited Rydberg atoms and a geometrically frustrated antiferromagnetic two-dimensional lattice model, and illustrate the accuracy of the approach. In particular, we find that the ground-state manifold can be approximated to sufficient accuracy with a moderate number of basis functions, which increases very mildly when the number of microscopic constituents grows-in stark contrast to the exponential growth of the Hilbert space needed to describe each of the few snapshots. A combination of the presented RBM approach with other numerical techniques circumventing even the latter big cost, e.g., tensor network methods, is a tantalizing outlook of this work.

Coarse-Graining ddCOSMO through an Interface between Tinker and the ddX Library.

The domain decomposition conductor-like screening model is an efficient way to compute the solvation energy of solutes within a polarizable continuum medium in a linear scaling computational time. Despite its efficiency, the application to very large systems is still challenging. A possibility to further accelerate the algorithm is resorting to coarse-graining strategies. In this paper we present a preliminary interface between the molecular dynamics package Tinker and the ddX library. The interface was used to test a united atom coarse-graining strategy that allowed us to push ddCOSMO to its limits by computing solvation energies on systems with up to 7 million atoms. We first present benchmarks to find an optimal discretization, and then, we discuss the performance and results obtained with fine- and coarse-grained solvation energy calculations.

Manipulating Interactions between Dielectric Particles with Electric Fields: A General Electrostatic Many-Body Framework.

We derive a rigorous analytical formalism and propose a numerical method for the quantitative evaluation of the electrostatic interactions between dielectric particles in an external electric field. This formalism also allows for inhomogeneous charge distributions, and, in particular, for the presence of pointlike charges on the particle surface. The theory is based on a boundary integral equation framework and yields analytical expressions for the interaction energy and net forces that can be computed in linear scaling cost, with respect to the number of interacting particles. We include numerical results that validate the proposed method and show the limitations of the fixed dipole approximation at small separation between interacting particles. The proposed method is also applied to study the stability and melting of ionic colloidal crystals in an external electric field.

Grassmann Extrapolation of Density Matrices for Born-Oppenheimer Molecular Dynamics.

Born-Oppenheimer molecular dynamics (BOMD) is a powerful but expensive technique. The main bottleneck in a density functional theory BOMD calculation is the solution to the Kohn-Sham (KS) equations that requires an iterative procedure that starts from a guess for the density matrix. Converged densities from previous points in the trajectory can be used to extrapolate a new guess; however, the nonlinear constraint that an idempotent density needs to satisfy makes the direct use of standard linear extrapolation techniques not possible. In this contribution, we introduce a locally bijective map between the manifold where the density is defined and its tangent space so that linear extrapolation can be performed in a vector space while, at the same time, retaining the correct physical properties of the extrapolated density using molecular descriptors. We apply the method to real-life, multiscale, polarizable QM/MM BOMD simulations, showing that sizeable performance gains can be achieved, especially when a tighter convergence to the KS equations is required.

Quantum Calculations in Solution of Energies, Structures, and Properties with a Domain Decomposition Polarizable Continuum Model.

In this work, we present the first implementation of the domain decomposition polarizable continuum model for a solute described at a quantum mechanical level of theory. After briefly recapitulating the theory, we discuss the coupling of ddPCM to a quantum mechanical level of theory based on the self-consistent field approach, i.e., Hartree-Fock, density functional theory, and semiempirical methods. We then present benchmarks of the new implementation, comparing it to a currently available state-of-the-art one, and use it to describe the structure and excitation properties of a large multichromophoric system.

Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields.

We present Tinker-HP, a massively MPI parallel package dedicated to classical molecular dynamics (MD) and to multiscale simulations, using advanced polarizable force fields (PFF) encompassing distributed multipoles electrostatics. Tinker-HP is an evolution of the popular Tinker package code that conserves its simplicity of use and its reference double precision implementation for CPUs. Grounded on interdisciplinary efforts with applied mathematics, Tinker-HP allows for long polarizable MD simulations on large systems up to millions of atoms. We detail in the paper the newly developed extension of massively parallel 3D spatial decomposition to point dipole polarizable models as well as their coupling to efficient Krylov iterative and non-iterative polarization solvers. The design of the code allows the use of various computer systems ranging from laboratory workstations to modern petascale supercomputers with thousands of cores. Tinker-HP proposes therefore the first high-performance scalable CPU computing environment for the development of next generation point dipole PFFs and for production simulations. Strategies linking Tinker-HP to Quantum Mechanics (QM) in the framework of multiscale polarizable self-consistent QM/MD simulations are also provided. The possibilities, performances and scalability of the software are demonstrated via benchmarks calculations using the polarizable AMOEBA force field on systems ranging from large water boxes of increasing size and ionic liquids to (very) large biosystems encompassing several proteins as well as the complete satellite tobacco mosaic virus and ribosome structures. For small systems, Tinker-HP appears to be competitive with the Tinker-OpenMM GPU implementation of Tinker. As the system size grows, Tinker-HP remains operational thanks to its access to distributed memory and takes advantage of its new algorithmic enabling for stable long timescale polarizable simulations. Overall, a several thousand-fold acceleration over a single-core computation is observed for the largest systems. The extension of the present CPU implementation of Tinker-HP to other computational platforms is discussed.

Meshing molecular surfaces based on analytical implicit representation.

We develop an algorithm for meshing molecular surfaces that is based on patch-wise meshing using an advancing-front method adapted to the particular case of molecular surfaces. We focus on the solvent accessible surface (SAS) and the solvent excluded surface (SES). The essential ingredient is a newly developed analysis of such surfaces in [18] that allows to describe all SES-singularities a priori and therefore a complete characterization of the SES. In addition, an algorithm for filling molecular inner holes is proposed based on the pre-computed data structures of molecular surfaces.