Coarse-Grained Water Model Development for Accurate Dynamics and Structure Prediction.

Several coarse-graining (CG) methods have been combined to develop a CG model of water capable of the accurate prediction of structure and dynamics properties. The multiscale coarse-graining (MS-CG) method based on force matching and the PDF-based coarse-graining method were used for accurate dynamics prediction. The iterative Boltzmann inversion (IBI) method was added for accurate structure representation. The approach is applied to bulk water, and the results show close reproduction of the CG structure when compared with the reference atomistic data. The combination of MS-CG and IBI methods facilitates the development of CG force fields at different temperatures based on a single MS-CG coarse-graining procedure. The dynamic properties of the CG water model closely match those obtained from the reference atomistic system. The general application of this approach to any existing coarse-graining methods is discussed.

Introduction

The derivation of accurate and efficient force fields for molecular systems are desirable if one is aiming to investigate the behavior of a system at the molecular or atomic scale. Classical molecular dynamics (MD) simulations are commonly used to predict the structure and dynamics of a system at the atomic scale.[1] MD is capable of handling large systems containing millions of atoms and running for up to microseconds. The accuracy of properties predicted with MD simulations depend on the goodness of the implemented atomic force field. For many systems of interest, such as bulk water, ionic liquids, and deep eutectic solvents, adequate accuracy in predicting structure, dynamics, and thermodynamic properties may be achieved only with multibody polarizable force fields.[2−6] Implementation of such force fields into an MD simulation can significantly increase the computational power demand and places constraints on the maximum size of the system that can fit within simulation time constraints.

To address these constraints, coarse-grained (CG) models can be used as an alternative.[7−16] In these models, atoms are combined into groups (CG site), thereby achieving a significant reduction in the number of degrees of freedom that must be included in the simulation. Once the CG sites have been defined, the effective nonbonded forces between these CG sites are derived from all-atom simulation trajectories or fit to experimental data. Bonded potentials between CG sites within the same molecule may be developed using various approaches such as fitting with harmonic potentials, applying a Boltzmann inversion method, or using coarse-graining methods.

Several coarse-graining methods have been used to derive effective nonbonded CG forces (or potentials): (i) fitting the free energy in the system,[7,8] (ii) structure-based methods that reproduce a predefined target structure,[9−13] and (iii) force matching approaches, where the instantaneous CG forces are fitted to the forces from all-atom molecular simulations.[14−17] There are a few common limitations when using CG approaches. The derived CG interactions are the many-body potentials of the mean forces and representation of such potentials in a form of pairwise interactions might eliminate the higher-order terms if no additional treatment is performed.[18] Moreover, a reduction in the degrees of freedom during the CG procedure eliminate the effective friction between CG sites and results in faster CG dynamics.[19] These limitations in the CG model development lead to incorrect dynamics for CG systems.

Discrepancies in CG dynamics are typically addressed by utilizing the Mori–Zwanzig (MZ) formalism.[20,21] (1) In the MZ approach, the equation of motion for the CG sites may be derived according to the generalized Langevin equation (GLE). Frictional forces in GLE may be determined through calculations of velocity–velocity time correlation functions (TCFs) of the CG sites[22] or through the incorporation of hydrodynamic interactions into the CG model[23] from simulations of a higher-resolution model. (2) With an assumption that the friction kernel in the MZ approach is only dependent on the CG site position and is pairwise additive, the GLE equation converts to the dissipative particle dynamics (DPD) method.[24,25] (3) Within the MZ approach, the relative entropy method has been formulated.[26−28] Other dynamic-matching methods such as scaling approaches[29−31] are also available. The PDF-based CG method[23] is used in the present work due to its easy implementation into multiscale coarse-graining (MS-CG) codes and its high accuracy in predicting dynamic properties for CG systems.

More than two dozen water models have been developed for use with computer simulations. Some of these models are simpler[32−34] (account for pair interactions with three atoms present), and some are more complex[35,36] (account for higher-order interactions with virtual atoms present). Each water model was designed to represent the desired properties at a reasonable computational cost. For the coarse-grained simulations, the following coarse-grained water model properties are desired: A CG water model needs to reproduce the structure and dynamic properties accurately; a developed coarse-grained force field must be computationally efficient; a coarse-grained water model should be simple to minimize computational time. Several research studies have been reported for developing CG force fields for a CG water molecule,[18,37−40] focusing on accurate structure prediction.

