Modeling Charge Flux by Interpolating Atomic Partial Charges

Article Cite This: J. Chem. Inf. Model. 2019, 59, 2837−2849

pubs.acs.org/jcim

Modeling Charge Flux by Interpolating Atomic Partial Charges Seung Soo Kim† and Young Min Rhee*,‡ †

Department of Chemistry, Pohang University of Science and Technology (POSTECH), Pohang 37673, Korea Department of Chemistry, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea

‡

Downloaded via NOTTINGHAM TRENT UNIV on August 13, 2019 at 20:40:58 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Although the charge flux effect or the geometric dependence of the atomic partial charges have been known for a long time, how it can be effectively handled is not yet established. Here, we present a charge interpolation scheme as a new general tool for representing the charge flux in an analytically well-defined manner. By applying it to the anionic GFP chromophore with the diabatically represented atomic charges, we show that the charge interpolation provides a substantial improvement on the accuracy of the geometrydependent changes in the molecular dipole moments in the gas phase. We also test the scheme toward describing the electrostatic term in the solvation energy in the aqueous environment and observe that it is also improved but that the extent of the improvement is somewhat limited. We show that the remaining errors can be largely corrected by introducing atomic polarizabilities. Overall, our results show that charge interpolation is an amenable approach for describing the charge flux effect and that its description in the condensed phase should be accompanied by proper treatments of polarization effects. studying dynamics of complex systems, such as a fluorescent protein19 and a light-harvesting complex.25 It was also extended to allow simulations of nonadiabatic processes26,27 by employing a diabatic representation of PES28 together with similarly defined atomic partial charges.29 Up to this point, however, the electrostatic interaction between the IM and the MM regions have been based on fixed atomic charges, which sometime induce substantial errors. Because the atomic charge of a molecule will change upon the variation in the surrounding electrostatic environment and upon the distortion of the molecular geometry, neglecting the charge response in the molecule will inevitably cause inaccuracies in the PES description.24,29 The first type of change, namely the polarization effect, has indeed been aimed as a primary task in futuristic force field developments.30−33 Various force fields concerning the polarization effect have been reported34−48 adopting the different polarization models: the fluctuating charge model,49,50 the Drude oscillator model,51,52 and the induced point dipole model.2,53,54 The nonadditive nature of the polarization effects sometimes makes it difficult to find the optimal parameters with these models. A previous analysis on the induced point dipole model actually pointed out that the accuracy of the parametrization was markedly lower than the maximally achievable level,55 and thus treating polarization effects in a generalized way still remains a challenge. Meanwhile, the geometric dependence of the electronic density in a given

1. INTRODUCTION Molecular dynamics (MD) is a universal simulation tool and helps gain detailed insights from atomistic motions of molecules in diverse systems. It indeed has been used to model a wide variety of phenomena such as protein folding, enzymatic catalysis, and even photobiological processes.1 Its accuracy critically depends on the quality of the mathematical expression of the involved potential energy surface (PES) containing all possible interactions in the system, namely the force field. Often, the terms in the force fields rely on simplified forms of interactions partly in order to reduce the associated computational costs. If necessary, sophisticated terms based on quantum chemical calculations can replace some parts of the force field. Of course, the direct quantum mechanics/ molecular mechanics (QM/MM) will be an exemplary approach.2,3 However, it is often hurdled by excessive computational costs. Accordingly, more efficient methods have been devised over the past few decades.4,5 One prominent example will be the empirical valence bond approach,4,6 together with its semiempirical7−12 or ab initio13−17 extensions. Our group has also developed an alternative approach named interpolation mechanics/molecular mechanics (IM/ MM)18,19 to circumvent the issue related to the choice of analytic functional forms. It employs local interpolations of quantum chemical data20−23 to describe the PES function of the important part of a system, in combination with the MM description of its environment. The advantage of IM/MM lies on its efficiency and its formal simplicity, as it avoids both the on-the-fly quantum calculations and the rather laborious PES fitting.24 IM/MM has indeed been successfully applied to © 2019 American Chemical Society

Received: April 11, 2019 Published: May 24, 2019 2837

DOI: 10.1021/acs.jcim.9b00307 J. Chem. Inf. Model. 2019, 59, 2837−2849

Article

Journal of Chemical Information and Modeling

2.1. IM/MM Approach. Suppose that one is interested in modeling a multistate PES with a number of adiabatic electronic states and that all the transition densities are spatially localized within the IM region. The IM/MM method basically employs the same partitioning scheme as the additive QM/MM approach.3 With this, the overall PES matrix Utot of the system can be written as

molecule accounts for another portion of electrostatic variations. This is sometimes referred to as the charge flux,56 and it has been shown to be important in modeling IR intensities,56−64 hydrogen bonding,64−66 and conformational energies.67−69 One study showed that the charge flux can contribute to 40% of the total electrostatic force in some cases.70 More recent studies also showed that the charge flux accounts for 10−15% of the geometric dependence of molecular dipole moments. 71,72 However, the general importance of the charge flux in condensed phase is still elusive because a universal model for treating its effect is still absent. Some force fields that can deal with this charge flux have been reported.45,66,70 As is usual with conventional force fields, they need to rely on parametric fitting, and how to effectively assign the charge flux parameters may become an issue. In principle, at least one parameter is needed for each internal coordinate and each atomic partial charge. Thus, as the number of parameters grows, redundancy and/or overfitting issues may arise. While such an issue can be circumvented by choosing a small number of atoms whose charge variations are important together with the coordinates associated with them, a recent analysis on the mode contribution to the charge flux in a set of organic molecules suggested that all types of internal coordinates should be treated within the same theoretical footage to produce a quantitatively reliable charge flux model.72 This aspect clearly calls for a more systematic and efficient way of handling charge flux effects. In fact, the IM/MM scheme can provide a natural way of treating charge fluxes. Namely, just like interpolating geometrydependent energies, one can also interpolate each atomic partial charge according to the simple Shepard’s method. This tactic of interpolating atomic partial charges can certainly be an attractive way of describing charge flux effects as it can avoid the arbitrariness of choosing some parametrized atoms or the hassles linked with the parametrization itself. Moreover, it should be applicable to a wide range of systems with minimal additional efforts as all required information naturally come from quantum chemical calculations that are needed for energy interpolation. Here, we construct this charge interpolation scheme and test the resulting charge flux model. As a test case, we will adopt the anionic green fluorescent protein (GFP) chromophore molecule in combination with its already constructed IM/MM potential data. How well the charge flux effect is reproduced will be checked with two limiting properties: (i) the geometry-dependent molecular dipole moments in vacuum and (ii) the electrostatic interaction energies in water. Indeed, we will find excellent reproductions of the geometry-dependent dipoles with the interpolation. In addition, we will observe some systematic errors in the interaction energies in water, which we will ascribe to the environmental polarization effect. This will be confirmed by considering atomic polarizabilities. We will conclude by considering future prospects for further advancing the IM/ MM technique by including polarization.

Utot = (UIM + Ucoul) + (Udisp + UMM + Uh)I

(1)

