Constant chemical potential-quantum mechanical-molecular dynamics simulations of the graphene-electrolyte double layer.

We present the coupling of two frameworks-the pseudo-open boundary simulation method known as constant potential molecular dynamics simulations (CµMD), combined with quantum mechanics/molecular dynamics (QMMD) calculations-to describe the properties of graphene electrodes in contact with electrolytes. The resulting CµQMMD model was then applied to three ionic solutions (LiCl, NaCl, and KCl in water) at bulk solution concentrations ranging from 0.5 M to 6 M in contact with a charged graphene electrode. The new approach we are describing here provides a simulation protocol to control the concentration of electrolyte solutions while including the effects of a fully polarizable electrode surface. Thanks to this coupling, we are able to accurately model both the electrode and solution side of the double layer and provide a thorough analysis of the properties of electrolytes at charged interfaces, such as the screening ability of the electrolyte and the electrostatic potential profile. We also report the calculation of the integral electrochemical double layer capacitance in the whole range of concentrations analyzed for each ionic species, while the quantum mechanical simulations provide access to the differential and integral quantum capacitance. We highlight how subtle features, such as the adsorption of potassium graphene or the tendency of the ions to form clusters contribute to the ability of graphene to store charge, and suggest implications for desalination.

Cleaving Method for Molecular Crystals and Its Application to Calculation of the Surface Free Energy of Crystalline ß-d-Mannitol at Room Temperature.

Calculation of the surface free energy (SFE) is an important application of the thermodynamic integration (TI) methodology, which was mainly employed for atomic crystals (such as Lennard-Jones or metals). In this work, we present the calculation of the SFE of a molecular crystal using the cleaving technique which we implemented in the LAMMPS molecular dynamics package. We apply this methodology to a crystal of ß-d-mannitol at room temperature and report the results for two types of force fields belonging to the GROMOS family: all atoms and united atoms. The results show strong dependence on the type of force field used, highlighting the need for the development of better force fields to model the surface properties of molecular crystals. In particular, we observe that the united-atoms force field, despite its higher degree of coarse graining compared to the all-atoms force field, produces SFE results in better agreement with the experimental results from inverse gas chromatography measurements.

Dynamically Polarizable Force Fields for Surface Simulations via Multi-output Classification Neural Networks.

We present a general procedure to introduce electronic polarization into classical Molecular Dynamics (MD) force fields using a Neural Network (NN) model. We apply this framework to the simulation of a solid-liquid interface where the polarization of the surface is essential to correctly capture the main features of the system. By introducing a multi-input, multi-output NN and treating the surface polarization as a discrete classification problem, we are able to obtain very good accuracy in terms of quality of predictions. Through the definition of a custom loss function we are able to impose a physically motivated constraint within the NN itself making this model extremely versatile, especially in the modeling of different surface charge states. The NN is validated considering the redistribution of electronic charge density within a graphene based electrode in contact with an aqueous electrolyte solution, a system highly relevant to the development of next generation low-cost supercapacitors. We compare the performances of our NN/MD model against Quantum Mechanics/Molecular Dynamics simulations where we obtain a most satisfactory agreement.

Shuttleworth equation: A molecular simulations perspective.

Even though the study of interfacial phenomena can be traced back to Laplace and was given a solid thermodynamic foundation by Gibbs, it appears that some concepts and relations among them are still causing some confusion and debates in the literature, particularly for interfaces involving solids. In particular, the definitions of the concepts of interfacial tension, free energy, and stress and the relationships between them sometimes lack clarity, and some authors even question their validity. So far, the debates about these relationships, in particular the Shuttleworth equation, have taken place within the framework of classical thermodynamics. In this work, we are offering to look at these concepts within the framework of statistical mechanics, which can be readily tested in Molecular Dynamics (MD) simulations. For a simple one component system of particles interacting via the Lennard-Jones potential, we calculate by the cleaving method the excess free energy of a solid-vacuum interface (solid surface) for systems in different states of tangential strain and compare the results to the calculation of surface stress via the difference of normal and tangential forces at the surface. As a result, we demonstrate consistency, within the statistical uncertainty, of the thermodynamic and statistical mechanical definitions of surface free energy and surface stress and how they are expressed via interaction-dependent quantities in MD simulations.