The present work introduces the development of a new CG water model that predicts structure and dynamic properties for a system of bulk water. The developed model is based on the combination of the MS-CG method[15] and a PDF-based CG method[23] for predicting correct dynamic properties. For accurate structure prediction, the IBI method[10] has been applied. A combination of these three coarse-graining methods delivers a simple CG water model with a pairwise interaction force field that accurately reproduces the structure and dynamic properties of a CG system. The applicability of the proposed approach to the CG force field development at different temperatures is also discussed.

Methodology

Molecular Dynamics Simulation Details

The reference all-atom simulations for bulk water were carried out with TINKER[41] using the AMOEBA polarizable atomic multipole water model.[42] The cutoff distance for nonbonded interactions was 12.0 Å. The cutoff radius for electrostatic interactions was 7.0 Å. As a starting point for the equilibration run, a system of 4000 water molecules provided by the TINKER developers[43] was used. The system was equilibrated using an NPT ensemble via a Martyna–Tuckerman–Klein integrator incorporating a Nose–-Hoover thermostat and barostat.[44] Two production runs were carried out at T = 298 and 363 K in a cubic box with a length size of 49.323 and 49.858 Å, respectively, using a Berendsen thermostat.[45] The production run was conducted for 1 ns with a time step of 1.0 fs. Configuration and force data were saved at every 1.0 ps.

Development of a Coarse-Grained Force Field

A coarse-grained model for the bulk water was developed by representing each water molecule with a single CG site calculated at the center of mass (COM) of the molecule. Coarse-grained nonbonded force fields for bulk water were derived from the all-atom data set of trajectories and forces using the MS-CG method based on force matching.[15,16] A detailed description of the force matching method can be found elsewhere.[15−17]

Dynamics for the CG water model were captured through the implementation of the PDF-based CG method. With this approach, probability distributions of coarse-grained forces were extracted with the MS-CG method in addition to the mean CG forces. With the extraction of these random fluctuation forces, the actual force, F̃(q,t), for each CG site due to interactions with all other CG sites can be written aswhere q is the position vector, F(q) is the CG mean force, and ΔF(q,t) is the probability distribution of the CG force. Equation may be compared to the Langevin equation with hydrodynamic interactionswhere M is the mass of the CG site, v(t) is the velocity vector, ζ is the friction tensor, and η(t) is the noise force vector. This leads to the main idea of the PDF-CG method, that the friction tensor, ζ, may be obtained from the probability distribution of CG forces. It can be shown[23] that the probability distribution force is represented by a normal distribution and can be written aswhere F0(q) is the standard deviation of the probability distribution force and N(0,1) is a normal distribution with a mean of zero and standard deviation of 1. According to Berendsen’s formulation[46] of the second fluctuation–dissipation theorem, the following relation between the friction and noise terms in eq and the fluctuation force in eq may be written aswhere kB is the Boltzmann constant, Tref is the reference temperature, and α and β are coefficients that are determined in an iterative manner. In eq , the coefficient α has units of time and adjusts the magnitude of the friction tensor for accurate dynamics prediction and the coefficient β is dimensionless and controls the magnitude of the fluctuation term for accurate temperature representation. By combining eq and eq , the final Langevin equation with hydrodynamic interaction is written asOnce coefficients α and β are determined, the PDF-CG method with hydrodynamic interactions may be implemented.

Iterative Boltzmann Inversion Method