Here, UIM is the potential energy matrix of the IM region and Ucoul is the state-specific electrostatic interaction energy between the IM region and the MM region. In addition, Udisp is the Lennard-Jones dispersion interaction energy between the two regions. UMM is the potential energy of the MM region and Uh is the MM-style stabilizer function18 for the IM region, which is identically applied to all electronic states. Of course, I is the identity matrix with the dimensionality of the considered electronic states. Note that in eq 1, the matrix form is adopted to highlight the difference between the statedependent and the state-independent terms. Although the dispersion term can be expressed in a state specific manner, it usually does not vary much depending on electronic states and is included as a state-independent term in the equation. Because the interpolation in the above handles multiple electronic states, the explicit definitions of UIM + Ucoul can vary depending on how we define the state basis. To avoid singularities around surface crossing points, it is advantageous to employ (quasi-)diabatic representations of PES’s for this purpose. The diabaitc potential energy matrix D can be expressed as D = RTVR

(2)

where V is the adiabatic potential energy matrix and R is the adiabatic-to-diabatic transformation (ADT) matrix. There are many methods to practically obtain R.73−76 Regardless of the choice, however, its gradient should satisfy the following condition: ∂R = −Fk R ∂Xk

(3)

Here Xk is a Cartesian component of the molecular geometry vector X and Fk is its corresponding nonadiabatic coupling matrix. By definition, its elements can be written as (Fk )ab = ⟨Ψa|

∂ |Ψb⟩ ∂Xk

(4)

with the wave functions Ψa and Ψb of the electronic states a and b. Thus, with properly defined R, the diabatic potential energy matrix and its derivatives can be computed from {Fk}, V, and the corresponding derivatives. We should note that conventional quantum chemical tools provide information for the adiabatic states in the direct sense. Now, we define the diabatic potential energy matrix of the IM region DIM as a weighted average of the Taylor-expanded diabatic potential energy matrices from a set of data point geometries:

2. METHODS Before explaining the charge interpolation scheme, let us first briefly restate the diabatic version of IM/MM. In the following two parts, we will illustrate how the current IM/MM approach is operated to represent PES with multiple electronic states. More detailed explanations can be found in earlier reports.18,23,28,29

DIM (X) =

∑ wn(X)Dn(X) n

(5)

Here Dn(X) is the second-order Taylor-expansion of the diabatic potential energy centered at an interpolating data geometry n. The interpolation weight wn(X) is defined as 2838


Article

Journal of Chemical Information and Modeling wn(X ) =

[1/dn(X)]2p 2p

∑m [1/dm(X)]

UIM + Ucoul as a whole should match V′IM. Again, all forces and nonadiabatic couplings can be computed by differentiating eq 12. One should note that the fixed diabatic charge model already contains some charge flux effect. This is because the adiabatic atomic partial charge matrix pα at an arbitrary molecular α geometry X is converted from qref with the geometry dependent ADT matrix R:

(6)

where dn(X) is the Cartesian distance function from n. The definition of p is not so crucial as long as it is large enough and we chose p = 6. More explicit formulations of eqs 5 and 6 can be found elsewhere.23,28 Once all Dn(X) is calculated quantum chemically, DIM at an arbitrary molecular geometry can be back-transformed into VIM by matrix diagonalization: VIM(X) = R(X)DIM (X)R(X)T

pα(X) = R(X)q αref R(X)T

(7)

Thus, one can derive the charge flux from eqs 3 and 13 as

Note here that the ADT matrix R(X) is found in an on-the-fly manner. The dependence on X will be omitted hereafter unless doing so induces an ambiguity. The adiabatic force can be derived by differentiating this equation. Also, the nonadiabatic coupling on the IM region without the MM part can be obtained as (Fk )ab = −

∂pα = −Fk pα + pαFk ∂Xk

(8)

2.2. Fixed Diabatic Charge Model. In this part, we will describe how Ucoul is treated in relation to the diabatic interpolation scheme. Within a fixed diabatic charge model, the charge distribution of the IM region is represented by a set of atomic partial charge matrices {qαref} determined on one representative molecular geometry Xref. Because the partial change information is obtained also from quantum chemical calculations, which normally provide numbers based on adiabatic representations, similar to eq 2, the diabatic atomic partial charge matrix of the αth atom can be defined as q αref = R Tref pαref R ref

(9)

q αintp(X) =

Here, we explicitly used subscripts to emphasize that the matrices are obtained at Xref. The adiabatic atomic partial charge matrix pαref can be expressed as

(10) α

where N is the total number of electrons, V is the volume spanned by the αth atom, and {ri} represents the electronic coordinates. Certainly, there is no unique and exact way to define Vα, and somewhat different approaches of atomically assigning charges can be employed to obtain pαref. Once the set {pαref} becomes available, converting it to {qαref} is trivial. This set of diabatic charges are then used to calculate the electrostatic perturbation on DIM due to the external MM charges. The perturbed diabatic potential energy matrix D′IM can be written as D′IM = DIM +

α q αref ∑ ϕMM α

The introduction of the charge flux effect with eq 15 or 16 affects the diabatic potential in eq 11 and subsequently the associated force and the nonadiabatic coupling. Calculating the adiabatic force will involve the charge flux computation, which can be attained with a back-transformation from qαintp(X) as

(11)

with ϕαMM denoting the Coulombic potential from the MM charges on the IM atom α. This combined diabatic potential energy matrix should be diagonalized to give the adiabatic potential matrix. Namely, V′IM = R′D′IM R′T

(15)

where qnα is the diabatic atomic partial charge matrix at a data point geometry n. One can think of this equation as the zerothorder version of the local interpolation.20 In principle, qαn can be a Talyor-expanded diabatic atomic partial charge just like Dn(X) in eq 5. However, the gradient of the diabatic atomic partial charge matrix is expected to be small because the charge distributions of diabatic states change slowly, especially in comparison with the adiabatic one. In a later section, we will indeed check the error arising from the neglect of ∂qα/∂X|n by adding this gradient to eq 15 such that the expansion order for the local interpolation becomes one. In this case, the geometrydependent charge is expressed as ÄÅ ÉÑ α ÅÅ ÑÑ ÅÅ α ∂q Ñ α ·(X − X n)ÑÑÑ q intp(X) = ∑ wn(X)ÅÅq n + ÅÅ ÑÑ ∂X n ÅÇ ÑÖ (16) n

α

× Ψ*b (r1, r2, ···, rN ; X ref )

∑ wn(X)q nα n

∫V dr1∫ dr2dr3 ··· drNΨa(r1, r2, ···, rN ; X ref )

α (p ref )ab = N

(14)

with the gas-phase nonadiabatic coupling matrix Fk expressed in eq 8. 2.3. Interpolated Diabatic Charge Model. Let us now move on to the main topic of this work. In principle, the fixed diabatic charge model should give the most reliable statedependent charge distribution within the realm of fixed charge approximation. However, a previous study29 showed that the diabatic atomic partial charge itself has a non-negligible geometric dependence. Of course, this will cause an error in both the adiabatic atomic partial charges and the corresponding charge fluxes. One simple solution to this limitation is to include the geometry dependence, and interpolation will again be an attractive tactic toward this end. In the interpolated diabatic charge model, the diabatic atomic partial charge matrix qαintp is allowed to vary at different molecular geometries as

[R(∂DIM /∂Xk)RT]ab (VIM)aa − (VIM)bb

(13)

∂pα ∂ = [Rq αintp(X)RT] ∂Xk ∂Xk

ÅÄÅ ∂q α (X) ÑÉÑ ÅÅ intp ÑÑ T ÑÑR = −Fk p + p Fk + RÅÅÅÅ Ñ ÅÅ ∂Xk ÑÑÑ ÅÇ ÑÖ