Systematic derivation of hybrid coarse-grained models.

Molecular dynamics represents a key enabling technology for applications ranging from biology to the development of new materials. However, many real-world applications remain inaccessible to fully resolved simulations due to their unsustainable computational costs and must therefore rely on semiempirical coarse-grained models. Significant efforts have been devoted in the last decade towards improving the predictivity of these coarse-grained models and providing a rigorous justification of their use, through a combination of theoretical studies and data-driven approaches. One of the most promising research efforts is the (re)discovery of the Mori-Zwanzig projection as a generic, yet systematic, theoretical tool for deriving coarse-grained models. Despite its clean mathematical formulation and generality, there are still many open questions about its applicability and assumptions. In this work, we propose a detailed derivation of a hybrid multiscale system, generalizing and further investigating the approach developed in Español [Europhys. Lett. 88, 40008 (2009)10.1209/0295-5075/88/40008]. Issues such as the general coexistence of atoms (fully resolved degrees of freedom) and beads (larger coarse-grained units), the role of the fine-to-coarse mapping chosen, and the approximation of effective potentials are discussed. The theoretical discussion is supported by numerical simulations of a monodimensional nonlinear periodic benchmark system with an open-source parallel Julia code, easily extensible to arbitrary potential models and fine-to-coarse mapping functions. The results presented highlight the importance of introducing, in the macroscopic model, nonconstant fluctuating and dissipative terms, given by the Mori-Zwanzig approach, to correctly reproduce the reference fine-grained results, without requiring ad hoc calibration of interaction potentials and thermostats.

The accuracy of ab initio calculations without ab initio calculations for charged systems: Kriging predictions of atomistic properties for ions in aqueous solutions.

Using the machine learning method kriging, we predict the energies of atoms in ion-water clusters, consisting of either Cl- or Na+ surrounded by a number of water molecules (i.e., without Na+Cl- interaction). These atomic energies are calculated following the topological energy partitioning method called Interacting Quantum Atoms (IQAs). Kriging predicts atomic properties (in this case IQA energies) by a model that has been trained over a small set of geometries with known property values. The results presented here are part of the development of an advanced type of force field, called FFLUX, which offers quantum mechanical information to molecular dynamics simulations without the limiting computational cost of ab initio calculations. The results reported for the prediction of the IQA components of the energy in the test set exhibit an accuracy of a few kJ/mol, corresponding to an average error of less than 5%, even when a large cluster of water molecules surrounding an ion is considered. Ions represent an important chemical system and this work shows that they can be correctly taken into account in the framework of the FFLUX force field.

Geometry Optimization with Machine Trained Topological Atoms.

The geometry optimization of a water molecule with a novel type of energy function called FFLUX is presented, which bypasses the traditional bonded potentials. Instead, topologically-partitioned atomic energies are trained by the machine learning method kriging to predict their IQA atomic energies for a previously unseen molecular geometry. Proof-of-concept that FFLUX's architecture is suitable for geometry optimization is rigorously demonstrated. It is found that accurate kriging models can optimize 2000 distorted geometries to within 0.28 kJ mol-1 of the corresponding ab initio energy, and 50% of those to within 0.05 kJ mol-1. Kriging models are robust enough to optimize the molecular geometry to sub-noise accuracy, when two thirds of the geometric inputs are outside the training range of that model. Finally, the individual components of the potential energy are analyzed, and chemical intuition is reflected in the independent behavior of the three energy terms [Formula: see text](intra-atomic), [Formula: see text] (electrostatic) and [Formula: see text] (exchange), in contrast to standard force fields.

