A polarizable force field for water using an artificial neural network

Abstract

A force field for liquid water including polarization effects has been constructed using an artificial neural network (ANN). It is essential to include a many-body polarization effect explicitly into a potential energy function in order to treat liquid water which is dense and highly polar. The new potential energy function is a combination of empirical and nonempirical potentials. The TIP4P model was used for the empirical part of the potential. For the nonempirical part, an ANN with a back-propagation of error algorithm (BPNN) was introduced to reproduce the complicated many-body interaction energy surface from ab initio quantum mechanical calculations. BPNN, described in terms of a matrix, provides enough flexibility to describe the complex potential energy surface (PES). The structural and thermodynamic properties, calculated by isobaric–isothermal (constant-NPT) Monte Carlo simulations with the new polarizable force field for water, are compatible with experimental results. Thus, the simulation establishes the validity of using our estimated PES with a polarization effect for accurate predictions of liquid state properties. Applications of this approach are simple and systematic so that it can easily be applied to the development of other force fields besides the water–water system.

Introduction

For studies of solution chemistry and biological systems, it is essential to understand the structure and physical properties of liquid water because water plays an essential role in biological systems. From a theoretical point of view, an accurate description of the water molecule is a crucial point in a treatment of the liquid phase [1].

As a monomer, the water molecule is structurally quite simple, containing one oxygen atom and two hydrogen atoms. At the same time, it is quite difficult to treat with molecular mechanics because of its ability to act as both an acceptor and a donor in hydrogen-bonding [1], and because of its high polarizability [2]. The high polarizability leads to a large change in the magnitude and orientation of the dipole of the water molecule when present in a liquid environment. The water molecule is one of the most challenging systems to model precisely, and a wide range of models of liquid water have been proposed and reviewed in the literature [1], [3], [4]. Since many of its structural and physical properties can easily be determined by computer simulation [1], [5], [6], [7], [8], [9], [10], [11], [12], [13] and easily compared with experimental data [14], [15], [16], [17], computer simulations of models of liquid water have become very valuable tools to test the accuracy of such models. Moreover, a very large number of water molecules are involved in the computer simulations, so that the computational efficiency with which the energy can be calculated for a given model is often an important factor.

Rigid geometry models of water are conventionally divided into two types. In the effective models of liquid water, rigid geometry is maintained for the water monomers, and the interactions between molecules are described with pairwise electrostatic and nonbonded interactions. Other models (known as polarizable water models) have been developed to include the effects of polarization explicitly.

The total potential energy, V, of a system of N particles can generally be expressed by summing over the contributions arising from interacting pairs, triplets, etc.

$$\text{V(r}{}_1\text{,r}{}_2\text{,r}{}_3\text{,…,r}{}_{\mathrm{N}}\text{)=}\text{∑}\text{i<j}\text{f}{}_{\mathrm{l}}\text{(r}{}_{\mathrm{i}}\text{,r}{}_{\mathrm{j}}\text{)+}\text{∑}\text{i<j<k}\text{f}{}_{\mathrm{m}}\text{(r}{}_{\mathrm{i}}\text{,r}{}_{\mathrm{j}}\text{,r}{}_{\mathrm{k}}\text{)+⋯}$$

where f_l,f_m,… indicate a two-body, three-body,… interaction, and r describes the coordinates of one molecule. If the contributions from the higher terms to the total potential energy are negligible compared to the two-body interaction terms, the total potential energy can be reduced to a simple form arising solely from pairwise interactions.

$$\text{V(r}{}_1\text{,…,r}{}_{\mathrm{N}}\text{)≈}\text{∑}\text{i<j}\text{f}{}_{\mathrm{l}}\text{(r}{}_{\mathrm{i}}\text{,r}{}_{\mathrm{j}}\text{).}$$

Since the water molecule has a large permanent dipole moment [18] and a high polarizability [2], and the intermolecular distances in condensed systems are very small, the polarization effect becomes a main contributor to f_m and higher terms in Eq. (1), and cannot be ignored in a model for polar and condensed water molecules. An efficient way to include the changes in the electronic charge distribution of the water molecule directly in response to its condensed environment is to assign an effective charge distribution, so that the dipole moment is increased to take account of the surrounding environment. The experimentally determined dipole moment of a water molecule in the gas phase is 1.85 D [18]. The dipole moment of an individual water molecule in the liquid, calculated with the effective models of liquid water, is higher compared to the values in the gas phase; e.g. the dipole moment of the SPC model is 2.27 D [19] and that for the TIP4P model is 2.18 D [20]. These values are closer than the gas phase value to the experimental effective dipole moment of monomeric water in the liquid, which is approximately 2.6 D [21]. The use of effective electrostatic and nonbonded parameters to incorporate many-body effects is more cost effective than explicitly parameterizing and calculating f_m and higher order contributions in Eq. (1). An effective pairwise potential is then evaluated as