(12)

It should be noted that the ADT matrix R′ here is different from the one in eq 7, as the diabatic potential includes the IMMM interaction. If eq 1 is represented in the adiabatic form,

α

2839

α

(17)


Article

Journal of Chemical Information and Modeling 2.4. Obtaining pαn. In the original proposition of the diabatic atomic partial charge formulation,29 due to the conceptual simplicity toward atomizing the one-particle density given in eq 10, the atomic charge from the natural bond orbital (NBO) analysis77 was adopted. However, the more widely adopted restrained electrostatic potential (RESP) fitting78 can still be applied if the electrostatic potential concept can be extended to multistate cases. There is one complication, however, in this case, as there is no electrostatic potential corresponding to the off-diagonal state with a ≠ b. For this, the “transition electrostatic potential” can be defined by adopting the transition density as if it is the regular oneparticle density. Then the same RESP procedures can be taken to obtain the transition atomic partial charges, on the condition that the sum of all transition atomic partial charges is zero. We note that a similar approach of obtaining transition atomic partial charges was previously reported,79 even though the purpose was somewhat different from ours. Indeed, it will be interesting to test which scheme of charge atomization works better for describing the charge flux effect. In order to actually verify this, we will adopt both RESP and NBO charges in the following. These tests will also illustrate the intrinsic differences in the two charge assigning schemes in terms of the qualities of the interpolated diabatic charges as well as the interpolated energies, which will be crucial items in performing dynamics simulations. As already mentioned, the order of Taylor expansion in the charge may or may not be an important parameter, and to test the extent of such importance, we considered the numerical gradients of the diabatic charges at all data point geometries. These gradients were transformed to fit with the diabatic formalism and added to the expansion according to eq 16. To distinguish the charge interpolation with the charge gradient terms explicitly, we will refer to the schemes as RESP-1 and NBO-1 models. Of course, the charge interpolation without the gradients will be designated as RESP-0 and NBO-0. When unambiguous, we will just refer to RESP-0/NBO-0 as interpolated RESP/NBO. 2.5. Computational Details. As a test case for our diabatic charge interpolation scheme, we will adopt the anionic form of the chromophore molecule in the green fluorescent protein (GFP). We chose this system as its charge indeed has a strong geometric dependence within the adiabatic picture.27 For this chromophore, the charge-localized diabatic representation can be generated by block diagonalization based on fragment-localized orbitals. 80 Following the earlier reports,26,29,80 we will label the three diabatic states as P, B, and I as they are dominated by electronic configurations with a negative charge on the phenoxy, bridge, and imidazolinone groups, respectively (Figure 1). All the interpolation data for the Talyor expansions of the diabatic potential energy matrices at 1500 geometries were taken from ref 26. These geometries are the data locations of the interpolation (n in eq 5), and for succinctness we will denote the set of these geometries as n1500. In generating the charges, one-particle density matrices at the SA(3)-CAS(4,3)SCF level were re-evaluated at these n1500 geometries together with the 6-31G(d) basis set. Each SA-CASSCF solution was checked carefully by confirming that the charge center of each fragment-localized orbital80 is located on the appropriate moiety. At each geometry, the RESP fitting and the NBO analysis were performed with the density, followed by the diabatization of the adiabatic atomic partial charges with the ADT matrix. Because we are handling a three-

Figure 1. (a) Chemical structure of the anionic GFP chromophore in water with the adopted atom and group designations. (b) Schematic descriptions of three diabatic states. Each chemical structure represents the dominant electronic configuration in the corresponding diabat.

dimensional electronic system with six independent matrix elements, in total, we will be using 9000 sets of diabatic atomic partial charges toward the charge interpolation. Next, to evaluate the errors associated with our charge interpolation scheme, 2000 test geometries were additionally taken from ref 29. These geometries were actually generated from 500 excited state IM/MM MD trajectories that contained solvent water as the MM region. In estimating the interpolation error with this set, because comparing atomic charges for all atoms will be cumbersome, we will employ the chromophore molecular dipoles. In comparing the dipoles, the SA-CASSCF level gas phase molecular dipole moments were calculated after omitting all MM water molecules. At 63 out of the 2000 test geometries, the SA-CASSCF solutions converged inappropriately, and instead of reattempting toward improved convergences, we simply eliminated them from the test set. Thus, 1937 geometries with proper SA-CASSCF solutions were left and adopted for the tests. For succinctness, this set will be referred to as T-1937. In addition to the dipole, for additional comparisons, the QM/MM interaction energies between the chromophore and the MM water were calculated at the SACASSCF level. In some cases (87 out of 1937 geometries), the inclusion of MM external charges caused breakdowns of SACASSCF solution arising from lowering the S2 state energy. In fact, this type of breakdown was already mentioned in an earlier report.80 For these cases, we estimated the QM/MM interaction energies using frozen inactive orbitals from the matching SA(2)-CAS(4,3)SCF solutions. All multireference quantum chemical calculations were performed with the MOLPRO package (version 2010.1),81 and atomic charge analyses were carried out with the help of the developers’ version of Q-Chem82 bridged with NBO 5.0.83 IM/MM calculations were performed through GROMACS 4.0.784 linked with the homemade libraries for executing diabatic potential energy and charge interpolations. 2840


Article


3. RESULTS AND DISCUSSION 3.1. RESP versus NBO. Quite naturally, by construction, the diabatic atomic charge depends on the molecular electronic structure. As already noted, in the original construction of the diabatic charge formalism,29 NBO concept was utilized. However, it is known that fitting against the electrostatic potential with RESP performs better than adopting natural population at least for describing a situation with a single electronic state and is thus more widely adopted. Therefore, we extend the RESP formalism here such that it can handle multiple diabatic electronic states as mentioned in the earlier section. Now we will demonstrate how the two schemes behave with the charge interpolation technique that we have developed. The diabatic RESP charges of the anionic GFP chromophore at the S0-optimized geometry, after grouping as in Figure 1, are listed in Table 1. The constituting atomic charges are

The intrinsic disparity in the NBO scheme is further noticed with the molecular dipole moment. Figure 2 compares the predicted diabatic state dipole moment and diabatic transition dipole moment against the reference values from the quantum chemical calculations. From the figure, one can clearly see that the RESP dipole moments are nearly flawless for all geometries in the n-1500 data set, while NBO dipoles are somewhat deviating. Indeed, the root-mean-squared errors in dipole moments and transition dipole moments are, respectively, 0.104 and 0.011 D with RESP, whereas the errors are 1.692 and 0.081 D with NBO. Interestingly, the NBO scheme displays a tendency of overestimating the diabatic state dipole moments (Figure 2a). This discrepancy appears to have arisen from incorrect bond polarization in NBO predominantly with the two C−C bonds near the bridge moiety (Figure S1). From the geometry-dependent partial charge information collected from the n-1500 set, we can also infer how RESP and NBO atomic charges respond to the geometry changes. Figure 3 shows the standard deviations of diabatic RESP and NBO charges obtained within the set. On most atoms, the RESP charges display larger standard deviations for all electronic states. Thus, we can infer that RESP responds better to the geometry change and this will be one of the reasons RESP behaves much better in reproducing molecular dipoles in Figure 2. In a later part, we will relate this difference in the ability of handling the charge changes at different molecular geometries to the fidelity of describing electrostatic interactions. 3.2. Charge Flux with Charge Interpolation. On the basis of the diabatic charge matrices obtained in the above, we built the interpolated diabatic charge model with the zeroth and the first-order Taylor expansions of the charges, as formulated in the previous section. In the case of the first-order expansion, calculating gradients of RESP charges was not pursued because of its high computational cost. Instead, we approximated the RESP charge gradient with the NBO charge gradient. We believe this is a reasonable approximation as the charge flux, or the geometric dependence of the charge change, from the two schemes are relatively well-correlated (Figure S2). For benchmarking, we will again use the molecular dipole moment, but we will momentarily switch to the dipoles in the adiabatic basis as the dipole itself changes more drastically with the geometry change for the adiabatic states.