FEREBUS: Highly parallelized engine for kriging training.

FFLUX is a novel force field based on quantum topological atoms, combining multipolar electrostatics with IQA intraatomic and interatomic energy terms. The program FEREBUS calculates the hyperparameters of models produced by the machine learning method kriging. Calculation of kriging hyperparameters (θ and p) requires the optimization of the concentrated log-likelihood L̂(θ,p). FEREBUS uses Particle Swarm Optimization (PSO) and Differential Evolution (DE) algorithms to find the maximum of L̂(θ,p). PSO and DE are two heuristic algorithms that each use a set of particles or vectors to explore the space in which L̂(θ,p) is defined, searching for the maximum. The log-likelihood is a computationally expensive function, which needs to be calculated several times during each optimization iteration. The cost scales quickly with the problem dimension and speed becomes critical in model generation. We present the strategy used to parallelize FEREBUS, and the optimization of L̂(θ,p) through PSO and DE. The code is parallelized in two ways. MPI parallelization distributes the particles or vectors among the different processes, whereas the OpenMP implementation takes care of the calculation of L̂(θ,p), which involves the calculation and inversion of a particular matrix, whose size increases quickly with the dimension of the problem. The run time shows a speed-up of 61 times going from single core to 90 cores with a saving, in one case, of ∼98% of the single core time. In fact, the parallelization scheme presented reduces computational time from 2871 s for a single core calculation, to 41 s for 90 cores calculation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Kriging atomic properties with a variable number of inputs.

A new force field called FFLUX uses the machine learning technique kriging to capture the link between the properties (energies and multipole moments) of topological atoms (i.e., output) and the coordinates of the surrounding atoms (i.e., input). Here we present a novel, general method of applying kriging to chemical systems that do not possess a fixed number of (geometrical) inputs. Unlike traditional kriging methods, which require an input system to be of fixed dimensionality, the method presented here can be readily applied to molecular simulation, where an interaction cutoff radius is commonly used and the number of atoms or molecules within the cutoff radius is not constant. The method described here is general and can be applied to any machine learning technique that normally operates under a fixed number of inputs. In particular, the method described here is also useful for interpolating methods other than kriging, which may suffer from difficulties stemming from identical sets of inputs corresponding to different outputs or input biasing. As a demonstration, the new method is used to predict 54 energetic and electrostatic properties of the central water molecule of a set of 5000, 4 Å radius water clusters, with a variable number of water molecules. The results are validated against equivalent models from a set of clusters composed of a fixed number of water molecules (set to ten, i.e., decamers) and against models created by using a naïve method of treating the variable number of inputs problem presented. Results show that the 4 Å water cluster models, utilising the method presented here, return similar or better kriging models than the decamer clusters for all properties considered and perform much better than the truncated models.

Incorporation of local structure into kriging models for the prediction of atomistic properties in the water decamer.

Machine learning algorithms have been demonstrated to predict atomistic properties approaching the accuracy of quantum chemical calculations at significantly less computational cost. Difficulties arise, however, when attempting to apply these techniques to large systems, or systems possessing excessive conformational freedom. In this article, the machine learning method kriging is applied to predict both the intra-atomic and interatomic energies, as well as the electrostatic multipole moments, of the atoms of a water molecule at the center of a 10 water molecule (decamer) cluster. Unlike previous work, where the properties of small water clusters were predicted using a molecular local frame, and where training set inputs (features) were based on atomic index, a variety of feature definitions and coordinate frames are considered here to increase prediction accuracy. It is shown that, for a water molecule at the center of a decamer, no single method of defining features or coordinate schemes is optimal for every property. However, explicitly accounting for the structure of the first solvation shell in the definition of the features of the kriging training set, and centring the coordinate frame on the atom-of-interest will, in general, return better predictions than models that apply the standard methods of feature definition, or a molecular coordinate frame. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Optimization Algorithms in Optimal Predictions of Atomistic Properties by Kriging.