$$\text{V(r}{}_1\text{,…,r}{}_{\mathrm{N}}\text{)≈}\text{∑}\text{ij}\text{N}\text{f}{}_{\mathrm{<mtext>eff</mtext>}}\text{(r}{}_{\mathrm{i}}\text{,r}{}_{\mathrm{j}}\text{)}$$

where f_eff now includes all contributions from the higher order terms in Eq. (1) by a ‘mean field’ approximation. Because an effective charge distribution defines all electrostatic moments of the molecule, the correction pertains to infinite order. This approximation is computationally efficient, which was an important consideration when the available computational facilities were limited.

Even though different effective water models have been suggested by several authors [19], [20], [22], [23], [24], [25], [26], the basic components of the models are similar to each other. The electrostatic interactions are described and modeled by an arrangement of point charges located in a well-defined molecular geometry, and the short-range dispersive and repulsive interaction terms are approximated by a Lennard–Jones (LJ) or a Morse potential. A description of the most commonly used effective water potential energy functions is listed in Table 1. The effective water models can be classified according to the number of interaction sites in the water molecule. The TIP3P [20] and SPC [19] models use a total of three sites for the electrostatic interactions; the partial positive charges on the hydrogen atoms are exactly balanced by an appropriate negative charge located on the oxygen atom for a given geometry of the water molecule. The nonbonded interaction between two water molecules is computed with a LJ function with only a single interaction point per molecule centered on the oxygen atom for computational efficiency. The four-site models such as the BF [22], TIP4P [20] and MCY [24] models shift the negative charge from the oxygen atom to a point along the bisector of the HOH bond angle towards the hydrogens and also have one nonbonded interaction site for each molecule. The ST2 [23] and EPEN [25], [26] potentials are five-site models in which the charges are placed on the hydrogen atoms and two lone-pair sites on the oxygen atom and only one nonbonded interaction site for each molecule.

Since the number of interaction sites of the effective models is less than the number which includes polarization or many-body effects explicitly, the simple effective models are computationally efficient, very robust, and have been used widely in biomolecular simulations. Many effective water models, such as SPC, TIP3P, and TIP4P, were developed by fitting liquid state properties, including internal energy, density, and structure to include the many-body effect of water under dense conditions. Thus, the polarization effect is included implicitly in these models. The simple models give good results for a wide range of properties of pure liquid water under normal conditions, viz. 1 atm and room temperature.

Because the parameters for the electrostatic and nonbonded interactions were developed by fitting liquid state properties, and polarization and many-body effects are not included explicitly in these models, their transferability to other physical states or environments such as the low-density gas phase or high-charge-density is limited because the environments of the simulation differ from those of the parameterization. Thus, nonpolarizable water models suffer from a lack of transferability, and, therefore, they cannot provide an accurate description of the phase behavior under conditions other than the one used in the parameterization.

Because of the limited transferability of the simple models, simulations of other thermodynamic properties, such as those at the critical point, and in mixtures with ions or with nonpolar species, can exhibit problems. It is obvious that the electronic configuration of a given water molecule depends explicitly on its environment, and this effect should be included in a water potential for transferability to various conditions.

Pure pair potentials cannot reproduce condensed-state properties for polar molecules because such potentials neglect the polarization effect beyond the level of pair interactions [19], and the effective potentials cannot represent the properties of water except in the state for which it is parameterized. In the case of water, where induction effects change dramatically from low-density states to high-density states, the total energy cannot be expressed as a sum of two-body interactions. Both experimental and computational studies of water from small clusters to bulk phases show that the strong many-body polarization effect is one of the most important features of aqueous systems [21], [27], [28]. Current research [8], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51] on water potentials tends to focus on incorporating explicit many-body polarization terms in the potential energy functions.

An explicit parameterization of f_l and f_m in Eq. (1) has been attempted for water using ab initio quantum mechanical calculations [29], [30], but it has been unsuccessful in simulations of liquid water. A potential parameterized from ab initio calculations is referred to as a nonempirical potential compared to the empirical potential parameterized from computer simulations of experimental properties. The major obstacle in parameterizing an ab initio potential is the choice of an appropriate set of analytical functions that allows for enough flexibility to describe the complex potential surface.

One way to incorporate many-body polarization effects in a potential energy function is to combine a nonempirical and an empirical potential. This can be constructed by explicitly incorporating many-body polarization terms. The polarization energy f_pol is added to the pairwise effective potential.