Table 1. Group Diabatic Charges at the S0-Optimized Geometry of the GFP Chromophore method

groupa

P statea

B statea

I statea

RESP

GP GB GI GP GB GI

−0.814 −0.206 0.022 −0.997 0.149 −0.153

−0.042 −1.067 0.109 −0.446 −0.147 −0.406

−0.064 −0.135 −0.803 −0.186 0.158 −0.972

NBO

a

Group and state designations can be found in Figure 1.

listed in Tables S1 and S2. Not surprisingly, the group charge distributions are in line with the charge localization pattern shown in Figure 1. Strikingly, the RESP charges are quite different from the NBO charges. The difference is the largest on the bridge moiety especially with the diabatic B state. For instance, the RESP group charge on GB is −1.067e, while it is only −0.147e with NBO. The same analyses were performed for all n-1500 elements, and they showed that the NBO scheme tends to assign less population on the bridge moiety in comparison with the RESP scheme (Figure S1). Apparently, the RESP group charges are in better accord with the charge localization pattern expected from our diabatization scheme.80

Figure 2. Comparison of adopting RESP (red) and NBO (black) atomic partial charges in predicting (a) diabatic state dipole moments and (b) diabatic transition dipole moments, with the n-1500 set. μchg is the dipole from RESP/NBO parametrizations while μqm is the reference values from the QM one-particle densities. 2841


Article


Figure 3. Standard deviations of the diabatic atomic partial charges for (a) P state, (b) B state, and (c) I state. The standard deviations were calculated based on the n-1500 set. The state and atom designations are given in Figure 1.

The error in the S1 state dipole moment as a function of the molecular geometry is shown in Figure 4. The errors were projected on the two-dimensional conformational space in the vicinity of the S0-optimized geometry. The space was spanned by a torsional angle and its coupled bond length as depicted in the figure. These two degrees of freedom actually change to large extents upon photoexcitation,85 and the adiabatic charge distribution also changes quite drastically with their distortions.27 From Figure 4a,b, we can indeed observe that the adiabatic dipole moments predicted from the fixed charge models without any charge interpolation display non-negligible geometry-dependent errors, whether the fixed charges come from RESP or NBO. At large distortions from the S0-optimized geometry, the errors generally increase. The RMS errors from the two schemes are 0.852 and 3.012 D, respectively, similar to what was observed with Figure 2. The error from the RESP scheme is much smaller simply due to the improvement in the accuracy of the diabatic atomic partial charges. When we adopted the charge interpolation scheme with the zeroth-order expansions in charge, the geometric dependence of the dipole was better reproduced as is visibly clear in Figure 4c,d, with the RMS errors estimated as 0.263 and 2.424 D, respectively. What is more important than the numerical improvement is the fact that the deviation at large displacements from the equilibrium geometry is almost eliminated with the RESP interpolation such that the uncertainty in the electrostatic nature is almost homogeneous in space (Figure 4c). When the expansion order is raised to one by adopting atomic charge derivatives, the improvement is actually unnoticeable as shown in Figure 4e,f compared to Figure 4c,d. The RMS errors from the two schemes are 0.475 and 2.164 D, respectively. In fact, in some cases, the deviation at large displacements from the equilibrium geometry becomes larger as with the charge derivatives. This is likely because of the extrapolation effect. Namely, when the conformation reaches a point where no interpolation data point is close enough, derivative components may dominate the interpolated

terms. For interpolated energies, there are mechanisms that can avoid the associated artifacts such as using stabilizer,18 but devising a similar mechanism is not needed for the charge if we omit derivatives. Similar behaviors are observed for the S0 state dipole moment (Figure S3). To more carefully check the potential extrapolation issue, we further evaluated the dipole error at geometries collected from MD simulations, namely with the T-1937 set. With dynamic simulations, there will be geometries at fairly large distances from all the data point geometries. As already mentioned, T1937 geometries were sampled from the excited state IM/MM simulations. For each sampled conformation, the distance from the closest data point was calculated to correlate with the interpolation error in the state-averaged dipole. This is pictorially displayed in Figure 5 and one can clearly see that the error can become large at a large distance. In most cases, the errors with the T-1937 set is much reduced with the adoption of the RESP interpolation (Table S3). This is basically in accord with what we have already observed with Figure 4. While the improvement in the average error listed in Table S3 may not appear drastic, this is because the average error is dominated by the conformations at large distances. However, when we consider only conformations that are within 1 Å distance from any data point in Figure 5b,c, the improvement is quite noticeable. One additional interesting aspect that we can observe from Table S3 is the fact that the average error in transition dipole moments actually increases with the first-order expansion in charge. This happens due to the extrapolating effect for largely deviating conformations and indeed shows that adopting RESP-1 can potentially be harmful as expected from Figure 4e,f. Overall, we note that the diabatic charge still depends on the molecular geometry, and thus considering this charge flux effect is actually important at least in reproducing the molecular dipole moment and the related electrostatic interactions. Including the zeroth-order expansions with the diabatic charge model was enough for reconstructing the 2842


Article


Figure 4. Error in the S1 state adiabatic dipole moment estimated with (a) the fixed RESP model, (b) the fixed NBO model, (c) the RESP-0 model, (d) the NBO-0 model, (e) the RESP-1 model, and (f) the NBO-1 model. The error is projected on the rp−ϕp plane around the S0-optimized geometry. The S0-optimized geometry is marked with arrows, located at (rp, ϕp) = (1.4 Å, 0 deg).

of modes should be considered with equal footing in treating the charge flux. 3.3. Electrostatic Interaction Energy with the Environment. Although the molecular dipole can serve as a useful metric for characterizing the charge distribution, the electrostatic interaction energy is also important as it ultimately affects the dynamics in the condensed phase. Therefore, let us now consider how the charge flux affects this interaction energy. To probe this, we calculated the IM/MM electrostatic interaction energies by adopting the series of test geometries of

dipole moment. Including the charge gradient in the expansion for the interpolation actually should be used with care when we consider the potentially extrapolating situations encountered during simulations. Perhaps considering higher-order derivatives of the charges beyond the gradient level could be helpful in improving the accuracy, but the associated increase in the computational cost will not justify this additional work. From Figure 4, we can also observe that the charge flux depends both on the bond length and the dihedral angle. This aspect is in accord with a previous report,72 which suggested that all types 2843


Article


Figure 5. Correlations between the minimum Cartesian distance from the interpolating data points in the n-1500 set and the error in the stateaveraged dipole moment. Each point represents error with a test geometry belonging to the T-1937 set. The dipoles are calculated (a) without the charge interpolation, (b) with the RESP-0 interpolation, and (c) with the RESP-1 interpolation.

