Modeling Charge Penetration Effects in Water-Water Interactions

Abstract

This report introduces Gaussian electrostatic models (GEMs) to account for charge penetration effects in water-water interactions, allowing electrostatic interactions to be accurately described. Three different Gaussian electrostatic models, GEM-3S, GEM-5S, and GEM-6S are designed with s-type Gaussian functions. The coefficients and exponents of the Gaussian functions are optimized using the electrostatic potential (ESP) fitting procedure based on that of the MP2/aug-cc-pVTZ method. The electrostatic energies of ten different water dimers that were calculated with GEM-6S agree well with the results of symmetry-adapted perturbation theory (SAPT), indicating that this designed model can be effectively applied to future water models.

Introduction

Over the last few years, a large number of water models have been continuously introduced due to their fundamental importance for implementing Monte Carlo and molecular dynamics simulations in chemistry and biology.1-6 The increasing accuracy with respect to the treatment of electrostatic interactions has emerged as a central point in the development of these potentials. Although several popular water models, such as TIP4P,1 DC (Dang-Chang),2 TTM (Thole-Type Model),3 and DPP (Distributed Point Polarizable),5 successfully described the relative energies of water clusters, none of them considered charge penetration effects7 because these models used a classical point-charge representation of electrostatics.

The intermolecular electrostatic contribution in the short range is considered by the overlap of the charge densities when two molecules are in close proximity. The electrostatic interaction between these two species is no longer well represented by methods such as Stone’s distributed multipolar analysis (DMA)7 as the nuclei on one molecule are no longer shielded by its own electron density, thus experience a greater attraction for the electron density associated with the other species. The energy resulting from this increased attraction is referred to as charge penetration.

Recently Kumar et al. introduced a second generation distributed point polarizable water model (DPP2) by applying explicit terms for charge penetration.8 This model has been successfully improved for describing interaction energies in small and large water clusters by applying the simple damping scheme for point-charge electrostatics. Introducing a screening function as a pre-factor to the Coulomb potential has been a popular method used to address charge penetration.9-12 This method does have the advantage of including the short-range electrostatic penetration effect at a very low cost, but it still does not account for all the nuclear-electron interactions. As an alternative method, the Gaussian electrostatic model (GEM) aims to fit electron densities to Gaussian auxiliary basis sets (ABSs) centered on specific molecular sites, thus can accurately generate the Coulomb interaction energies and explicitly considers nuclear charge and electron densities. This method has been developed and studied in depth by Cisneros and co-workers.13-17 As a recent work, they introduced distributed multipoles based on the GEM (GEMDM) applied in the AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Applications) force field to improve the treatment of electrostatics.17 The ab initio calculated electron density can be fit with a number of basis functions of differing angular momentum, but ABSs are usually restricted to s-type Gaussian functions for the computational efficiency.

In the present study, Gaussian electrostatics models on specific molecular sites have been introduced to describe the charge penetration effects. Three different GEMs have been suggested based on the DPP and SAPT-5S water models:18 GEM-3S has one s-Gaussian type function (s-GTF) at each of three sites, GEM-5S has one s-GTF located at each of five sites, and GEM-6S has two s-GTFs at each of three sites. The coefficients and exponents of the Gaussian functions of the ABSs are fitted to theMP2/aug-ccpVTZ electrostatic potential. The validities of the three models have been examined by calculating the intermolecular Coulomb interactions for several stationary points on the water dimer potential energy surface. These results are compared with those of symmetry-adapted perturbation theory (SAPT), absolutely localized molecular orbitals energy decomposition analysis.8,26

 

Methodology

Models. Figure 1 compares the DPP water model and the three different GEMs considered in this work. GEM-3S locates two s-GTFs at the hydrogens and one at the M site, which is displaced along the HOH bisector. Cisneros et al. had tested a similar three-site model, where the GTFs were located at the oxygen and two hydrogen sites.16 This resulted in a water monomer dipole moment of about 1.9-2.0 Debye, which is much higher than the experimental result (1.855 Debye). Thus, GEM-3S suggested in this study can be applicable as a better descriptor of the dipole moment. GEM-5S employs five sites, one on each atom and two D sites located above and below the HOH plane. The geometry of GEM-5S is adapted from the SAPT-5S model,18 which is a rigid model based on the vibrational averaged geometry. The positions of the D sites in GEM-5S have been optimized to better fit the electrostatic energy. GEM-6S places two s-GTFs at each of the three sites of GEM-3S for a total of six s-GTFs. Nuclear point charges of +6 and +1 are placed on M and each H, respectively, for both GEM-3S and GEM-6S, whereas those of GEM-5S are located on O and each H.

Figure 1.Gaussian electrostatic models: GEM-3S, GEM-5S, and GEM-6S have 3, 5, and 6 s-types.

The electrostatic interaction energy for the classical model using the point-charge (PC) representation is given by

where QAi and QBj represent the point charges at i and j sites of two different monomers A and B, respectively. For the GEMs, a more detailed description of the intermolecular Coulomb energy is required and nucleus-nucleus (N−N), nucleus-electron (N−e), and electron-electron (e−e) contributions must be considered.