The CG force fields developed in the previous section were adjusted using an iterative Boltzmann inversion method[10] for the accurate prediction of structural properties. A detailed description of the IBI method can be found elsewhere.[9,10] Briefly, the IBI method is based on a self-consistent adjustment of the current potential to achieve convergence with the reference structure. For nonbonded interactions, the IBI method allows one to derive potentials/forces to match the reference radial distribution function (RDF). Potentials are updated in an iterative manner according towhere Un(r) is the new updated CG pair potential, Uc(r) is the current CG pair potential, r is the separation between CG sites, γ is the damping coefficient to reduce large fluctuations in the potential typically varying from 0.01 to 0.1, kB is the Boltzmann constant, Tref is the reference temperature, g(r) is the RDF obtained from the simulation of Uc(r), and gref(r) is the reference or target RDF from the reference. Then, the potential may be converted to a force, F(q), asFor every iteration, the CG MD simulation was conducted for 2.0 ns with a time step of 1.0 fs using LAMMPS.

Results and Discussion

CG Bulk Water

Coarse-grained molecular dynamics simulations for bulk water were conducted with LAMMPS[42] for 4000 CG water molecules. The cutoff distance for nonbonded interactions was 12.0 Å. The system was equilibrated via a Nose–Hoover thermostat for 1.0 ns. The production run was conducted at T = 298 and 363 K in a cubic box with a length side of 49.323 and 49.858 Å, respectively, for 20 ns with a time step of 4.0 fs. Configuration data were saved every 1.0 ps.

The MS-CG method was applied to the all-atom configurations to develop the CG force field. The RDF predicted from the all-atom MD and MS-CG force fields were compared, as shown in Figure a. The location and magnitude for the first peak at 2.82 Å and the third peak at 6.85 Å closely match the reference all-atom data. The magnitude of the second peak is predicted reasonably well; however, the shape of the second peak, as well as the locations and magnitudes of the first and second minima are not accurately predicted. Such a discrepancy in the structure prediction occurs because MS-CG is not designed to fit the two-body RDF. The MS-CG method is based on the force field optimization process, thus delivering the optimal potential for the force residual.

Figure 1

Radial distribution function for COM water molecules at 298 K from all-atom molecular dynamics data (red solid line), single-site CG model (black dashed line), and corresponding CG forces and potentials: (a) by applying the MS-CG method; (b) by applying the MS-CG+IBI approach; (c) CG force field developed with the MS-CG method only (red solid line) and with the MS-CG+IBI approach (black dashed line); (d) CG potentials developed with the MS-CG method only (red solid line) and with the MS-CG+IBI approach (black dashed line).

To improve the structure representation for the CG system, the IBI method has been applied to the already developed MS-CG force field After only two iterations of the IBI method (eq ), the RDF converges to the reference all-atom data as shown in Figure b. As a result, the CG force field slightly changed when compared to the coarse-grained force field before application of IBI, as shown in Figure c. The changes in the coarse-grained potential are insignificant and within a thickness of the line as shown in Figure d. Thus, a small change in the nonbonded force/potential leads to a significant change in the structure of the system.

The dynamic properties of the CG system were recovered by implementing the PDF-based CG method (eq ). The standard deviation data used in eq was obtained from the MS-CG method and coefficients α and β were found after less than 10 iterations: α = 77.8 fs and β = 201.5. The diffusion coefficient, calculated for CG bulk water, was compared to the experimentally measured value.[46,47]Figure shows the mean square displacement (MSD) versus time and the diffusion coefficient, D, calculated from the MSD data according to Einstein’s relation[48]where ⟨(q(0) – q(t))2⟩ is the MSD in d-dimensional space (here, d = 3) during time interval t.

Figure 2

MSD for bulk water at 298 K reconstructed from the experimentally measured value[49] (red solid line) and calculated from the CG MD simulation (black dashed line) using eq .

The MSD of every CG water molecule was computed for each of these time intervals. The production runs were split into t = 1.0 ns time intervals, and the final MSD data are obtained by averaging over all the time intervals. The diffusion coefficient was calculated by fitting MSD versus time with a linear function. The self-diffusion coefficient calculated for water in the CG MD simulation is found to be DCG = 2.27 × 10–9 m2/s, which is in good agreement with the experimentally measured value.[49]

It is important to highlight the differences between the approach used in this work and an IBI method. First, with the current approach the IBI method is applied to an already developed MS-CG coarse-grained force field. In contrast, in the IBI method the reference target property, such as the RDF or the probability distribution function, is used as a starting point. Such a target property is converted into a potential by using a Boltzmann inversion[9,50]where gref is the reference RDF and Tref is the constrained temperature.