important in governing nonadiabatic molecular dynamics, the surface hopping behaviors described with the NBO charges earlier26,27 will likely have had similar reliabilities as what could have been attained with the RESP charges. Of course, this may not be a general event and applying RESP should likely be preferred even for nonadiabatic molecular simulations with multiple electronic states. There is one additional important aspect that we should note from Table 2. The improvement by adopting the charge interpolation is now only marginal and is not as visible as the one provided with Figures 4 and 5. One may argue that this is a consequence of the extrapolation effect and the limitation of the interpolation scheme itself. However, when we inspected the correlation between the error in the electrostatic interaction energy and the distance between the test geometry and its most closely neighboring data point, the error was rather insensitive to the distance (Figure 6a,b). Indeed, there were cases where the error in the interaction energy was quite large even when the molecular dipole from charge interpolation was quite accurate. This is indeed striking and maybe puzzling when we consider the success of the quality of the dipole moments with RESP shown in Figure 2. Thus, again, we can infer that the source of error resides in the external part of the chromophore molecule. In fact, unlike QM/MM, our interpolation scheme does not have an external polarization mechanism. Namely, the charge distribution in the

the chromophore (T-1937) solvated by water. These were compared against the reference values obtained directly from QM/MM calculations. The test geometries contained MM water molecules from the beginning without any additional conformational adjustments for the test. For fair comparison between the IM/MM and the QM/MM electrostatics, the periodic boundary condition and the cutoff in long-range electrostatics were removed in the IM/MM calculations to match with the conditions in the QM/MM calculations. The RMS errors of the electrostatic interaction energies with various IM/MM schemes are presented in Table 2. Interestingly, IM/MM schemes with RESP charges show similar or even slightly lower reliabilities as the ones with NBO charges. This is in contrary to what is expected from the results in Figures 2 and 4, suggesting that some cancellation of errors occurs with NBO. Also, the gap energies are reproduced slightly better with NBO. Because the gap energy is primarily Table 2. RMS Errors in the IM/MM Electrostatic Interaction Energies in kJ/mol for the T-1937 Set method

S0

S1

S2

S0-S1 gap

S0-S2 gap

fixed RESP RESP-0 fixed NBO NBO-0

51.6 50.8 54.4 52.1

54.0 52.2 51.6 51.9

53.2 54.0 50.0 49.1

15.1 12.9 12.2 11.3

17.1 15.3 12.7 12.3 2844


Article


Figure 6. Correlations between the minimum Cartesian distance from interpolation data set and the error in the state-averaged IM/MM electrostatic interaction energy with (a) the RESP-0 charge interpolation and (b) the NBO-0 charge interpolation, based on the T-1937 set. The same correlations after the polarization correction with (c) the RESP-0 charge interpolation and (d) the NBO-0 charge interpolation.

interpolated chromophore does not change with the field from the aqueous environment, and we can imagine this lack of physicality can contribute importantly to the errors. To confirm this conjecture, we calculated the isotropic atomic polarizabilities at the S0-optimized geometry following a simplified approach based on the scheme in ref 86. Then, the energies from the interpolated charge model were corrected by including the polarization energies from the water environment. To reduce the complexity in modeling, we restricted ourselves to the additive induced point dipole model that neglects all interactions between the induced dipoles. The details of how we treated the polarization effect are described in the Supporting Information. The change in the electrostatic interaction energies with this added polarization is schematically provided in Figure 6c,d. In the case of the RESP-based IM/MM, certainly, much of the systematic error can be fixed with this rather simple correction. This proves that the polarization of the IM part by the MM environment indeed contributes to the errors in RESP-based IM/MM observed in Table 2. In contrast, the same correction applied to the NBObased IM/MM scheme did not result in any improvement (Figure 6d). This deterioration suggests that the apparent accuracy of the NBO-based IM/MM shown in Table 2 actually came from the cancellation of errors between the intrinsic dipole deviation in NBO and the missing contribution of polarization. Overall, we anticipate that including external

polarization effect should be the key task for the future development of the IM/MM approach. Basically, considering the charge flux effect should be accompanied by properly describing molecular polarizations when handling molecular complexes or solvation.

4. CONCLUDING REMARKS In this work, we developed the diabatic charge interpolation model to probe the importance of charge flux effect in developing the IM/MM potential. By applying it to the anionic GFP chromophore, we checked the accuracies of the gas phase dipole moments and of the aqueous phase electrostatic interaction energies. In the case of the dipole moments, the diabaic charge interpolation has a substantial effect on the accuracy, demonstrating its utility toward describing the charge flux effect. The electrostatic interactions were also improved, but the extent of this improvement was only marginal. In other words, even with the much improved gas phase molecular dipoles, there were still significant errors in the solvation energies. We showed that the external polarization indeed was the main source of the remaining errors. Thus, the future development of the IM/MM scheme will have to be on implementing the external polarization capability to the scheme. Our approach of interpolating the atomic partial charges toward describing the charge flux effect has the advantage of 2845


Article

Journal of Chemical Information and Modeling being simple toward modeling the charge flux effect itself. We do not have to assume any functional forms, and the algorithm can be seamlessly implemented to the existing interpolation code set. While we tested it with a multistate electronic system within the diabatic realm of representation, there is no reason the model will not be applicable to simpler systems for which only single electronic states are considered. Namely, interpolating adiabatic charges should be equally feasible. Indeed, even with the GFP chromophore that we have tested, the ground state charge flux is relatively well represented with the interpolation of adiabatic atomic partial charges as shown in Figure S4. (We note that because the charge distribution in the GFP chromophore strongly depends on the molecular geometry within the adiabatic representation, RESP-1 performed much better than RESP-0 in this case.) Thus, we anticipate that our method can become a generally useful tool for carefully treating electrostatic interactions of flexible systems for both diabatic and adiabatic representations. The current interpolating scheme has required 1500 reference data points to properly model potential energy and the atomic partial charges of the GFP chromophore. This may seem somewhat undesirable considering the relatively small size of the chromophore. However, we emphasize that the amount of calculations will still be much less than what will be required for direct dynamics simulations at the QM level when lengthy simulations are required especially for handling condensed phase situations. Moreover, the wall clock time for the QM calculations can be effectively reduced by performing the QM calculations in a parallel manner. In our experience, this aspect indeed reduces the total computational time for QM calculations at least by an order of magnitude. Of course, in comparison to IM/MM without charge interpolation, the added computational cost results from treating the IM region. For the IM region with N atoms modeled by M interpolation data points, the processor FLOP count for the PES interpolation with its gradient computation can be shown to be ∼4M(3N)2 for each potential energy element. Actually, the most time-consuming part comes from treating Hessians from all reference data points for both energy and force evaluations. In the case of the charge interpolation with RESP0 or NBO-0, most of the additional time is spent in calculating {∂qintpα/∂Xk} and scales as ∼2M(3N)N in terms of the FLOP counts. Thus, the charge interpolation itself is not that burdensome and the FLOP count remains as O(MN2). To further check this, we performed 20000 steps MD simulation with the GFP chromophore in 2000 explicit water molecules and confirmed that the computational overhead by additionally performing charge interpolation was ∼10%. This was ∼5% of all IM/MM time, and as the size of the MM region grows the fractional cost will of course diminish. In the end, the cost for IM/MM is negligible compared to what is needed for constructing the interpolation data set. Thus, we conclude that the current interpolation scheme may again provide an efficient route for performing simulations at the QM/MM-like level at the pure MM-like speed, to be in line with a previous report.26 Including polarization will further complete the scheme, and we hope to report on this in due course.