The following expression gives the intermolecular Coulomb interactions in the GEM model.

The first term represents the nucleus-nucleus interaction, the second and third terms are the nucleus−electron contributions, and the last term is the electron−electron interaction. ZAi and ZBj represent the nuclear charges of atom i on molecule A and atom j on molecule B, respectively. Only valence nucleus charges are considered, i.e., +6 and +1 are used for O and H, respectively. ρAand ρB represent the Gaussian electron densities of molecules A and B, respectively.

Fitting Procedure. Gaussian electron densities can be found through fits to several molecular properties: density, electrostatic potential (ESP), and electric field. In this work, the ESP fitting procedure has been applied based on the reference ESP of the MP2/aug-cc-pVTZ method. The densities in this calculation are based on Z-vector27 and hence yield multipole moments which are correct analytical derivative energy. The least-squares method is used to optimize the coefficients and exponents of GTFs in addition to locating the D sites by minimizing the following fitting function:

where yα(ri) corresponds to the ab initio molecular property of interest at point i, yα(ri, ck) is the same property evaluated with the kth ABS element at the same point on the grid, and α denotes the molecular properties of interest. wα is the relative weight for property α, and W(r) is the weighting function for the point on the grid as suggested by Hu et al.19 and defined as:

in which ρ(r) is the predefined electron density, and ρref is a reference electron density that was chosen in conjunction with σ to ensure the weighting function behaves as expected. Eq. (4) is minimized by using both the genetic algorithm (GA) and quasi-Newton optimization methods, which employ analytical gradients with respect to the coefficients and exponents of GTFs as well as the position of the D sites. Due to the symmetry of the water monomer, the fitting points span a quarter of the space. As suggested in Reference 16, the minimum cutoff values of 1.5 and 1.7 bohr were used for hydrogen and oxygen, respectively, and a maximum cutoff value of 10 bohr was used with the grid point 0.1 Å. This brings the total number of data points to 157412.

 

Results

The D sites of GEM-5S, which are located beyond the HOH plane, were optimized and fitted to ABS, so that the positions of the D sites are different from original locations based on the SAPT-5S model. Figure 2 compares the geometry of the SAPT-5S model with that of GEM-5S. The positions of the D sites for GEM-5S have been optimized, and it resulted in bigger inner angles than those of SAPT-5S. Figure 2 also shows the geometries of GEM-3S and GEM-6S, which include the positions of the M sites. The M sites have been placed at 0.10 Å away from oxygen along the HOH bisector by the trial-and-error method.

Figure 2.The original SAPT-5S model and the optimal geometry of GEM-5S, and the geometries of GEM-3S and GEM-6S.

The optimized coefficients and exponents of the Gaussian type functions (GTFs) of GEM-3S, GEM-5S, and GEM-6S are shown in Table 1. The dipole moments of the monomer calculated with the GTFs in each model are also indicated in this table. The dipole moment of GEM-5S is very close to the experimental value, 1.855 Debye,20 whereas those of GEM-3S and GEM-6S are slightly underestimated and overestimated, respectively.

Table 1.Optimized coefficients (a.u.) and exponents of GTFs for the electron density and the dipole moments (Debye) of GEM-3S, GEM-5S, and GEM-6S, respectively

Electrostatic Potential. To assess the fitting performance of the three different models for representing the electrostatic potential of each water monomer, I report the difference between the electrostatic potential associated with each model and that calculated at the MP2/aug-cc-pVTZ level of theory as a function of distance. Figure 3 shows the difference contours at the ± 0.5 kcal/mol level at points in the plane of the water monomer. The ± 0.5 kcal/mol contours were chosen for display because they represent a sizable portion (10%) of the net binding energy of the water dimer. Also I added the test charge (extra oxygen) in the figure to indicate the location of the second water molecule in the optimized geometry of the water dimer.

Figure 3.Electrostatic potential difference plots for the water monomer. In each case the plot gives the difference of the electrostatic potentials calculated using the model and the MP2/aug-ccpVTZ method. An extra oxygen atom is located in the optimized geometry of the water dimer.

The electrostatic potential differences for the DPP model were greater than 0.5 kcal/mol when the test charge is brought within near 3.5 Å of the center. The difference between the contours of the DPP and DPP2 models was not found to be within ± 0.5 kcal/mol. The 0.5 kcal/mol contour lines of the ESP differences in GEM-5S were located near ± 3.0 Å, whereas those in GEM-3S and GEM-6S were found to be near ± 4.5 Å. Water molecules are usually bound with other waters at near ROO = 3.0 Å, so this figure indicates GEM-5S shows accurate electrostatics potentials within ± 0.5 kcal/mol errors when the test charge is located in the region of longer than 3.0 Å whereas GEM-3S and GEM-6S indicate this distance as 4.5 Å. Compared to the DPP and DPP2 water models, only GEM-5S exhibited better performance. GEM-3S and GEM-6S performed worse than the DPP water model during the analysis of the electrostatic potential difference of the water monomer. However, this contour analysis represents the information of only HOH plane.