Another significant difference is in the number of IBI iterations needed to develop CG effective potentials using the different methods. With the MS-CG+IBI approach, an initial coarse-grained effective force field is developed in a semiautomatic mode. Then, only a few iterative Boltzmann inversions are needed to complete the CG force field development. On the other hand, when implementing a pure IBI method, a significant number of iterations are required to develop a force field even for a relatively simple system.[37]

The IBI approach may be used to develop a CG force field to accurately match the reference all-atom radial distribution function. Ruhle et al. used the VOTCA package[37] to automate the development of a CG force field for bulk water that accurately reproduces the structure of the reference all-atom system. The method introduced in this work uses a different approach in which one would start with MS-CG and PDF-based CG methods to develop an effective potential/force field with an accurate diffusion coefficient and then apply a few iterations of the IBI method to accurately tune the structure of CG system.[38,39] As a result, a CG water model that accurately predicts dynamic and structural properties has been developed.

CG Force Field Development for Various Temperatures

It is desirable for a newly developed force field to be applicable at different thermodynamic state points. Generally, a CG model may deliver reliable results at the thermodynamic conditions of a reference system only. If one expects the CG model to reproduce thermodynamic properties at different state points, the corresponding corrections to the mean-field potentials need to be incorporated.[51] Although multiple approaches address the temperature transferability in CG models,[40,52−56] they are not fully predictive or generally applicable. The temperature- and phase-transferable ultra-coarse-grained (UCG) models were developed for a variety of systems.[56−62] Wide implementation of UCG models is currently limited due to the high computational demand of the method.

In the present work, a CG model that reproduces a structure at different temperatures has been developed by modifying a coarse-grained potential with the IBI approach. Understandably, a newly developed potential would accurately represent solely the structural properties of the CG system. Additional modifications to the coarse-grained potential would be required for the representation of other properties.

To demonstrate a key aspect of this approach, a CG potential/force field at 363 K is constructed from the CG potential, developed at 298 K. For this purpose, all-atom MD simulations of bulk water were performed at 298 and 363 K, and radial distribution functions for COM of water molecules were calculated as shown in Figure a. A CG potential derived at 298 K (see Figure d, black dashed line) has been used as a starting point in matching the structure of the all-atom system at 363 K; with five IBI iterations a new CG effective potential was derived for the system at 363 K (Figure d, red solid line). The derived CG potential slightly differs from those obtained for the CG system at 298 K. The new CG force field/potential accurately predicts the RDF at 363 K when compared to the all-atom data as shown in Figure b. With the MS-CG+IBI approach, a rapid development of the coarse-grained potential that accurately captures the structure of the CG systems at different temperatures is achieved based on a single MS-CG coarse-graining procedure. The diffusivity of the CG system at 363 K was matched by adjusting coefficients α and β in eq in less than five iterations: α = 42 fs and β = 57.2. The self-diffusion coefficient calculated for water in the CG MD simulation at 363 K is found to be DCG = 7.45 × 10–9 m2/s, which is in good agreement with experimental data.[63]

Figure 3

Results of a CG force field development at 363 K by applying the IBI approach to the CG force field derived at 298 K: (a) radial distribution functions for the COM water molecules from all-atom MD simulation at 363 K (red solid line) and all-atom MD simulation at 298 K (black dashed line); (b) radial distribution function for the COM water molecules from an all-atom MD simulation at 363 K (red solid line) and an MS-CG+IBI approach at 363 K (black dashed line); (c) comparison of the MS-CG+IBI force fields at 363 K (red solid line) and at 298 K (black dashed line); (d) comparison of the MS-CG+IBI potentials at 363 K (red solid line) and at 298 K (black dashed line).

Discussion

The development of a new CG model for bulk water with accurate structure and dynamic representation is possible due to a combination of several coarse-graining techniques. Application of the MS-CG and IBI methods allows CG potential/force field development for accurate structure predictions (see Figure b and Figure b). Implementation of the PDF-based CG method allows accurate prediction of diffusion coefficients for CG systems (see Figure ).