■

■

Procedures for treating the polarization effect; RESP and NBO diabatic charges at the S0 geometry (Tables S1 and S2); RMS errors in diabatic dipole moments and transition dipole moments for the T-1937 set (Table S3); the diabatic atomic partial charges averaged over the n-1500 set (Figure S1); correlations between RESP and NBO diabatic charge displacement in the n-1500 set (Figure S2); errors in the S0 state adiabatic dipole moment estimated with various charge treating schemes (Figure S3); errors in the S0 state dipole moment when the charge interpolation is performed with adiabatic atomic partial charges without any diabatization on quantum chemical data (Figure S4) (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Seung Soo Kim: 0000-0002-1179-9765 Young Min Rhee: 0000-0002-2392-3962 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was financially supported by the Midcareer Researcher Program (Grant 2017R1A2B3004946) and by the Creative Materials Discovery Program (Grant 2018M3D1A1058813) through National Research Foundation, funded by Ministry of Science and ICT of Korea.

■

REFERENCES

(1) Karplus, M.; McCammon, J. A. Molecular Dynamics Simulations of Biomolecules. Nat. Struct. Biol. 2002, 9, 646−652. (2) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (3) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (4) Kamerlin, S. C. L.; Warshel, A. The Empirical Valence Bond Model: Theory and Applications. WIREs Comput. Mol. Sci. 2011, 1, 30−45. (5) Shurki, A.; Sharir-Ivry, A. Valence Bond-Based Hybrid Quantum Mechanics Molecular Mechanics Approaches and Proper Inclusion of the Effect of the Surroundings. Isr. J. Chem. 2014, 54, 1189−1204. (6) Åqvist, J.; Warshel, A. Simulation of Enzyme Reactions Using Valence Bond Force Fields and Other Hybrid Quantum/Classical Approaches. Chem. Rev. 1993, 93, 2523−2544. (7) Kim, Y.; Corchado, J. C.; Villà, J.; Xing, J.; Truhlar, D. G. Multiconfiguration Molecular Mechanics Algorithm for Potential Energy Surfaces of Chemical Reactions. J. Chem. Phys. 2000, 112, 2718−2735. (8) Higashi, M.; Truhlar, D. G. Electrostatically Embedded Multiconfiguration Molecular Mechanics Based on the Combined Density Functional and Molecular Mechanical Method. J. Chem. Theory Comput. 2008, 4, 790−803. (9) Higashi, M.; Truhlar, D. G. Combined Electrostatically Embedded Multiconfiguration Molecular Mechanics and Molecular Mechanical Method: Application to Molecular Dynamics Simulation of a Chemical Reaction in Aqueous Solution with Hybrid Density Functional Theory. J. Chem. Theory Comput. 2008, 4, 1032−1039. (10) Vuilleumier, R.; Borgis, D. Molecular Dynamics of an Excess Proton in Water Using a Non-Additive Valence Bond Force Field. J. Mol. Struct. 1997, 436−437, 555−565.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.9b00307. 2846


Article

Journal of Chemical Information and Modeling (11) Schmitt, U. W.; Voth, G. A. Multistate Empirical Valence Bond Model for Proton Transport in Water. J. Phys. Chem. B 1998, 102, 5547−5551. (12) Chakraborty, A.; Zhao, Y.; Lin, H.; Truhlar, D. G. Combined Valence Bond-Molecular Mechanics Potential-Energy Surface and Direct Dynamics Study of Rate Constants and Kinetic Isotope Effects for the H+C2H6 Reaction. J. Chem. Phys. 2006, 124, 044315. (13) Shurki, A.; Crown, H. A. Hybrid Ab Initio VB/MM Method - A Valence Bond Ride through Classical Landscapes. J. Phys. Chem. B 2005, 109, 23638−23644. (14) Sharir-Ivry, A.; Crown, H. A.; Wu, W.; Shurki, A. Density Embedded VB/MM: A Hybrid Ab Initio VB/MM with Electrostatic Embedding. J. Phys. Chem. A 2008, 112, 2489−2496. (15) Mo, Y.; Gao, J. Ab Initio QM/MM Simulations with a Molecular Orbital-Valence Bond (MOVB) Method: Application to an SN2 Reaction in Water. J. Comput. Chem. 2000, 21, 1458−1469. (16) Song, L.; Mo, Y.; Gao, J. An Effective Hamiltonian Molecular Orbital-Valence Bond (MOVB) Approach for Chemical Reactions as Applied to the Nucleophilic Substitution Reaction of Hydrosulfide Ion and Chloromethane. J. Chem. Theory Comput. 2009, 5, 174−185. (17) Ying, F.; Chang, X.; Su, P.; Wu, W. VBEFP: A Valence Bond Approach That Incorporates Effective Fragment Potential Method. J. Phys. Chem. A 2012, 116, 1846−1853. (18) Park, J. W.; Kim, H. W.; Song, C.-I.; Rhee, Y. M. Condensed Phase Molecular Dynamics Using Interpolated Potential Energy Surfaces with Application to the Resolvation Process of Coumarin 153. J. Chem. Phys. 2011, 135, 014107. (19) Park, J. W.; Rhee, Y. M. Interpolated Mechanics−Molecular Mechanics Study of Internal Rotation Dynamics of the Chromophore Unit in Blue Fluorescent Protein and Its Variants. J. Phys. Chem. B 2012, 116, 11137−11147. (20) Ischtwan, J.; Collins, M. A. Molecular Potential Energy Surfaces by Interpolation. J. Chem. Phys. 1994, 100, 8080−8088. (21) Thompson, K. C.; Jordan, M. J. T.; Collins, M. A. Polyatomic Molecular Potential Energy Surfaces by Interpolation in Local Internal Coordinates. J. Chem. Phys. 1998, 108, 8302−8316. (22) Bettens, R. P. A.; Collins, M. A. Learning to Interpolate Molecular Potential Energy Surfaces with Confidence: A Bayesian Approach. J. Chem. Phys. 1999, 111, 816−826. (23) Rhee, Y. M. Construction of an Accurate Potential Energy Surface by Interpolation with Cartesian Weighting Coordinates. J. Chem. Phys. 2000, 113, 6021−6024. (24) Rhee, Y. M.; Park, J. W. Interpolation for Molecular Dynamics Simulations: From Ions in Gas Phase to Proteins in Solution. Int. J. Quantum Chem. 2016, 116, 573−577. (25) Kim, C. W.; Choi, B.; Rhee, Y. M. Excited State Energy Fluctuations in the Fenna−Matthews−Olson Complex from Molecular Dynamics Simulations with Interpolated Chromophore Potentials. Phys. Chem. Chem. Phys. 2018, 20, 3310−3319. (26) Park, J. W.; Rhee, Y. M. Towards the Realization of Ab Initio Dynamics at the Speed of Molecular Mechanics: Simulations with Interpolated Diabatic Hamiltonian. ChemPhysChem 2014, 15, 3183− 3193. (27) Park, J. W.; Rhee, Y. M. Electric Field Keeps Chromophore Planar and Produces High Yield Fluorescence in Green Fluorescent Protein. J. Am. Chem. Soc. 2016, 138, 13619−13629. (28) Park, J. W.; Rhee, Y. M. Constructing Polyatomic Potential Energy Surfaces by Interpolating Diabatic Hamiltonian Matrices with Demonstration on Green Fluorescent Protein Chromophore. J. Chem. Phys. 2014, 140, 164112. (29) Park, J. W.; Rhee, Y. M. Diabatic Population Matrix Formalism for Performing Molecular Mechanics Style Simulations with Multiple Electronic States. J. Chem. Theory Comput. 2014, 10, 5238−5253. (30) Ponder, J. W.; Case, D. A. Force Fields for Protein Simulations. Adv. Protein Chem. 2003, 66, 27−85. (31) Cieplak, P.; Dupradeau, F.-Y.; Duan, Y.; Wang, J. Polarization Effects in Molecular Mechanical Force Fields. J. Phys.: Condens. Matter 2009, 21, 333102.

(32) Lopes, P. E. M.; Roux, B.; MacKerell, A. D. Molecular Modeling and Dynamics Studies with Explicit Inclusion of Electronic Polarizability: Theory and Applications. Theor. Chem. Acc. 2009, 124, 11−28. (33) Demerdash, O.; Yap, E.-H.; Head-Gordon, T. Advanced Potential Energy Surfaces for Condensed Phase Simulation. Annu. Rev. Phys. Chem. 2014, 65, 149−174. (34) Banks, J. L.; Kaminski, G. A.; Zhou, R.; Mainz, D. T.; Berne, B. J.; Friesner, R. A. Parametrizing a Polarizable Force Field from Ab Initio Data. I. The Fluctuating Point Charge Model. J. Chem. Phys. 1999, 110, 741−754. (35) Stern, H. A.; Kaminski, G. A.; Banks, J. L.; Zhou, R.; Berne, B. J.; Friesner, R. A. Fluctuating Charge, Polarizable Dipole, and Combined Models: Parameterization from Ab Initio Quantum Chemistry. J. Phys. Chem. B 1999, 103, 4730−4737. (36) Kaminski, G. A.; Stern, H. A.; Berne, B. J.; Friesner, R. A.; Cao, Y. X.; Murphy, R. B.; Zhou, R.; Halgren, T. A. Development of a Polarizable Force Field for Proteins via Ab Initio Quantum Chemistry: First Generation Model and Gas Phase Tests. J. Comput. Chem. 2002, 23, 1515−1531. (37) Patel, S.; Brooks III, C. L. CHARMM Fluctuating Charge Force Field for Proteins: I Parameterization and Application to Bulk Organic Liquid Simulations. J. Comput. Chem. 2004, 25, 1−16. (38) Patel, S.; Mackerell, A. D.; Brooks III, C. L. CHARMM Fluctuating Charge Force Field for Proteins: II Protein/Solvent Properties from Molecular Dynamics Simulations Using a Nonadditive Electrostatic Model. J. Comput. Chem. 2004, 25, 1504−1514. (39) Lamoureux, G.; MacKerell, A. D.; Roux, B. A Simple Polarizable Model of Water Based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185−5197. (40) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D. Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430−5449. (41) Yu, H.; Hansson, T.; van Gunsteren, W. F. Development of a Simple, Self-Consistent Polarizable Model for Liquid Water. J. Chem. Phys. 2003, 118, 221−234. (42) Geerke, D. P.; van Gunsteren, W. F. On the Calculation of Atomic Forces in Classical Simulation Using the Charge-on-Spring Method to Explicitly Treat Electronic Polarization. J. Chem. Theory Comput. 2007, 3, 2128−2137. (43) Cieplak, P.; Caldwell, J.; Kollman, P. Molecular Mechanical Models for Organic and Biological Systems Going Beyond the Atom Centered Two Body Additive Approximation: Aqueous Solution Free Energies of Methanol and N-Methyl Acetamide, Nucleic Acid Base, and Amide Hydrogen Bonding and Chloroform/Water Partition Coefficients of the Nucleic Acid Bases. J. Comput. Chem. 2001, 22, 1048−1057. (44) Wang, Z.-X.; Zhang, W.; Wu, C.; Lei, H.; Cieplak, P.; Duan, Y. Strike a Balance: Optimization of Backbone Torsion Parameters of AMBER Polarizable Force Field for Simulations of Proteins and Peptides. J. Comput. Chem. 2006, 27, 781−790. (45) Palmo, K.; Mannfors, B.; Mirkin, N. G.; Krimm, S. Potential Energy Functions: From Consistent Force Fields to Spectroscopically Determined Polarizable Force Fields. Biopolymers 2003, 68, 383−394. (46) Ren, P.; Ponder, J. W. Consistent Treatment of Inter- and Intramolecular Polarization in Molecular Mechanics Calculations. J. Comput. Chem. 2002, 23, 1497−1506. (47) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549−2564. (48) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046−4063. (49) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358− 3363. 2847


Article

Journal of Chemical Information and Modeling (50) Rick, S. W.; Stuart, S. J.; Berne, B. J. Dynamical Fluctuating Charge Force Fields: Application to Liquid Water. J. Chem. Phys. 1994, 101, 6141−6156. (51) Jacucci, G.; McDonald, I. R.; Singer, K. Introduction of the Shell Model of Ionic Polarizability into Molecular Dynamics Calculations. Phys. Lett. A 1974, 50, 141−143. (52) Straatsma, T. P.; McCammon, J. A. Molecular Dynamics Simulations with Interaction Potentials Including Polarization Development of a Noniterative Method and Application to Water. Mol. Simul. 1990, 5, 181−192. (53) Applequist, J.; Carl, J. R.; Fung, K.-K. Atom Dipole Interaction Model for Molecular Polarizability. Application to Polyatomic Molecules and Determination of Atom Polarizabilities. J. Am. Chem. Soc. 1972, 94, 2952−2960. (54) Vesely, F. J. N-Particle Dynamics of Polarizable StockmayerType Molecules. J. Comput. Phys. 1977, 24, 361−371. (55) Li, A.; Voronin, A.; Fenley, A. T.; Gilson, M. K. Evaluation of Representations and Response Models for Polarizable Force Fields. J. Phys. Chem. B 2016, 120, 8668−8684. (56) Decius, J. C. An Effective Atomic Charge Model for Infrared Intensities. J. Mol. Spectrosc. 1975, 57, 348−362. (57) van Straten, A. J.; Smit, W. M. A. Bond Charge Parameters from Integrated Infrared Intensities. J. Mol. Spectrosc. 1976, 62, 297− 312. (58) King, W. T.; Mast, G. B. Infrared Intensities, Polar Tensors, and Atomic Population Densities in Molecules. J. Phys. Chem. 1976, 80, 2521−2525. (59) Gussoni, M.; Ramos, M. N.; Castiglioni, C.; Zerbi, G. Ab Initio Counterpart of Infrared Atomic Charges. Chem. Phys. Lett. 1987, 142, 515−518. (60) Gussoni, M.; Castiglioni, C.; Ramos, M. N.; Rui, M.; Zerbi, G. Infrared Intensities: From Intensity Parameters to an Overall Understanding of the Spectrum. J. Mol. Struct. 1990, 224, 445−470. (61) Haiduke, R. L. A.; Bruns, R. E. An Atomic Charge-Charge FluxDipole Flux Atom-in-Molecule Decomposition for Molecular DipoleMoment Derivatives and Infrared Fundamental Intensities. J. Phys. Chem. A 2005, 109, 2680−2688. (62) Silva, A. F.; Richter, W. E.; Meneses, H. G. C.; Faria, S. H. D. M.; Bruns, R. E. How Accessible Is Atomic Charge Information from Infrared Intensities? A QTAIM/CCFDF Interpretation. J. Phys. Chem. A 2012, 116, 8238−8249. (63) Torii, H. Mechanism of the Secondary Structure Dependence of the Infrared Intensity of the Amide II Mode of Peptide Chains. J. Phys. Chem. Lett. 2012, 3, 112−116. (64) Baranović, G. Intramolecular Resonance-Assisted Hydrogen Bonds: A Theoretical Description by Means of Atomic Charges and Charge Fluxes. Spectrochim. Acta, Part A 2014, 117, 465−472. (65) Dinur, U.; Hagler, A. T. The Role of Nonbond and Charge Flux in Hydrogen Bond Interactions. The Effect on Structural Changes and Spectral Shifts in Water Dimer. J. Chem. Phys. 1992, 97, 9161−9172. (66) Mannfors, B.; Palmo, K.; Krimm, S. Spectroscopically Determined Force Field for Water Dimer: Physically Enhanced Treatment of Hydrogen Bonding in Molecular Mechanics Energy Functions. J. Phys. Chem. A 2008, 112, 12667−12678. (67) Mannfors, B. E.; Mirkin, N. G.; Palmo, K.; Krimm, S. Analysis of the Pyramidalization of the Peptide Group Nitrogen: Implications for Molecular Mechanics Energy Functions. J. Phys. Chem. A 2003, 107, 1825−1832. (68) Palmo, K.; Mannfors, B.; Mirkin, N. G.; Krimm, S. Inclusion of Charge and Polarizability Fluxes Provides Needed Physical Accuracy in Molecular Mechanics Force Fields. Chem. Phys. Lett. 2006, 429, 628−632. (69) Rasmussen, T. D.; Ren, P.; Ponder, J. W.; Jensen, F. Force Field Modeling of Conformational Energies: Importance of Multipole Moments and Intramolecular Polarization. Int. J. Quantum Chem. 2007, 107, 1390−1395. (70) Dinur, U.; Hagler, A. T. Geometry-Dependent Atomic Charges: Methodology and Application to Alkanes, Aldehydes, Ketones, and Amides. J. Comput. Chem. 1995, 16, 154−170.

(71) Sedghamiz, E.; Nagy, B.; Jensen, F. Probing the Importance of Charge Flux in Force Field Modeling. J. Chem. Theory Comput. 2017, 13, 3715−3721. (72) Sedghamiz, E.; Ghalami, F. Evaluating the Effects of Geometry and Charge Flux in Force Field Modeling. J. Phys. Chem. A 2018, 122, 4647−4653. (73) Pacher, T.; Cederbaum, L. S.; Köppel, H. Approximately Diabatic States from Block Diagonalization of the Electronic Hamiltonian. J. Chem. Phys. 1988, 89, 7367−7381. (74) Subotnik, J. E.; Yeganeh, S.; Cave, R. J.; Ratner, M. A. Constructing Diabatic States from Adiabatic States: Extending Generalized Mulliken−Hush to Multiple Charge Centers with Boys Localization. J. Chem. Phys. 2008, 129, 244101. (75) Cave, R. J.; Stanton, J. F. A Simple Quasi-Diabatization Scheme Suitable for Spectroscopic Problems Based on One-Electron Properties of Interacting States. J. Chem. Phys. 2016, 144, 054110. (76) Hoyer, C. E.; Parker, K.; Gagliardi, L.; Truhlar, D. G. The DQ and DQΦ Electronic Structure Diabatization Methods: Validation for General Applications. J. Chem. Phys. 2016, 144, 194101. (77) Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural Population Analysis. J. Chem. Phys. 1985, 83, 735−746. (78) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269−10280. (79) Madjet, M. E.; Abdurahman, A.; Renger, T. Intermolecular Coulomb Couplings from Ab Initio Electrostatic Potentials: Application to Optical Transitions of Strongly Coupled Pigments in Photosynthetic Antennae and Reaction Centers. J. Phys. Chem. B 2006, 110, 17268−17281. (80) Olsen, S.; McKenzie, R. H. A Diabatic Three-State Representation of Photoisomerization in the Green Fluorescent Protein Chromophore. J. Chem. Phys. 2009, 130, 184302. (81) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Györffy, W.; Kats, D.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bennie, S. J.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Lee, S. J. R.; Liu, Y.; Lloyd, A. W.; Ma, Q.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Miller, III, T. F.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Welborn, M. MOLPRO, version 2010.1, a Package of Ab Initio Programs; Cardiff University: Cardiff, UK, 2010. (82) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; Kuś, T.; Landau, A.; Liu, J.; Proynov, E. I.; Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J.; Woodcock, H. L.; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chang, C.-M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M.; DiStasio, R. A.; Do, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; FustiMolnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W. D.; Harbach, P. H. P.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T.-C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A. D.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S.-P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Neuscamman, E.; Oana, C. M.; Olivares-Amaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Russ, N. J.; Sharada, S. M.; Sharma, S.; Small, D. W.; Sodt, A.; Stein, T.; Stück, D.; Su, Y.-C.; Thom, A. J. W.; Tsuchimochi, T.; Vanovschi, V.; Vogt, L.; Vydrov, O.; Wang, T.; Watson, M. A.; Wenzel, J.; White, A.; 2848


Article

Journal of Chemical Information and Modeling Williams, C. F.; Yang, J.; Yeganeh, S.; Yost, S. R.; You, Z.-Q.; Zhang, I. Y.; Zhang, X.; Zhao, Y.; Brooks, B. R.; Chan, G. K. L.; Chipman, D. M.; Cramer, C. J.; Goddard, W. A.; Gordon, M. S.; Hehre, W. J.; Klamt, A.; Schaefer, H. F.; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xu, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J.-D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C.-P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; van Voorhis, T.; Herbert, J. M.; Krylov, A. I.; Gill, P. M. W.; HeadGordon, M. Advances in Molecular Quantum Chemistry Contained in the Q-Chem 4 Program Package. Mol. Phys. 2015, 113, 184−215. (83) Glendening, E. D.; Badenhoop, J. K.; Reed, A. E.; Carpenter, J. E.; Bohmann, J. A.; Morales, C. M.; Weinhold, F. NBO 5.0; Theoretical Chemistry Institute, University of Wisconsin: Madison, WI, 2001. (84) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (85) Zhao, L.; Zhou, P.-W.; Li, B.; Gao, A.-H.; Han, K.-L. NonAdiabatic Dynamics of Isolated Green Fluorescent Protein Chromophore Anion. J. Chem. Phys. 2014, 141, 235101. (86) Kim, H. W.; Rhee, Y. M. Molecule-Specific Determination of Atomic Polarizabilities with the Polarizable Atomic Multipole Model. J. Comput. Chem. 2012, 33, 1662−1672.

2849


Modeling Charge Flux by Interpolating Atomic Partial Charges

Recommend Documents