The machine learning method kriging is an attractive tool to construct next-generation force fields. Kriging can accurately predict atomistic properties, which involves optimization of the so-called concentrated log-likelihood function (i.e., fitness function). The difficulty of this optimization problem quickly escalates in response to an increase in either the number of dimensions of the system considered or the size of the training set. In this article, we demonstrate and compare the use of two search algorithms, namely, particle swarm optimization (PSO) and differential evolution (DE), to rapidly obtain the maximum of this fitness function. The ability of these two algorithms to find a stationary point is assessed by using the first derivative of the fitness function. Finally, the converged position obtained by PSO and DE is refined through the limited-memory Broyden-Fletcher-Goldfarb-Shanno bounded (L-BFGS-B) algorithm, which belongs to the class of quasi-Newton algorithms. We show that both PSO and DE are able to come close to the stationary point, even in high-dimensional problems. They do so in a reasonable amount of time, compared to that with the Newton and quasi-Newton algorithms, regardless of the starting position in the search space of kriging hyperparameters. The refinement through L-BFGS-B is able to give the position of the maximum with whichever precision is desired.

Solvent structuring and its effect on the polymer structure and processability: the case of water-acetone poly-ε-caprolactone mixtures.

One of the most common processes to produce polymer nanoparticles is the solvent-displacement method, in which the polymer is dissolved in a "good" solvent and the solution is then mixed with an "anti-solvent". The polymer processability is therefore determined by its structural and transport properties in solutions of the pure solvents and at the intermediate compositions. In this work, we focus on poly-ε-caprolactone (PCL) which is a biocompatible polymer that finds widespread application in the pharmaceutical and biomedical fields, performing full atomistic molecular dynamics simulations of one PCL chain of different molecular weight in a solution of pure acetone (good solvent), of pure water (antisolvent), and their mixtures. Our simulations reveal that the nanostructuring of one of the solvents in the mixture leads to an unexpected identical polymer structure irrespectively of the concentration of the two solvents. In particular, although in pure solvents the behavior of the polymer is, as expected, very different, at intermediate compositions, the PCL chain shows properties very similar to those found in pure acetone as a result of the clustering of the acetone molecules in the vicinity of the polymer chain. We derive an analytical expression to predict the polymer structural properties in solution at different solvent compositions and show that the solvent clustering affects in an unpredictable way the polymer diffusion coefficient. These findings have important consequences on the optimization of the nanoparticle production process and in the implementation of continuum modeling techniques to model it.

A multiple time step scheme for multiresolved models of macromolecules.

In hybrid particle models where coarse-grained beads and atoms are used simultaneously, two clearly separate time scales are mixed. If such models are used in molecular dynamics simulations, a multiple time step (MTS) scheme can therefore be used. In this manuscript, we propose a simple MTS algorithm which approximates for a specific number of integration steps the slow coarse-grained bead-bead interactions with a Taylor series approximation while the atom-atom ones are integrated every time step. The procedure is applied to a previously developed hybrid model of a melt of atactic polystyrene (di Pasquale, Marchisio, and Carbone, J. Chem. Phys. 2012, 137, 164111). The results show that structure, local dynamics, and free diffusion of the model are not altered by the application of the integration scheme which can confidently be used to simulate multiresolved models of polymer melts.

Mixing atoms and coarse-grained beads in modelling polymer melts.

We present a simple hybrid model for macromolecules where the single molecules are modelled with both atoms and coarse-grained beads. We apply our approach to two different polymer melts, polystyrene and polyethylene, for which the coarse-grained potential has been developed using the iterative Boltzmann inversion procedure. Our results show that it is possible to couple the two potentials without modifying them and that the mixed model preserves the local and the global structure of the melts in each of the case presented. The degree of resolution present in each single molecule seems to not affect the robustness of the model. The mixed potential does not show any bias and no cluster of particles of different resolution has been observed.