There are two reasons for the improvement in structure predictions with the MS-CG+IBI approach: (1) the use of the IBI structure matching method that was developed specifically for matching structures of the coarse-grained model to the reference structure and (2) the application of the IBI method on top of predeveloped MS-CG force fields requires few iterations when developing nonbonded force fields. The MS-CG method delivers an optimal potential for the force residual and produces reasonable RDF data for a variety of systems. Thus, fast convergence in the RDF data is partially due to the good initial guess of the force field. This assertion is supported by ref (38). Such a quick and accurate development of a CG force field facilitates the development of a CG force field at different temperatures (see Figure ) based on a single MS-CG coarse-graining procedure.

A simple pairwise CG force field cannot fully reproduce the three-body correlations and hydrogen bonds present in the water system. As a result, the MS-CG force field is not able to reproduce an all-atomistic RDF. The results can be significantly improved with an explicit addition of three-body terms in the CG potential.[39] The introduction of explicit three-body correction terms increases the computational expense, potentially making such CG simulations infeasible for many complex systems. Therefore, it is beneficial for some CG systems to find an effective force field that accurately reproduces the all-atomistic RDF to improve computational efficiency. Application of the IBI approach on top of an already developed coarse-grained potential using the MS-CG approach fine-tunes the CG potential for accurate structural property predictions in bulk water (see Figure b and Figure b).

It is understandable that the combination of MS-CG and IBI coarse-graining approaches to achieve accurate structure prediction introduces additional complexity into the development of a coarse-grained potential. If the pure IBI approach was used instead, it could deliver a similar CG potential in accord with Henderson’s theorem.[64] Moreover, the MS-CG approach using pairwise basis sets has connections to the liquid-state theory represented with the Yvon–Bonn–Green equation.[65] When the MS-CG approach is mixed with the IBI approach, the liquid-state theory principle is no longer valid. However, the focus of the present work is to develop a CG water model for an accurate representation of the structure and dynamics of bulk water, whereas the main purpose of the MS-CG approach is to generate the necessary inputs for the dynamic-matching PDF-based CG method. Thus, the loss of a connection of the MS-CG approach to the liquid-state theory does not negatively impact the focus of this work.

When the structure of the system can be obtained from the experiment, for example, with the X-ray scattering technique, the method described in this work may be used to develop a coarse-grained potential that matches structures obtained from the experiment. For example, a number of research groups performed X-ray scattering experiments on bulk water.[66−70] When an oxygen–oxygen radial distribution function predicted from a reference all-atom simulation is similar to the experimental scattering plot, the coarse-grained potential may be derived in the following manner. First, the CG potential is derived using the MS-CG approach based on the all-atomistic data. Next, the IBI approach is applied to the developed CG potential to derive an improved potential that matches the CG model structure predictions with experimental results. In this manner, a bottom-up CG approach is used to develop a CG potential with a fine adjustment of the potential to match experimental data.

Conclusions

In summary, this work introduces a derivation of coarse-grained force fields for accurate structure and dynamic predictions of a CG water model. In this approach, an iterative Boltzmann inversion method is applied on top of an already developed MS-CG force field to improve RDFs for the CG water model when compared to the reference all-atomistic system. The dynamics of the CG water system were recovered with the application of a PDF-based CG method.

The described MS-CG+IBI approach combines the best features of three different coarse-graining methods: (1) systematic derivation of the effective coarse-grained force from explicit atomistic molecular simulation with the MS-CG approach based on a force matching method; (2) application of the IBI method designed to match the structure of the reference atomistic system; (3) implementation of the PDF-based CG method to accurately capture the diffusion coefficient of the coarse-grained model. With such a combination of CG methods an efficient CG water model with accurate structure and dynamics properties has been developed.

The ability of the described method to quickly and accurately modify a CG effective force field to match a difference in the structure properties allow the development of a CG force field at different temperatures based on a single set of atomistic configurations.

In general, the application of iterative methods need not be based solely on the MS-CG method; therefore, the IBI approach may be applied in conjunction with any coarse-graining method. Implementation of the IBI approach is based on coarse-graining methods other than MS-CG is a topic for a future study.