J. Phys. Chem. B 2000, 104, 2117-2122
2117
Charge Flow between Ions and a Dielectric Continuum. 2. Variational Method for Distributing Charge into the Dielectric Valentin Gogonea and Kenneth M. Merz, Jr.* Department of Chemistry, 152 DaVey Laboratory, The PennsylVania State UniVersity, UniVersity Park, PennsylVania 16802 ReceiVed: NoVember 5, 1999; In Final Form: December 15, 1999
The paper presents a variational method to distribute the charge transferred (CT) between an ion and a dielectric continuum onto a boundary element (BE) grid used to represent the dielectric continuum media. The CT charge is transferred into the dielectric using a methodology that combines a semiempirical effective Hamiltonian with a modified Poisson-Boltzmann equation, which includes CT in the form of a surface charge density positioned at the dielectric interface. The variational method eliminates the bias of our earlier method used to distribute CT charge into the dielectric, i.e., the need to specify the solute atoms used in constructing the BE grid. The surface “pattern” of the magnitudes of the CT charges placed onto the BE grid is obtained by minimizing an energy functional defined as the sum of the CT charge self-interaction energy and the energies due to CT charge interaction with the reaction field (RF) and solute partial atomic charges. The total amount of charge, distributed on the BE grid, is constrained to the amount of charge transferred from/to the solute HOMO/LUMO molecular orbital, by using the method of Lagrange multipliers.
Introduction The recent progress in developing linear-scaling algorithms for matrix diagonalization1-5 has facilitated the study of charge transfer (CT) interactions between biomolecules in either the gas phase or solution. For example, Nadig et al.6 demonstrated that about 2 eu are transferred between the major cold shock protein (CspA) and the surrounding water. This shows that interfacial CT interactions could be important for the stability and reactivity of biomolecular systems.7-9 In paper I,10 we showed that the flow of charge between ions and a dielectric continuum medium can be modeled by mixing a quantum mechanical Hamiltonian with a modified PoissonBoltzmann equation. The charge transferred into the dielectric is represented as a surface charge density located at the dielectric interface, which interacts with the surface charge density due to the reaction field generated by solvent polarization. Our previous calculations (paper I) indicated that the magnitude of CT obtained with the dielectric continuum method is in good agreement with the full quantum mechanical calculations on equilibrated molecular dynamic configurations of organic ions embedded into a box of water molecules. In paper I, we used the boundary element method (BEM) to implement the CT modified Poisson-Boltzmann equation. The CT charge was distributed only on that part of the boundary element (BE) grid which “covered” solute atoms believed to be involved in CT. For example, in the case of the acetate anion we distribute the CT charge on the grid covering only the carboxylate oxygen atoms. Here we present a variational method (VCT) for distributing the CT charge on the BE grid, which eliminates the bias found in the uniformly spread-on-the-grid charge method (UCT) used in paper I. This new approach for distributing CT onto BE grid is the topic of this paper. Below we give a brief description of the formalism used in the CT dielectric continuum methodology, and introduce the energy functional for CT charge interaction used in the VCT
method. Then we investigate the dependence of the amount of CT charge on the BE grid resolution, compare the results obtained with the UCT and VCT methods, and comment on the efficiency of the two methods. Method CT-Modified Poisson-Boltzmann Equation. In paper I,10 we have shown that the CT between charged molecules and a dielectric continuum can be modeled using a modification of the Poisson-Boltzmann equation which includes the CT charge in the form of a surface charge density, σCT, located at the dielectric interface. The CT charge equals a fraction of one spincoupled electron pair removed from the highest occupied molecular orbital (HOMO, for solute to solvent CT), or placed in the LUMO orbital. The modified Poisson-Boltzmann equation has the form
∇‚[(r)∇φ(r)] - κ(r)2 φ(r) ) -4π[F(r) + δ(r - r′)σCT(r)], r′ ∈ S (1) and obeys the following boundary conditions (at the dielectric interface) for the dielectric displacement11 (eq 2) and the electric field (eq 3):
(oEo - iEi)‚n ) 4πσCT (Eo - Ei)‚n ) 4π
(
σCT + σRF o
(2)
)
(3)
φ(r) is the electrostatic potential, Eo and Ei are the electric fields at the dielectric interface, o and i are the solvent and solute dielectric constants, κ2 is the Debye-Hu¨ckel parameter, F(r) is the solute charge density, and n is the normal vector at the dielectric interface S. The delta function δ(r - r′) in eq 1 confines σCT to the dielectric interface (vide infra).
10.1021/jp9939064 CCC: $19.00 © 2000 American Chemical Society Published on Web 02/16/2000
2118 J. Phys. Chem. B, Vol. 104, No. 9, 2000
Gogonea and Merz
The solvation free energy, Gsol, will also contain a term which accounts for all energy changes due to transferred charge, GCT:
Gsol ) GRF + GCT + Gwfd + Gnp
(4)
Gnp is the hydrophobic (cavity formation and dispersionrepulsion) term,12 and Gwfd so-called wave function distortion term,13,14 CT Variational Spread-on-the-Grid Charge Method for Boundary Element Grid. As we explained in paper I, we decided to use the boundary element method (BEM)15-17 to model charge flow between ions and dielectrics because CT is restricted mainly to the first solvation layer.18 For completeness we give here the new BEM equations, which include CT:
(
2πI -
)
[
i - o i - o 1 Kσ j RF ) T - (2πI - K)σ j CT i + o i + o o
]
(5)
σ j RF and σ j CT are vector quantities which contain the set of reaction field (RF) and CT charge densities, respectively, and T is a vector whose elements (eq 6, Ti ) (Es‚n)i) are normal components of the electric field generated by solute charge density F on boundary elements:
(Es‚n)i )
F(r)(ri - r)‚n
∫V dr
(6)
i|ri - r|3
Kij terms depend only on the shape of the dielectric interface and have the following form:
Kij )
(r - r )‚n
∫S |ri - rj |3 dSj
(7)
j
i
j
The amount of charge transferred (qCT) between solute and the dielectric continuum is found by minimizing the total semiempirical energy of the solute in solution. The energy is obtained as the solution of the Schro¨dinger equation which uses, in this particular case, an effective semiempirical Hamiltonian including energy terms due to mutual interaction of RF18 and CT charges, and with the solute wave function. The expectation value of this Hamiltonian is
E ) 〈Ψ|H|Ψ〉 + Enucl +
Nsc Nat Z
∑i ∑A
i A(qRF
+ qiCT)
|ri - rA|
+
∫VERF(r)‚ECT(r) dr
(8)
The last two terms in eq 8 are the electrostatic interaction of solute cores with RF (qRF) and CT (qCT) charges, and the interaction between RF and CT charge densities, respectively. The pattern of CT charge density magnitudes distributed onto the BE grid is determined by minimizing the functional E (eq 9), assuming that a pattern of RF charge densities already exist on the BE grid:
E)
∫VECT‚(ECT + ERF + Es) dr + λ(∫sσCT ds - qCT)
(9)
ECT and ERF are the electric fields generated by CT and RF surface charge densities and Es is the electric field due to solute wave function. λ is the Lagrange multiplier.19 The CT charge densities on BE's are obtained as the solution of the following system of linear equations:
[(8π2 + ai)I + PCT]σ j CT ) (4πI + PRF)σ j RF - Psq - λIs (10) s and q are vector quantities whose elements are the areas (si) of BE's and solute partial atomic charges (qi), respectively, and ai is a geometric factor which accounts for the curvature of BE i
ai )
ds ds
∫s |rj -j rkk| i
(11)
where the j and k indices run over surface area subelements belonging to BE i. The elements of the Ps, PRF, and PCT matrices are as follows:
(Ps)ij )
si |ri - rj|
(PRF)ij ) (Ps)ijsj (PRF)ij 2
(PCT)ij ) (1 - δij)
(12)
The algorithm for variationally calculating the CT charge densities is as follows: (1) The RF charge densities are obtained from the BEM equations20 by assuming no CT. (2) Then the CT charge densities are obtained my minimizing the energy functional in eq 9, i.e., by solving the system of linear equations given in eq 10, under the constraint that the total CT surface charge equals the amount of charge “taken” from the solute HOMO or “placed” in the LUMO orbital. (3) The expectation value (the SCF energy) of the effective Hamiltonian, which includes terms due to the interaction of the RF and CT surface charges with the solute wave function, is calculated after removing from the solute HOMO or adding into the LUMO the amount of charge created in the dielectric. Because of ion discharging, the density matrix used as a “guess” in the SCF calculation contains a larger (anion) or smaller (cation) number of electrons in the full charged ion. Consequently, this step is iterated until convergence of the SCF energy is achieved. Results and Discussion We carried out CT-dielectric continuum calculations with the semiempirical AM1/PM3 Hamiltonians and charges were obtained using the Mulliken, CM1,21 and CM222 approaches. The test set consisted of 36 ions (Table 1) and 24 neutral molecules (Table 2). The results for the shift in total semiempirical energy, solvation free energy, and the amount of CT are given in Tables 1 and 2. The dielectric continuum was assigned a dielectric constant of 80 (highly polar solvent) and the solute dielectric constant was set to 1. The BEM grid was generated and triangulated with Connolly’s dot surface program.23 A finer subgrid was also generated to take into account the curvature dependence of the surface charge density.20 Because the accuracy of calculating physical quantities (e.g., solvation free energy) with grid-based methods usually depends on grid resolution, we first investigated the change in the amount of CT as a function of the BE density. Table 3 shows that accurate quantities can be obtained with a grid resolution higher than 4 BE/Å2. Beyond a resolution of 4 BE/Å2 the calculations become prohibitively expensive as can be seen from the trends in CPU time given in Table 4. By using a variational method to determine the CT charge density on the BE, we were expecting that the pattern of the
Charge Flow between Ions and a Dielectric Continuum
J. Phys. Chem. B, Vol. 104, No. 9, 2000 2119
TABLE 1: Computational Results for Semiempirical AM1/PM3 Calculations of CT between Small Ions and a Polar Continuum Dielectric Media. Comparison between the Uniformly (UCT) and Variational (VCT) Spread-on-the-Grid Charge Methods AM1 compda CNMullikenf UCTg VCTh CM1i UCT VCT CM2j UCT VCT NH2 Mulliken UCT VCT CM1 UCT VCT CM2 UCT VCT HOMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT N3Mulliken UCT VCT CM1 UCT VCT CM2 UCT VCT CH3OMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT C6H5OMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT HSMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT CH3SMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT CH3CH2SMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT C6H5SMulliken UCT VCT CM1 UCT VCT CM2 UCT VCT CH3COOMulliken UCT VCT CM1 UCT
∆Eb
GCT/ Gexpc
-0.210 -0.370 -0.191 -0.755 -0.187 -0.883
PM3
qCTd
GRF/ Gexpe
1.18 1.20 1.27 1.35 1.31 1.40
-0.10 -0.10 -0.09 -0.10 -0.09 -0.10
-1.674 -1.577 -5.354 -4.746 -3.604 -2.386
1.15 1.10 1.76 1.67 1.48 1.28
-1.661 -1.370 -2.353 -2.230 -2.369 -2.243
AM1
∆E
GCT/ Gexp
qCT
1.09 1.09 1.18 1.18 1.21 1.21
-0.085 -0.261 -0.074 -0.678 -0.073 -0.698
1.17 1.19 1.29 1.39 1.29 1.40
-0.06 -0.06 -0.06 -0.09 -0.06 -0.09
1.10 CM2 1.10 1.20 H3O+ 1.20 Mulliken 1.20 1.20 CM1
-0.25 -0.25 -0.37 -0.38 -0.37 -0.29
1.02 1.02 1.33 1.33 1.10 1.10
-1.255 -1.235 -3.767 -3.465 -3.741 -3.446
1.13 1.10 1.59 1.54 1.59 1.53
-0.22 -0.23 -0.38 -0.39 -0.38 -0.39
1.02 CM2 1.02 1.14 CH3OH2+ 1.14 Mulliken 1.14 1.14 CM1
1.09 1.03 1.22 1.21 1.22 1.21
-0.24 -0.22 -0.29 -0.32 -0.29 -0.32
0.91 0.91 0.94 0.94 0.94 0.94
-1.513 -1.186 -2.255 -2.126 -2.239 -2.112
1.10 1.04 1.24 1.22 1.24 1.22
-0.23 -0.20 -0.28 -0.27 -0.28 -0.28
0.91 CM2 0.91 0.95 NH4+ 0.95 Mulliken 0.95 0.95 CM1
-0.319 -0.370 -0.971 -1.914 -0.319 -0.370
1.09 1.11 1.38 1.58 1.09 1.11
-0.12 -0.15 -0.14 -0.27 -0.12 -0.15
1.04 1.04 1.27 1.27 1.04 1.04
-0.456 -0.604 -0.456 -0.604 -0.456 -0.604
1.21 1.26 1.21 1.26 1.21 1.26
-0.14 -0.18 -0.14 -0.18 -0.14 -0.18
1.06 CM2 1.06 1.06 CH3NH3+ 1.06 Mulliken 1.06 1.06 CM1
-0.985 -0.431 -1.187 -0.523 -1.168 -0.507
1.04 0.92 1.10 0.96 1.09 0.95
-0.21 -0.07 -0.23 -0.08 -0.23 -0.07
0.84 0.84 0.86 0.86 0.85 0.85
-1.043 -0.455 -1.320 -0.568 -1.290 -0.550
1.06 0.93 1.14 0.98 1.13 0.97
-0.21 -0.07 -0.23 -0.07 -0.23 -0.07
0.85 CM2 0.85 0.88 (CH3)2NH2+ 0.88 Mulliken 0.87 0.87 CM1
-0.441 -0.300 -0.591 -0.360 -0.691 -0.385
1.07 1.01 1.12 1.05 1.15 1.04
-0.15 -0.08 -0.17 -0.08 -0.19 -0.08
0.88 0.88 0.90 0.90 0.89 0.89
-0.375 -0.295 -0.592 -0.416 -0.607 -0.421
1.06 1.01 1.15 1.08 1.16 1.08
-0.13 -0.07 -0.17 -0.08 -0.17 -0.08
0.87 CM2 0.87 0.90 (CH3)3NH+ 0.90 Mulliken 0.90 0.90 CM1
-0.370 -0.380 -0.341 -0.364 -0.433 -0.419
1.23 1.23 1.22 1.21 1.26 1.25
-0.13 -0.13 -0.13 -0.13 -0.14 -0.13
1.12 1.12 1.11 1.11 1.13 1.13
-0.247 -0.340 -0.275 -0.353 -0.380 -0.392
1.18 1.19 1.19 1.19 1.23 1.22
-0.11 -0.12 -0.12 -0.13 -0.13 -0.13
1.10 CM2 1.10 1.10 C6H5NH3+ 1.10 Mulliken 1.12 1.12 CM1
-0.726 -0.597 -0.703 -0.564 -1.046 -0.770
1.22 1.17 1.22 1.15 1.34 1.24
-0.17 -0.10 -0.17 -0.10 -0.20 -0.11
1.04 1.04 1.03 1.03 1.09 1.09
-0.908 -0.601 -0.974 -0.618 -1.152 -0.709
1.24 1.14 1.27 1.15 1.33 1.19
-0.18 -0.10 -0.19 -0.10 -0.21 -0.10
1.02 CM2 1.02 1.03 imidazole (H+) 1.03 Mulliken 1.06 1.06 CM1
-0.680 -0.573 -0.660 -0.544 -0.967 -0.743
1.27 1.22 1.26 1.20 1.39 1.31
-0.17 -0.11 -0.17 -0.10 -0.20 -0.11
1.08 1.08 1.06 1.06 1.14 1.14
-0.868 -0.573 -0.934 -0.589 -1.133 -0.690
1.25 1.14 1.28 1.15 1.35 1.21
-0.18 -0.09 -0.19 -0.09 -0.21 -0.10
-0.186 -0.324 -0.213 -0.336 -0.474 -0.478
1.23 1.26 1.25 1.27 1.37 1.35
-0.09 -0.09 -0.10 -0.10 -0.14 -0.10
1.07 1.07 1.07 1.07 1.12 1.12
-0.192 -0.305 -0.252 -0.340 -0.371 -0.419
1.27 1.28 1.31 1.31 1.38 1.37
-0.09 -0.09 -0.11 -0.10 -0.13 -0.10
1.02 CM2 1.02 1.02 pyridinium (H+) 1.02 Mulliken UCT -1.176 0.90 1.06 VCT -1.144 0.73 1.06 CM1 UCT -1.353 0.96 VCT -1.339 0.72 1.09 CM2 UCT -1.302 0.92 1.09 VCT -1.286 0.75 + 1.10 CH3SH2 1.10 Mulliken UCT -0.708 0.87 1.14 VCT -1.177 0.81 1.14 CM1 UCT -1.145 0.82 VCT -0.779 0.73 1.01 CM2 UCT -0.638 0.83 1.01 VCT -0.904 0.75 1.04
-0.247 1.19 -0.12 1.01 -0.147 1.15 -0.09 -0.248 1.15 -0.09 1.01 -0.194 1.15 -0.09 -0.295 1.22 -0.13 1.03 -0.210 1.22 -0.11
GRF/ Gexp
compda
∆Eb
GCT/ Gexpc
PM3
qCTd
GRF/ Gexpe
∆E
GCT/ Gexp
qCT
GRF/ Gexp
VCT -0.301 1.19 -0.10 1.03 -0.277 1.21 -0.10 1.04 UCT -0.311 1.24 -0.14 1.03 -0.216 1.22 -0.11 1.04 VCT -0.334 1.21 -0.11 1.03 -0.286 1.22 -0.10 1.04 UCT VCT UCT VCT UCT VCT
-1.889 -1.868 -2.211 -2.406 -2.206 -2.399
0.78 0.78 0.82 0.85 0.82 0.85
0.29 0.28 0.31 0.31 0.31 0.31
0.90 0.90 0.95 0.95 0.95 0.95
-2.017 -2.010 -2.327 -2.473 -2.308 -2.421
0.76 0.76 0.80 0.82 0.80 0.80
0.30 0.29 0.32 0.32 0.32 0.34
0.90 0.90 0.94 0.94 0.94 0.94
UCT VCT UCT VCT UCT VCT
-1.321 -1.102 -1.574 -1.422 -1.592 -1.404
0.73 0.70 0.77 0.76 0.77 0.76
0.25 0.23 0.27 0.24 0.27 0.23
0.88 0.88 0.91 0.91 0.91 0.91
-1.662 -1.120 -1.943 -1.357 -1.918 -1.331
0.73 0.69 0.77 0.76 0.76 0.75
0.28 0.23 0.30 0.22 0.30 0.22
0.90 0.90 0.93 0.93 0.92 0.92
UCT VCT UCT VCT UCT VCT
-1.175 -1.225 -1.270 -1.198 -1.366 -1.270
1.02 1.02 1.03 1.02 1.04 1.04
0.23 0.24 0.24 0.23 0.25 0.22
1.10 1.10 1.10 1.10 1.10 1.10
-1.169 -1.797 -1.392 -1.411 -1.389 -1.417
1.02 1.08 1.00 1.00 1.00 1.00
0.23 0.26 0.25 0.26 0.25 0.26
1.15 1.15 1.12 1.12 1.12 1.12
UCT VCT UCT VCT UCT VCT
-0.817 -0.636 -0.915 -0.765 -0.936 -0.874
0.99 0.98 1.02 1.04 1.03 1.06
0.20 0.15 0.21 0.15 0.21 0.14
1.05 1.05 1.08 1.08 1.08 1.08
-1.118 -0.907 -1.405 -1.034 -1.250 -0.909
1.01 1.06 1.06 1.13 1.03 1.07
0.24 0.14 0.26 0.13 0.24 0.14
1.11 1.11 1.14 1.14 1.11 1.11
UCT VCT UCT VCT UCT VCT
-0.668 -0.411 -0.906 -0.541 -0.946 -0.592
0.96 0.95 1.00 1.00 1.01 1.02
0.17 0.11 0.20 0.11 0.21 0.11
1.01 1.01 1.03 1.03 1.03 1.02
-1.075 -0.737 -1.696 -0.890 -1.351 -0.717
1.02 1.05 1.13 1.14 1.05 1.05
0.21 0.12 0.26 0.11 0.23 0.11
1.07 1.07 1.12 1.12 1.08 1.08
UCT VCT UCT VCT UCT VCT
-0.566 -0.291 -0.903 -0.348 -0.919 -0.348
0.96 0.87 1.01 0.91 1.01 0.91
0.14 0.13 0.17 0.12 0.18 0.12
0.98 0.98 1.00 1.00 1.00 1.00
-0.983 -0.471 -1.819 -0.554 -1.336 -0.450
1.01 0.94 1.15 0.99 1.06 0.93
0.18 0.12 0.23 0.12 0.20 0.13
1.03 1.03 1.07 1.07 1.03 1.03
UCT VCT UCT VCT UCT VCT
-1.729 -1.329 -1.955 -1.520 -1.904 -1.416
1.17 1.04 1.25 1.09 1.21 1.07
0.34 0.30 0.36 0.43 0.35 0.31
1.00 1.00 1.04 1.04 1.01 1.01
-2.090 -1.409 -2.452 -1.613 -2.193 -1.439
1.18 1.03 1.30 1.10 1.22 1.01
0.37 0.29 0.40 0.39 0.37 0.41
1.02 1.02 1.07 1.07 1.03 1.03
UCT VCT UCT VCT UCT VCT
-0.763 -1.041 -0.815 -1.445 -0.860 -1.471
0.87 0.84 0.94 0.94 0.93 0.93
0.23 0.27 0.24 0.33 0.25 0.35
0.98 0.98 1.03 1.02 0.99 1.03
-1.153 -0.803 -1.801 -1.398 -1.354 -1.229
0.85 0.78 1.00 0.90 0.92 0.80
0.27 0.24 0.34 0.31 0.30 0.35
0.99 0.99 1.04 1.04 1.01 1.01
0.23 0.31 0.24 0.39 0.25 0.34
1.00 1.00 1.03 1.03 1.01 1.01
-1.517 -1.038 -2.239 -1.275 -1.694 -1.185
0.92 0.73 1.02 0.77 0.93 0.70
0.25 0.28 0.31 0.31 0.27 0.32
1.03 1.03 1.05 1.05 1.01 1.01
0.17 0.24 0.21 0.21 0.16 0.22
1.01 1.01 0.95 0.95 0.97 0.97
-1.560 -1.406 -1.482 -1.302 -1.503 -1.307
0.75 0.65 0.73 0.61 0.74 0.62
0.27 0.29 0.26 0.29 0.26 0.29
0.96 0.96 0.93 0.93 0.93 0.93
a Single-point calculations on gas phase optimized geometry at the semiempirical level. The boundary elements (BE) grid density for the reaction field/CT calculations were set to 4 BE’s/Å2 (see text). A second grid with a density of 20 BE's/Å2 was used to account for the curvature dependence of the surface charge density. b The difference in the total semiempirical energy (eV) between the partially discharged and fully charged states of the ion in solution. c The ratio of the solvation free energy of the partially discharged ion to the experimental solvation free energy. d The amount of CT charge in eu. e The ratio of the solvation free energy of the fully charged ion to the experimental solvation free energy. f Mulliken (Coulson) charges. g Results obtained with the uniform charge grid method. h Results obtained with the variational method. i CM1 charges.21 j CM2 charges.22
2120 J. Phys. Chem. B, Vol. 104, No. 9, 2000
Gogonea and Merz
TABLE 2: Semiempirical AM1/PM3 Calculations of CT between Neutral Molecules and a Polar Continuum Dielectric Media: Results for the Variational Spread-on-the-Grid Charge Method AM1 compounda benzene Mullikenf CMg CM2h water Mulliken CM1 CM2 methanol Mulliken CM1 CM2 phenol Mulliken CM1 CM2 acetone Mulliken CM1 CM2 acetic acid Mulliken CM1 CM2 methylamine Mulliken CM1 CM2 dimethylamine Mulliken CM1 CM2 trimethylamine Mulliken CM1 CM2 3,5-dimethylpyridine Mulliken CM1 CM2 4-methylpyridine Mulliken CM1 CM2 acetamide Mulliken CM1 CM2 cis-N-methylacetamide Mulliken CM1 CM2 trans-N-methylacetamide Mulliken CM1 CM2
∆Eb
GCT/Gexpc
PM3 qCTd
GRF/Gexpe
∆E
GCT/Gexp
qCT
GRF/Gexp
-0.074 -0.074 -0.034
8.12 8.12 4.64
0.09 0.09 0.06
2.52 2.52 1.20
-0.054 -0.054 -0.060
4.64 4.64 5.09
0.08 0.08 0.08
0.59 0.59 0.77
0.0 -0.058 -0.060
0.15 1.93 1.99
0.0 0.04 0.04
0.13 1.04 1.06
0.0 -0.107 -0.098
0.10 2.32 2.21
0.0 0.05 0.05
0.07 1.04 1.00
0.0 -0.068 -0.068
0.14 1.42 1.40
0.0 0.02 0.02
0.10 0.79 0.77
0.0 -0.100 -0.092
0.01 1.46 1.38
0.0 0.02 0.02
-0.03 0.74 0.69
-0.171 -0.421 -0.346
1.67 3.18 2.57
0.12 0.17 0.14
0.44 0.95 0.81
-0.133 -0.387 -0.380
1.13 2.50 2.47
0.11 0.15 0.15
0.18 0.71 0.69
-0.001 -0.005 -0.005
1.00 1.54 1.81
0.01 0.01 0.01
0.64 0.99 1.23
-0.001 -0.009 -0.012
0.70 1.49 1.72
0.0 0.01 0.01
0.40 0.94 1.10
-0.028 -0.227 -0.205
1.46 2.24 2.28
0.09 0.09 0.09
0.51 0.93 1.00
-0.027 -0.263 -0.244
1.32 2.44 2.43
0.09 0.10 0.10
0.47 0.97 1.00
0.0 -0.124 -0.062
-0.11 2.96 2.01
0.0 -0.03 -0.02
-0.14 1.28 0.89
0.0 -0.028 -0.035
-0.41 1.33 1.37
0.0 -0.01 -0.01
-0.41 0.66 0.62
0.0 -0.015 -0.021
-0.15 1.11 0.82
0.0 0.01 0.01
-0.17 0.79 0.47
0.0 -0.052 -0.028
-0.45 0.72 0.66
0.0 0.01 0.01
-0.45 0.28 0.33
-0.001 -0.068 -0.032
-0.08 1.86 1.15
0.0 -0.02 -0.01
-0.21 0.64 0.23
0.001 -0.001 -0.016
-0.60 0.0 0.81
0.0 0.0 -0.01
-0.62 -0.24 0.05
-0.072 -0.101 -0.077
1.23 1.99 1.85
0.05 0.06 0.05
0.58 1.07 1.05
-0.045 -0.076 -0.094
0.46 1.27 1.81
0.04 0.06 0.07
0.11 0.62 0.96
-0.113 -0.150 -0.111
1.51 2.54 2.24
0.06 0.08 0.07
0.61 1.23 1.14
-0.077 -0.120 -0.145
0.74 1.80 2.43
0.06 0.09 0.10
0.14 0.73 1.11
0.0 -0.017 -0.023
0.79 1.48 1.99
0.0 -0.01 -0.02
0.61 0.98 1.18
0.0 -0.015 -0.013
0.60 1.72 1.59
0.0 -0.02 -0.02
0.41 1.02 0.94
0.0 -0.015 -0.026
0.75 1.21 1.81
0. 0 -0.01 -0.03
0.57 0.81 0.93
0.0 -0.003 -0.005
0.51 1.07 1.13
0.0 -0.01 -0.01
0.35 0.64 0.67
-0.028 -0.046 -0.049
1.27 1.58 2.00
-0.03 -0.03 -0.03
0.52 0.85 1.02
0.0 -0.018 -0.018
0.48 1.53 1.46
0.0 -0.02 -0.02
0.31 0.85 0.77
a
Single-point calculations on gas phase optimized geometry at the semiempirical level. The boundary elements (BE) grid density for the RF/CT calculations were set to 4 BE's/Å2 (see text). A second grid with a density of 20 BE's/Å2 was used to account for the curvature dependence of the surface charge density b The difference in the total semiempirical energy (eV) between the partially charged and neutral states of the solute in solution. c The ratio of the solvation free energy of the partially charged molecule to the experimental solvation free energy. d The amount of CT charge in eu. e The ratio between the solvation free energy of the neutral solute and the experimental solvation free energy. f Mulliken (Coulson) charges. g CM1 charges.21 h CM2 charges.22
CT charge density to be complementary to the pattern of the RF charge density. This is because the largest term in the energy functional defined in eq 9 is the one due to CT and RF charge interactions. Figure 1 shows the magnitudes of both RF and CT charge densities calculated for the hydroxyl anion with AM1/ CM1 at 4 BE/Å2 grid resolution. The hydroxyl ion transfers about -0.3 eu to the dielectric continuum and the CT charge is
distributed on BE's placed mainly on the surface of the oxygen atom. Figure 1 clearly indicates that the pattern of the CT charge densities complements the pattern of positive RF charge densities and that there is no negative charge transferred on BE's with negative RF charge density (mainly belonging to the hydrogen atom). This is because the interaction between the CT and RF charges would be repulsive.
Charge Flow between Ions and a Dielectric Continuum
J. Phys. Chem. B, Vol. 104, No. 9, 2000 2121
TABLE 3: Dependence of CT Amounta on the Density of the Boundary Element Grid grid densityb compound
1.2
2.0
4.0
0.6
8.0
10.0
water ammonia benzene phenol methanol acetic acid acetamide CH3COOCH3NH3+
0.10 0.13 0.26 0.42 0.09 0.14 0.18 -0.12 0.20
0.05 0.03 0.09 0.20 0.03 0.08 0.09 -0.14 0.17
0.04 0.01 0.08 0.17 0.02 0.08 0.05 -0.10 0.15
0.04 0.02 0.08 0.13 0.01 0.08 0.05 -0.09 0.15
0.03 0.02 0.06 0.11 0.01 0.08 0.04 -0.06 0.15
0.03 0.01 0.04 0.11 0.01 0.08 0.04 -0.05 0.15
a
In eu. b In boundary elements (BE) per Å2.
TABLE 4: Comparison of CPU Times for the Uniformly (UCT) and Variationally Spread-on-the-Grid (VCT) Charge Methodsa CPUb
CPUb
compound
UCT
VCT
compound
UCT
VCT
water CH3NH3+
200 504
230 2657
CH3COOC6H5O-
914 2366
4036 22425
a Single-point calculations at gas phase geometry. The reaction field calculations were performed at a grid density of 4 BE/Å2. b In seconds.
Figure 1. CT calculation (AM1/CM1) on the hydroxyl ion in a polar dielectric media ( ) 80): the magnitude of the reaction field (full line) and CT (dotted line) charges as a function of the position of the boundary element on the molecular surface. The boundary elements, which have CT charges, are placed on the surface of the oxygen atom.
CT between Ions and the Dielectric. In paper I we showed that the CT between ions and the dielectric continuum obtained with the UCT method ranges from -0.07 eu for C2H- (PM3) to -0.39 eu for NH2- (PM3/CM1/CM2). For cations CT ranges from 0.10 eu for (CH3)2SH+ (AM1/CM2) to 0.40 eu for C6H5NH3+ (PM3/CM2). For 18 ions out of 32 in our test set, the amount of CT charge changes little (maximum (0.04 eu) by using the VCT method. For seven ions the amount of CT charge is larger and for the final 11 smaller (average (0.10 eu) than the values obtained with the UCT method (Table 1). A significant drop in the amount of CT charge (AM1/CM1) was found, e.g., for CH3O- (-0.23 f -0.08 eu) and C6H5O(-0.17 f -0.08 eu). It is worth noting that the effect on aliphatic and aromatic ammonium ions is opposite. For example, the amount of CT for CH3NH3+ drops from 0.21 to 0.15 eu, while for C6H5NH3+ increases from 0.36 to 0.43 eu. A large increase in CT is observed for ions with two functional groups, for example, 0.13 eu for HC(OH)NH2+, 0.11 eu for imidazo-
lium, and 0.15 eu for pyridinium. The change in solvation free energy due to CT follows the same trend in both methods; i.e., the solvation free energy increases for anions and decreases for cations (Table 1). CT between Neutral Molecules and Dielectric. We did not obtain CT with the UCT method for neutral molecules (e.g., water), but we did determine it using the VCT method. For example, the amount of CT charge for water (AM1/CM1) is 0.04 eu at 4 BE/Å2 (Table 2) and 0.03 eu at 8 BE/Å2 grid resolution. Interestingly, a full quantum mechanical calculation on an equilibrated MD snapshot configuration of a box of water molecules shows that the charge on water molecules fluctuates between -0.04 and 0.06 eu, but the bulk of the water molecules have charges between (0.02 eu.24 The amount of CT for the majority of the neutral molecules in our set is small (less than 0.05 eu). Larger CT amounts were found (AM1/CM1 at 4 BE/Å2, Table 2) for benzene (0.09 eu), phenol (0.17, and 0.11 eu at 8 BE/Å2, Table 3), acetic acid (0.09 eu), and 4-methylpyridine (0.08 eu). Four out of 24 neutral molecules in our set show (small) negative CT (Table 2). The solvation free energy for all neutral molecules increases with the amount of CT (Table 2). Unlike ions, which are discharged, neutral molecules are charged by CT. Consequently, the reaction field energy (the largest component of solvation free energy) increases for neutral molecules for either direction of CT (i.e., from molecule to dielectric or vice versa). Table 4 shows a comparison between CPU times required by the UCT and VCT methods. For water the difference in CPU time is small but increases rapidly with the size of the molecule. Thus, for CH3NH3+ and CH3COO- the CPU time required by VCT method is 4 times larger than that required by the UCT method, while for C6H5O- it is 10 times larger. The reason for the steep increase of the CPU time for the VCT method comes from the way the CT grid is constructed. After the CT charge densities are calculated, the algorithm removes all boundary elements from the CT grid which have CT charge with the same sign as the RF charge. This step is iterated until no further BE's are removed. As the grid enlarges for larger molecules, it appears that constructing the CT grid becomes the calculation bottleneck. Future work will focus on increasing the efficiency of this step. Finally, we would like to comment on the earlier work of Chipman25,26 and Klamt and Jonas27 on charge penetration and dielectric volume polarization which may bear some resemblance with this work. Solute charge penetration into the dielectric continuum is observed in methods which mix the dielectric continuum method with a quantum mechanical Hamiltonian. The charge penetration arises due to the spread of solute wave function tail into the dielectric and the extent of charge penetration is closely related to the nature of the set of basis functions used in the ab initio calculation. The use of diffuse basis functions will, for example, considerably extend the penetration of solute charge into the dielectric. The collapse of the wave function into atomic charges eliminates the charge penetration effect, but Chipman25 argues that the volume polarization will be neglected. In Figure 3 of paper I,10 we showed that an energy decomposition of the total electrostatic energy of the solute in the continuum also displays a minimum at about the same amount of CT as that given by all-atom semiempirical calculations. Hence, CT can be modeled entirely within the framework of classical electrostatics. Our method was devised to simulate a quantum effect, i.e., the transfer of charge between two quantum systems and we believe that it has little connection with the charge renormalization issue discussed above. For example, the charge penetration into the
2122 J. Phys. Chem. B, Vol. 104, No. 9, 2000 dielectric still belongs to the solute while in our method the charge is transferred to the solvent by modifying the number of electrons in the solute frontier orbitals. Another aspect concerns the way in which the effect of charge penetration is treated. Chipman introduces an additional apparent surface charge to account for volume polarization of the dielectric, while in the case of our method the charge transferred to the dielectric is real and this is reflected in the new boundary conditions of the Poisson-Boltzmann equation. In other words, the CT in our method charges the dielectric while the penetration charge in the Chipman approach only additionally polarizes the dielectric. Conclusion We have presented in this paper an algorithm for distributing CT charge into a dielectric, which complements our effective Hamiltonian-Poisson-Boltzmann equation approach for studying CT between a quantum mechanical solute and a dielectric continuum.10 The present method eliminates the need to specify which atoms of the solute will “host” the CT grid, and the algorithm is especially useful for cases where the selection of “appropriate” atoms for constructing the CT grid is not straightforward. This is definitely the case for large molecules such as proteins or DNA where the polarization of the electronic charge density is significant and can render the biased scheme used in the UCT method ineffective. Comparison between the calculations performed on a set of 36 ions with the UCT and VCT methods shows that the amount of CT charge changes little for more than half of ions. For the rest of them the amount of CT either increases or decreases. In contrast to our former UCT calculations we found CT for neutral molecules with the VCT method. The amount of CT charge for neutral molecules is generally small with a few notable exceptions (phenol, acetic acid). Finally, we remark that the VCT method is much more expensive than the UCT method, as the size of the molecule increases, and its practical usage for large systems depends on optimization of the VCT algorithm. Work to find such solutions is currently underway in our group. Acknowledgment. We thank the DOE (DE-FGO296ER62270) for supporting this research. We also thank the
Gogonea and Merz Pittsburgh Supercomputer Center, the San Diego Supercomputer Center, the National Center for Supercomputer Applications, and the Cornell Theory Center for generous allocations of supercomputer time. References and Notes (1) Dixon, S. L.; Merz, K. M. Jr. J. Chem. Phys. 1996, 104, 66436649. (2) Dixon, S. L.; Merz, K. M. Jr. J. Chem. Phys. 1997, 107, 879-93. (3) Dixon, S. L.; Merz, K. M. Jr. The Divide and Conquer Methodology Applied to Semiempirical SCF MO Calculations. In Encyclopedia of Computational Chemistry Schleyer, P. v. R., Ed.; Wiley: Baffins Lane, Chichester, UK, 1998; Vol. 1, pp 762-776. (4) Millam, J. M.; Scuseria, G. E. J. Chem. Phys. 1997, 106, 55695577. (5) Yang, W.; Lee, T.-S. J. Chem. Phys. 1995, 103, 5674-5678. (6) Nadig, G.; van Zant, L. C.; Dixon, S. L.; Merz, K. M. Jr. J. Am. Chem. Soc. 1998, 120, 5593-5594. (7) van der Vaart, A.; Merz, K. M. Jr. Int. J. Quantum Chem., in press. (8) van der Vaart, A.; Merz, K. M. Jr. J. Phys. Chemistry A 1999, 103, 3321-3329. (9) van der Vaart, A.; Merz, K. M. Jr. J. Am. Chem. Soc., submitted for publication. (10) Gogonea, V.; Merz, K. M. Jr. J. Chem. Phys., in press. (11) Jackson, J. D. Classical Electrodynamics; John Wiley: New York, 1962; Vol. 1. (12) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027-2094. (13) Tannor, D. J.; Marten, B.; Murphy, R.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Ringnalda, M.; William, A. Goddard, I.; Honig, B. J. Am. Chem. Soc. 1994, 116, 11875-82. (14) Gogonea, V.; Merz, K. M. Jr. J. Phys. Chem. A 1999, 63, 517188. (15) Zauhar, R. J.; Morgan, R. S. J. Mol. Biol. 1985, 186, 115. (16) Zauhar, R. J.; Morgan, R. S. J. Comput. Chem. 1988, 9, 171-187. (17) Zauhar, R. J.; Varnek, A. J. Comput. Chem. 1996, 17, 864-877. (18) van der Vaart, A.; Merz, K. M. Jr. J. Am. Chem. Soc. 1999, 121, 9182-90. (19) Press, W. H.; Teokolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes. The Art of Scientific Computing, IInd ed.; Cambridge University Press: Cambridge, UK, 1992. (20) Rashin, A. A. J. Phys. Chem. 1990, 94, 1725-1733. (21) Storer, J. W.; Giesen, D. J.; Cramer, C. J.; Truhlar, D. G. J. Comput.Aided Mol. Des. 1995, 9, 87-110. (22) Li, J.; Zhu, T.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1998, 102, 1820-1831. (23) Connolly, M. L. J. Appl. Crystallogr. 1985, 18, 499. (24) Gogonea, V.; Merz, K. M. Jr., unpublished results. (25) Chipman, D. M. J. Chem. Phys. 1997, 106, 10194-206. (26) Zhan, C.-G.; Bentley, J.; Chipman, D. M. J. Chem. Phys. 1998, 108, 177-92. (27) Klamt, A.; Jonas, V. J. Chem. Phys. 1996, 105, 9972-81.