Water Dimer. Figure 4 displays the electrostatic interaction energy of the water dimer as a function of the oxygenoxygen separation between these two water monomers (ROO) for each considered model. This plot was generated while the monomers were kept rigid by fixing the OH bond lengths and HOH angles fixed at the experimentally determined values for the gas phase monomer. The flap angles that determine the orientation of these two monomers with respect to the OO axis were frozen at the values they have for the MP2/aug-cc-pVTZ optimized structure. More attractive electrostatic interactions are shown with more negative energies in this figure.

Figure 4.Electrostatic interaction energy for the water dimer, as a function of O−O distance, assuming rigid monomers and fixed flap angles.

The electrostatic interaction energies of the water dimer are compared with the values obtained from symmetryadapted perturbation theory (SAPT)21,22 with the aug-ccpVTZ basis set as a reference. The energy of the DPP water model less attractive at short OO distances. Most classical water models, such as TTM2-R, DC, and AMOEBA, give nearly the same electrostatic interaction energies as the DPP water model.5 These classical models are limited in their ability to express electrostatic energies at short distances since they do not account for charge penetration effects. In the case of the DPP2 model, the electrostatic energies are perfectly matched to the SAPT results because the parameter of the screening function was obtained by a least-square fitting of the electrostatic energies from the SAPT calculations.

The electrostatic energy of GEM-3S is more attractive in the region of ROO ≤ 2.8 Å, but less attractive in the region of ROO > 2.8 Å. GEM-5S and GEM-6S provide better electrostatic energy calculations of the dimer than GEM-3S because these models used five and six Gaussian type functions to describe the electron density, respectively. In the region of ROO > 2.7 Å, the electrostatic energies of GEM-5S are almost the same as those of the SAPT calculation.

The Ten Water Dimers. As a more general test for the water dimer, ten stationary points of the water dimer’s potential energy surface (PES) were considered. Figure 5 shows the geometries of the ten stationary points originally found by Smith et al.23 (the so-called Smith dimer set). In this work, they have been characterized by employing secondorder Moller-Pleset perturbation theory (MP2), and the OH bond lengths and HOH angles were fixed so that they conform to the rigid geometries of the DPP model.

Figure 5.The Smith dimer sets.

The reference electrostatic energies of the ten stationary points of the dimer were determined from the SAPT decomposition results.8,26 Figure 6 shows the errors in intermolecular Coulomb energy for the ten water dimers using the three investigated models, which are compared with those of the DPP and DPP2 water models. The DPP water model has smaller electrostatics energies than the SAPT results because of the lack of charge penetration. The electrostatic energies of GEM-3S, GEM-5S, GEM-6S, and DPP2 models of dimers 1−3 were almost the same as the SAPT results. However, for the stationary points of dimers 4−10, the electrostatic energies of all models greatly deviated from the SAPT results. For dimers 4−10, GEM-5S overestimated the electrostatic energies and gave values with inaccuracies comparable to the DPP water model. GEM-3S produced more accurate results than the DPP water model for all ten dimers, though the errors of the stationary points 7 and 9 were greater than 1 kcal/mol. Compared to the DPP2 model, GEM-3S give better performances for stationary points 5, 6, 9, and 10. GEM-6S reproduced the electrostatic energies for the ten stationary points quite well. It outperformed the other considered models of this study. Quite notably, GEM-6S gave significant improvements for dimers 7, 9, and 10, and the errors for all ten stationary points were less than 1 kcal/mol using this model.

Figure 6.Errors in the electrostatic interaction energies for the ten dimers.

 

Conclusion

Three different GEM models have been designed and implemented to reproduce ESP using ABSs. The ABS sites of GEM-3S, GEM-5S, and GEM-6S were limited to three, five, and six sites with s-type Gaussian functions, respectively, The electrostatic energies of the ten dimers were calculated to determine how well these models reproduce charge penetration.

GEM-5S analysis caused the electrostatic energies to be more attractive for seven out of the ten dimers, though it did exhibit the best electrostatic potential among the three tested models based on the difference of the electrostatic potential in the water monomer. The contour diagram of the water monomer only provides accurate information for the electrostatic energies of the HOH plane. Therefore, only three dimers, which have the interactions between the HOH planes, were well described with this model. In contrast to GEM-5S, GEM-6S performed more poorly at describing the electrostatic potential in the contour diagram of the water monomer than the DPP water model. However, the electrostatic energies of the ten dimers calculated with the GEM-6S model were in good agreement with the SAPT results, which indicated better performance than those of even the DPP2 model.

Unexpectedly, GEM-5S did not perform well in the ten water dimers, but GEM-6S agreed well with the SAPT results, which indicates that the M site is a suitable position to put auxiliary basis sets. GEM-3S can also be applied to water models with less computational cost. It is about twice faster than GEM-6S for calculating the electrostatic energy of a certain point.