$$\text{V(r}{}_1\text{,…,r}{}_{\mathrm{N}}\text{)≈}\text{∑}\text{ij}\text{N}\text{f′}{}_{\mathrm{<mtext>eff</mtext>}}\text{(r}{}_{\mathrm{i}}\text{,r}{}_{\mathrm{j}}\text{)+}\text{∑}\text{i}\text{N}\text{f}{}_{\mathrm{<mtext>pol</mtext>}}\text{(r}{}_1\text{,…,r}{}_{\mathrm{N}}\text{).}$$

The effective potential f′_eff differs from f_eff in Eq. (3). In this manner, the effective pair potential energy f′_eff is derived to reproduce physical properties in computer simulations with the polarization energy f_pol which is parameterized with data from ab initio calculations.

Polarizable force fields can be divided into two main categories [51]. One is the dipole-polarizable models [8], [31], [34], [35], [36], [37], [40], [41], [43], [45], [50] in which the induction effects are described by polarizable dipoles centered on atomic sites. The other are the fluctuating-charge models [32], [33], [38], [39], [42], [46], [48], [49], [51] in which atomic charges are changed according to their environments.

Most of the approaches that include polarization and many-body effects explicitly have been implemented to devise more sophisticated potentials [8], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [45], [46], [47], [48], [49], [50], [51] because polarization and many-body effects are difficult to represent by simple analytical functions. The sophisticated potentials have been abandoned in favor of simple effective potentials because they do not lead to better results compared to the computational time disadvantage arising from the additional terms.

Since polarization is the largest contribution to the many-body effect, inclusion of polarization is essential for the development of a polarizable force field for water molecules. The major obstacle in parameterizing the many-body polarization effect is the choice of an appropriate set of analytical functions that allows for enough flexibility to describe the complex potential surface.

Polarization is not an atomic but rather a group property. This means that polarization should be described by a molecular (many-body) interaction rather than by single-atom pair interactions. The magnitude of the polarization can vary with the interacting distance and the orientation of the interacting molecules. This characteristic feature of polarization leads to a complicated functional form for its description. One example of this type is the NCC model [34] which combines a two-molecule potential with a polarization term. The functional form of the NCC model demonstrates the complexity of some empirical models. Because of the lack of analytical functions to describe the complexity of polarization, efforts to include polarization have still not been successful.

In this paper, we suggest a new force field for water, including the polarization effect. The new potential is a combination of empirical and non empirical potentials. The TIP4P model is used for the empirical potential part, which represents Coulombic and LJ interactions. For the nonempirical potential part, an artificial neural network (ANN) is introduced to reproduce the complex many-body interaction-energy surface from ab initio calculations. Possible double counting of the polarization effect, because of its implicit inclusion in the TIP4P model, is rationalized at the end of Section 2.1.

ANNs have been used to describe the whole potential energy surface (PES) of water [44], [47] rather than only the polarization part, the latter being the part computed in the present paper. When we tried to use these methods [44], [47] in MC simulations, we could not reproduce the physical properties of liquid water. Apparently, the source of the problem was the failure of these methods [44], [47] to avoid configurations with close contacts of two water molecules in training the weight matrix. Such close contacts were avoided here by using TIP4P to represent the electrostatic and nonbonded interactions, instead of using an ANN to describe the whole PES. Our procedure uses the repulsive energy of TIP4P to avoid close contacts, and an ANN only for the polarization contribution.

ANNs are mathematical models and algorithms that have been designed to mimic the information processing and knowledge acquisition methods of the human brain. ANNs are nonlinear regression methods to predict nonlinear behavior of natural phenomena. In this work, ANNs, instead of analytical functions, were introduced to include the influence of the higher-order terms of the water dimer PES. For training the weight matrix in the ANN, use was made of the back-propagation of error (BPE) algorithm [52]; we designate the ANN trained with BPE as BPNN. BPE is not the name of a specific neural network architecture, but the name of a learning algorithm, a strategy for the correction of weights. The procedure uses a supervised learning method, i.e. it requires the answers (targets) to the inputs to be known in advance in order to train the neural network. Applications of ANNs in the field of chemistry, and the computational details of BPNN, have been summarized by Zupan and Gasteiger [52].

We have determined some structural and physical properties of liquid water by Metropolis Monte Carlo (MC) simulation [53], [54], [55], using the new water model including many-body polarization effects with BPNN. MC simulations were carried out in the NPT ensemble, and the results were compared with experimental data and with simulation results from other potential energy functions. The simulation conditions and details are presented in Section 2.

Applications of this approach for the development of other force fields, besides water–water interactions, and computational time issues are addressed in Section 3.