Article pubs.acs.org/JCTC
Probing the Importance of Charge Flux in Force Field Modeling Elaheh Sedghamiz,†,‡ Balazs Nagy,† and Frank Jensen*,† †
Department of Chemistry, Aarhus University, Langelandsgade 140, DK-8000 Aarhus, Denmark Department of Chemistry, University of Isfahan, Isfahan 81746-73441, Iran
‡
S Supporting Information *
ABSTRACT: We analyze the conformational dependence of atomic charges and molecular dipole moments for a selection of ∼900 conformations of peptide models of the 20 neutral amino acids. Based on a set of reference density functional theory calculations, we partition the changes into effects due to changes in bond distances, bond angles, and torsional angles and into geometry and charge flux contributions. This allows an assessment of the limitations of fixed charge force fields and indications for how to design improved force fields. The torsional degrees of freedom are the main contribution to conformational changes of atomic charges and molecular dipole moments, but indirect effects due to change in bond distances and angles account for ∼25% of the variation. Charge flux effects dominate for changes in bond distances and are also the main component of the variation in bond angles, while they are ∼25% compared to the geometry variations for torsional degrees of freedom. The geometry and charge flux contributions to some extent produce compensating effects.
1. INTRODUCTION Modeling the time-dependent behavior of biomolecular systems is an important element in many areas of research.1 The fundamental time step for atomistic models is ∼1 fs, and simulations covering micro- or milliseconds therefore necessitate a computational efficient energy function. Essentially all current production type molecular dynamics simulations rely on parametrized energy functions, usually denoted force fields, such as CHARMM,2 AMBER,3 and GROMOS.4 These force fields were proposed in the mid-1980s and have since undergone a series of reparameterizations. The bonded energy terms are taken as harmonic expressions or low-order Fourier series, while the nonbonded energy terms are partitioned into a nonpolar van der Waals term, usually taken as a Lennard-Jones potential, and a polar electrostatic term, described as fractional atomic charges interacting by a Coulomb expression.5 The continuing increase in computational resources allows longer simulations that reduce sampling errors, and the accuracy of the force field may in some cases become the limiting factor. Modern force fields, of which AMOEBA6 is one of the more widely parametrized and tested, employ terms beyond harmonic for the bonded terms and include up to quadrupole moments on all atoms for the electrostatic term, as well as including atomic dipole polarizability. The AMBER and CHARMM force fields have also been extended to include polarization, by means of atomic dipole polarizability7 and the Drude particle model,8 respectively. The molecular electrostatic potential (ESP) quantifies the nonpolarizable electrostatic interaction between molecules and can be calculated quite accurately by relatively low-level electronic structure methods, such as density functional theory (DFT) with a medium sized basis set. A popular approach is to use such electronic structure data as the reference for assigning © XXXX American Chemical Society
atomic charges, dipoles, and quadrupole moments within a force field environment.9 The atomic parameters can be obtained either by partitioning the ESP into atomic contributions, for example by the distributed multipole analysis (DMA)10 or the LoProp11 methods, or by fitting the parameters to the reference ESP sampled at a suitable number of points in space.12 The DMA and LoProp methods typically require up to atomic quadrupole moments in order to achieve a satisfactory reproduction of the reference ESP. The fitting approach determines the atomic parameters by minimizing the error in reproducing the reference data and will therefore necessarily provide a better representation of the reference data than the DMA or LoProp approaches for the same number of parameters. We have shown, however, that the parameter space for the fitting procedure rapidly becomes redundant, and roughly half of the parameters in an atomic charge, dipole, and quadrupole model can be removed without affecting the total fitting error.13 There are furthermore substantial degrees of freedom in choosing a set of nonredundant atomic parameters, i.e. a specific combination of charge, dipole, and quadrupole moments.14 The parameter redundancy raises the issue of how to select a consistent set of parameters. It is clear that different choices of atomic charge, dipole, and quadrupole parameters, which reproduce the ESP with essentially the same accuracy, will lead to different atomic forces when interacting with another molecule and thus potentially lead to different dynamics. Dinur has suggested that a consistent description of both the interaction energy and forces requires that a set of atomic multipole parameters to a given order should be constrained to Received: March 20, 2017 Published: June 21, 2017 A
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
20 natural amino acids capped with acetyl and N-methyl groups. The analyses are done using a density functional theory method with an augmented polarized double-ζ type basis set, but since the ESP is relatively insensitive to the level of theory, the general conclusions are unlikely to depend on this particular choice.
reproduce the molecular multipole moment to one order higher.15 A model with only atomic charges should thus reproduce the molecular dipole moment; a model with atomic charges and dipole moments should reproduce also the molecular quadrupole moment, etc. This implies that the sum of all atomic charges is equal to the total molecular charge, while the sum of all atomic multipole moments of a given order is zero. We note that in the quantum chemical topology approach,16 the sum of atomic multipole moments up to a given order only reproduces the molecular multipole moment to the same order, i.e. atomic charges reproduce the molecular charge, but atomic dipoles are required for also reproducing the molecular dipole moment. In some cases this leads to a large negative correlation between atomic parameters at different multipole orders,17 which is likely to be reduced by the present ESP fitting approach. Production type force fields describe the electrostatic energy by a Coulomb interaction of fixed fractional atomic charges, but it is well-known that the atomic charges depend on the molecular conformation,18 and this can be considered as an intramolecular polarization effect. From an order-expansion it would appear natural to describe the geometry dependence of the molecular dipole moment by allowing the atomic charges to depend on the (relative) atomic positions, but an efficient parametrization of this geometry dependence is an open question. The fluctuating charge method19 models the geometry dependence in terms of atomic electronegativity and hardness parameters and a Coulomb type interaction between all atom pairs. While this may be suitable for describing the intermolecular polarization, and intramolecular polarization between moieties that are far apart in terms of bonding, this may not be optimum for describing changes due to differences in bond distances or angles. The present work will attempt to quantify the relative importance of pure geometry and charge flux effects for different degrees of freedom in internal coordinates. The employed model is at the lowest expansion order in the Dinur classification,15a where atomic sites are only assigned a charge parameter, and these are constrained to reproduce the total molecular dipole moment. This charge-only model has inherent limitations, such as being unable to describe out-of-plane polarization in planar molecules. The next level would include also a set of atomic dipole and dipole flux parameters constrained to reproduce the molecular quadrupole moment. Our previous experience13,14 suggests that this may run into the problem of parameter redundancy, i.e. the ESP does not have sufficient information for a statistically valid fitting of both atomic charges and dipoles. We will in the present work assess limitations of fixed charge force fields and therefore work within a charge-only model. Molecular conformations can be characterized by a set of torsional angles, but it should be recognized that each conformation has its own set of optimum bond distances and angles. Within an atomic charge model constrained to reproduce the geometry dependent molecular dipole moment, we can ask to what extent the conformational dependent atomic charges is due to a direct effect, i.e. changes in the torsional angles, or due to an indirect effect, i.e. changes in bond distances and angles between each conformation. We will in the present work quantify the conformational dependence of atomic charges and partition the changes into contributions from differences in bond distances, bond angles, and torsional angles for a database consisting of ∼900 conformations of the
2. COMPUTATIONAL DETAILS Peptide geometries have been taken from ref 20 and consist of the 20 natural amino acids terminated by ACE and NME groups to create peptide models. All calculations have been done at the wB97XD/aug-pcseg-1 level21 using the Gaussian-09 program package.22 Atomic charges have been determined by the HLY method12c and have been constrained to reproduce the molecular dipole moment. The HLY method12c derives atomic charges by fitting to the molecular ESP and differs only from a number of similar methods12a,b by the selection of fitting points, where HLY employs a rather dense grid that ensures rotational invariance. The increase in energy by replacing optimized bond lengths by conformational averaged values is in all cases small (0.6 kJ/mol on average, standard deviation = 0.5 kJ/mol, maximum = 3.4 kJ/mol) which reflects mean absolute bond length changes of 0.0024 Å (standard deviation = 0.0027 Å, maximum = 0.038 Å). The increase in energy by replacing optimized bond angles by conformational averaged values in some cases leads to increases in energy by ∼100 kJ/mol due to steric clashes introduced by the geometry modifications. Conformations which lead to energy increases by more than 26 kJ/mol (0.01 hartree) were removed from the analysis, and average bond lengths and angles were recalculated from the reduced set of 896 conformations. The energy increase by employing average bond angles in this reduced set of conformations is 10.4 kJ/mol on average with standard deviation = 4.9 kJ/mol and maximum = 32.8 kJ/mol. These values reflect mean absolute bond angle changes of 1.2° (standard deviation = 1.1°, maximum = 11.7°). 3. THEORY We will in the following be working within an atomic chargeonly model, where the molecular dipole moment μk for a given conformation k can be calculated from the atomic charges qik and positions Rik as shown in eq 1. Natom
μk =
∑ qik R ik(rk, θk , ωk ) i=1
(1)
The atomic positions Rik depend on internal coordinates, where rk denotes (all) bond distances, θk denotes (all) bond angles, and ωk denotes (all) torsional angles. The atomic charges are determined by fitting to the DFT calculated ESP subject to the constraints that the total sum is zero and the molecular dipole moment is reproduced. The dipole moment modeled by fixed charge force fields has conformational independent atomic charges, which we will take as the conformational averaged values for each peptide model, and this provides conformational dependent dipole moments given by eq 2. μkq =
Natom
∑ qi R ik(rk, θk , ωk ) i=1
B
(2) DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation Nconf
1 qi = Nconf
∑ qik k=1
The difference between the dipole moments in eqs 1 and 2 can with eqs 7 and 11 be expanded as shown in eq 13.
(3)
qik = qi +
Δqikrθω
Natom
i=1
μkrθ =
i=1
(4)
Natom i=1 Natom
=
(8)
qikrθ = qirθ + Δqikω
(9)
The Δμ term represents a change in dipole moment arising from fixed atomic charges due to a change in bond distances (geometry change), Δμqr describes a change in atomic charge (charge flux) at a constant geometry, while Δμgqr describes a residual second order coupling between charges and bond distances. The first term on the right-hand side in eq 13 can be further expanded with the use of eqs 8 and 12. Natom
∑
=
− qi )
Natom
=
∑
[(qikrθ + Δqikθ − qi )(R ikrθ( r ̅ , θ ̅ , ωk ) + ΔR ikθ ( r ̅ , θk , ωk ))]
i=1 Natom
=
∑ i=1 Natom
=
∑
⎡(q rθ − q )R rθ( r , θ ̅ , ω ) + (q rθ − q )ΔR θ ( r , θ , ω ) +⎤ ik ̅ k ik ̅ k k i ik i ⎢ ik ⎥ ⎢ θ rθ ⎥ θ θ ⎣Δqik R ik ( r ̅ , θ ̅ , ωk ) + Δqik ΔR ik( r ̅ , θk , ωk ) ⎦ [(qikrθ − qi )R ikrθ( r ̅ , θ ̅ , ωk )] + Δμkgθ + Δμkqθ + Δμkgqθ
i=1
(14)
The Δμ , Δμ , and Δμ terms are geometry and charge flux terms associated with differences in bond angles and a residual second order coupling. The first term on the right-hand side in eq 14 can be further expanded with the use of eq 9. gθ
qθ
gqθ
Natom
∑ [(qikrθ − qi )R ikrθ( r ̅ , θ ̅ , ωk )] i=1
(10)
Natom
Since charge differences can be either positive or negative, we quantify the relative importance of each of the terms in eq 10 as an average of absolute values. The atomic charges determine the molecular electric moments, of which the dipole moment usually is the leading term for uncharged molecules. The conformational dependence of the dipole moment in eq 1 has contributions from changes in both the atomic charges and in the atomic coordinates, and these can analogously to the atomic charges be decomposed into contributions from bond distances, bond angles, and torsional angles. We define the geometry differences due to changes in bond lengths and angles (only) analogous to eqs 7 and 8.
=
∑ [(qirθ
+ Δqikω − qi)R ikrθ( r ̅ , θ ̅ , ωk )]
i=1 Natom
=
∑ [(qirθ
− qi )R ikrθ( r ̅ , θ ̅ , ωk ) + Δqikω R ikrθ( r ̅ , θ ̅ , ωk )]
i=1
= Δμkqqω + Δμkqω
(15)
The Δμqω term represents a charge flux due to changes in torsional angles (only), while Δμqqω can be considered as a residual term arising from the chosen decomposition involving conformational independent charge differences defined in eq 10. With the above-defined terms, the difference in dipole moments in eqs 1 and 2 can be written as shown in eq 16.
R ik(rk , θk , ωk ) = R ikr( r ̅ , θk , ωk ) + ΔR ikr(rk , θk , ωk )
μk − μkq = Δμkgr + Δμkqr + Δμkgqr + Δμkgθ + Δμkqθ
(11)
R ikr( r ̅ , θk , ωk ) = R ikrθ( r ̅ , θ ̅ , ωk ) + ΔR ikθ ( r ̅ , θk , ωk )
[(qikr − qi )R ikr( r ̅ , θk , ωk )]
i=1
Δqikrθω = Δqikr + Δqikθ + Δqikω + Δqiqq (qirθ
[(qikr − qi )R ikr( r ̅ , θk , ωk )] + Δμkgr + Δμkqr + Δμkgqr
(13)
From eqs 7−9 it follows that the conformational dependence of the atomic charges in eq 4 can be decomposed into terms arising from changes in bond distance, bond angles, torsional angles, and a residual conformational independent term.
Δqiqq
∑
⎡(q r − q )R ikr( r , θk , ωk ) + (q r − q )ΔR ikr(rk , θk , ωk ) +⎤ ̅ i ik i ⎢ ik ⎥ ⎢ Δq r R r ( r , θ , ω ) + Δq r ΔR r (r , θ , ω ) ⎥ ik k k k ⎣ ik ik ̅ k k ⎦ ik
gr
We can express the change in atomic charges due to changes in the bond distances and bond angles for a given conformation as shown in eqs 7 and 8, while eq 9 defines the change due to torsional angles (only), where the conformational averaged values are defined analogous to eq 3.
qikr = qikrθ + Δqikθ
[(qikr + Δqikr − qi )(R ikr( r ̅ , θk , ωk ) + ΔR ikr(rk , θk , ωk ))]
i=1
(6)
(7)
∑
=
(5)
qik = qikr + Δqikr
[(qik − qi )R ik(rk , θk , ωk )]
i=1
Natom
∑ qikrθ R ikrθ( r ̅ , θ ̅ , ωk )
∑
=
Natom
∑ qikr R ikr( r ̅ , θk , ωk )
∑ i=1
The goal of the present work is to analyze the conformational dependence of the charges in eq 4 and the difference between the dipole moments in eqs 1 and 2, i.e. the coupling between the leading electrostatic term and geometry. For this purpose, we will consider atomic charges and dipole moments calculated at geometries where all bond distances are constrained to their conformational averaged values, and where both bond distances and angles are constrained to their conformational averaged values. μkr =
Natom
μk − μkq =
The conformational dependence of the atomic charges as a correction relative to the average values is defined by eq 4, and this is the effect that is neglected in fixed charge force fields.
+ Δμkgqθ + Δμkqω + Δμkqqω
(12) C
(16)
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation Eq 16 allows a decomposition of the difference between a reference dipole moment and that modeled by conformational independent atomic charges into contributions from changes in bond length, bond angle, and torsional angle, and these can further be partitioned into geometry and charge flux terms. The dipole moment is a vector quantity and changes can be quantified in terms of a change in direction and a change in magnitude. The directional change depends on the alignment of the two sets of atomic coordinates, and it is always possible to find an alignment where the angle between the two dipole moments is identical zero. With a standard mass-weighted rootmean-squared alignment, the directional changes for the bond distance and angle dependence in eqs 13 and 14, and the torsional charge flux term in eq 15, are small with typical values of a few degrees, but since this choice of alignment is by no means unique, we will in the following focus on the changes in the magnitude of the dipole moment. The decomposition in eq 16 is exact, and the difference vector on the left-hand side can be written as a sum of eight vector terms on the right-hand side. When quantifying these by their magnitude, the equal sign no longer holds, i.e. the magnitude of the difference vector on the left-hand side is not the sum of the magnitudes of the eight terms on the right-hand side. The average of the magnitudes of the terms on the right-hand side should thus only be taken as an indication of the relative importance. The bond distance and angle second order terms Δμgqr and Δμgqθ are always significantly smaller (few percent) than the first order terms (Δμgr, Δμqr, Δμgθ and Δμqθ) and will not be discussed. Eq 16 allows an assessment of the μ - μq difference for a given conformation, but it should be recognized that the major part of the conformational variation of the dipole moment is correctly reproduced by the μq fixed charge model. In order to provide a measure of this pure geometrical effect, we define a scalar measure of the geometry torsional dependence of the dipole moment as shown in eq 17, where the last term is the conformational average length of the dipole moments using conformational average atomic charges. Δμkgω =
Natom
∑ i=1
qirθ R ikrθ( r ̅ , θ ̅ , ωk ) −
Natom
∑
Table 1. Atomic Charges for the Non-Hydrogen Atoms in a Representative Energy Conformation for the Alanine Peptide and Charge Differences Defined in Eq 10a atom
q
Δqrθω
Δqr
Δqθ
Δqω
C1 C2 O3 N4 C5 C6 C7 O8 N9 C10 av (|q|) %
−0.6393 0.8501 −0.5874 −0.6929 0.1856 −0.4758 0.5948 −0.5580 −0.3383 −0.3772 0.5299
−0.0383 0.0445 0.0106 0.0194 −0.1808 −0.0140 0.1514 −0.0258 0.0383 −0.1011 0.0624 100
−0.0015 0.0037 0.0003 −0.0060 −0.0015 0.0009 0.0068 0.0017 −0.0086 0.0016 0.0033 5
0.0138 0.0044 −0.0008 −0.0245 0.0405 0.0021 −0.0140 0.0051 −0.0045 −0.0040 0.0113 18
−0.0411 0.0393 0.0104 0.0411 −0.2136 −0.0191 0.1606 −0.0292 0.0431 −0.1034 0.0700 112
a Atoms are labeled starting from the ACE terminus. Superscripts r, θ, and ω indicate changes due to bond lengths, bond angles, and torsional angles, respectively.
Table 2. Atom and Conformational Averaged Values for the Absolute Charge Differences in Eq 10 for the 20 Amino Acid Peptide Modelsa
qirθ R ikrθ(r, θ , ωk )
i=1
(17)
This term will be denoted a torsional geometry term, in analogy with the bond distance and angle terms.
4. RESULTS AND DISCUSSION Table 1 illustrates the decomposition of atomic charges according to eq 10 for the non-hydrogen atoms in a representative conformation of the alanine peptide model. The conformational dependence of the charges is on average ∼10% of the absolute value but varies between 2% and 25% for the different atoms. The largest variation is due to the torsional degrees of freedom, but that effect alone in most cases overestimates the magnitude, with contributions from bond angle differences, and to a smaller extent bond distances, providing compensating effects. Table 2 shows values averaged over conformations and atoms for all of the 20 peptide models of the uncharged amino acids. All peptide models behave very similarly, and we will focus on the average over all 896 conformations. The torsional variation (Δqω) is the largest by 103% (as measured relative to the Δqrθω value), followed by the angle variation (Δqθ) with 20%, while the distance dependence (Δqr) of 4% is slightly
peptide
Nconf
Δqrθω
Δqr
Δqθ
Δqω
Δqqq
Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val av %
22 44 34 63 44 48 54 22 40 58 62 74 87 37 12 53 25 44 39 34 896
0.0494 0.0563 0.0524 0.0527 0.0696 0.0520 0.0489 0.0626 0.0488 0.0611 0.0514 0.0592 0.0597 0.0459 0.0498 0.0541 0.0478 0.0440 0.0483 0.0553 0.0535 100
0.0017 0.0017 0.0037 0.0026 0.0019 0.0025 0.0020 0.0020 0.0019 0.0015 0.0015 0.0017 0.0017 0.0016 0.0018 0.0024 0.0023 0.0015 0.0017 0.0016 0.0020 4
0.0086 0.0114 0.0108 0.0118 0.0126 0.0123 0.0101 0.0117 0.0088 0.0107 0.0109 0.0128 0.0119 0.0089 0.0109 0.0122 0.0129 0.0084 0.0094 0.0112 0.0109 20
0.0517 0.0579 0.0537 0.0540 0.0713 0.0538 0.0501 0.0631 0.0498 0.0634 0.0545 0.0589 0.0618 0.0476 0.0546 0.0557 0.0494 0.0455 0.0492 0.0566 0.0551 103
0.0037 0.0049 0.0052 0.0037 0.0055 0.0063 0.0024 0.0085 0.0024 0.0051 0.0047 0.0076 0.0032 0.0029 0.0065 0.0047 0.0059 0.0038 0.0047 0.0042 0.0048 9
a Superscripts r, θ, and ω indicate changes due to bond lengths, bond angles, and torsional angles, respectively.
smaller than the conformational independent term accounting for 9%. The charge decomposition thus suggests that the indirect effects of bond distance and angle changes with conformation are ∼25% of the total effect, and that these to some extent are in the opposite direction as the pure torsional effects, as illustrated in Table 1. Table 3 shows the Cartesian dipole components for the decomposition according to eq 16 for the same representative conformation for the alanine peptide model as in Table 1. The charge fluctuation component due to torsional angles (Δμqω) is the largest and in magnitude slightly larger than the total difference. The bond angle contributions are the second most D
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Table 3. Changes in the Cartesian Molecular Dipole Moment (Debye) for a Representative Conformation of the Alanine Peptidea component
μ
μ - μq
Δμgr
Δμqr
Δμgθ
Δμqθ
Δμqω
Δμqqω
X Y Z length %
−1.496 3.329 3.911 5.349
0.509 −0.335 −0.424 0.743 100
0.018 −0.006 −0.015 0.024 3
−0.019 −0.031 −0.032 0.049 7
0.011 −0.093 −0.104 0.140 19
−0.027 0.114 0.098 0.153 21
0.587 −0.353 −0.382 0.784 106
−0.066 0.033 0.014 0.075 10
a Superscripts g and q indicate geometry and charge flux terms, respectively, for changes in bond lengths (r), bond angles (θ), and torsional angles (ω), as defined in eqs 13−16.
Table 4. Changes in the Magnitude of Molecular Dipole Moments (Debye) Averaged over All Conformations of the 20 Amino Acid Peptide Modelsa peptide
Nconf
μ
μ - μq
Δμgω
Δμgr
Δμqr
Δμgθ
Δμqθ
Δμqω
Δμqqω
Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val av %
22 44 34 63 44 48 54 22 40 58 62 74 87 37 12 53 25 44 39 34 896
4.219 5.000 4.012 4.411 4.317 4.707 4.104 4.616 5.159 3.816 3.602 3.342 4.306 3.938 4.131 4.491 4.647 4.527 4.421 4.164 4.296
0.357 0.298 0.436 0.370 0.381 0.410 0.444 0.316 0.378 0.461 0.479 0.335 0.420 0.445 0.297 0.348 0.488 0.365 0.584 0.513 0.406 100
2.003 1.570 1.576 1.940 1.957 1.915 1.874 2.177 1.389 1.985 1.773 1.724 1.743 1.796 2.027 2.053 2.238 1.811 1.839 2.145 1.877
0.004 0.002 0.005 0.003 0.003 0.003 0.002 0.003 0.002 0.003 0.002 0.003 0.003 0.002 0.004 0.004 0.002 0.002 0.003 0.002 0.003 1
0.034 0.038 0.044 0.057 0.034 0.045 0.038 0.036 0.026 0.029 0.031 0.029 0.031 0.026 0.041 0.039 0.043 0.028 0.029 0.031 0.035 9
0.030 0.031 0.030 0.039 0.031 0.039 0.028 0.032 0.018 0.040 0.032 0.029 0.029 0.030 0.024 0.032 0.024 0.028 0.038 0.030 0.031 8
0.068 0.076 0.067 0.104 0.077 0.065 0.057 0.085 0.062 0.055 0.071 0.050 0.055 0.076 0.048 0.081 0.055 0.063 0.058 0.068 0.067 17
0.361 0.306 0.427 0.388 0.429 0.425 0.434 0.360 0.392 0.485 0.505 0.352 0.422 0.399 0.294 0.358 0.492 0.391 0.588 0.549 0.418 103
0.030 0.025 0.058 0.024 0.054 0.050 0.017 0.048 0.024 0.057 0.042 0.052 0.022 0.014 0.035 0.060 0.046 0.025 0.054 0.042 0.039 10
a Superscripts g and q indicate geometry and charge flux terms, respectively, for changes in bond lengths (r), bond angles (θ), and torsional angles (ω), as defined in eqs 13−17.
The decomposition into direct (torsional) and indirect (bond length and angle) effects from eq 16 is consistent with the results from atomic charge alone (Table 1), with bond angle effects (Δμgθ + Δμqθ) being ∼25% and bond length effects (Δμgr + Δμqr) being ∼10% of the total effect (μ - μq), while the torsional charge flux term (Δμqω) of almost equal magnitude as the total effect, indicating substantial compensations due to sign differences. Considering the two bond distance terms, Δμgr and Δμqr, Table 4 shows that the charge flux effect is an order of magnitude larger than the geometry effect. The charge flux effect (Δμqθ) also dominates over the geometry effect (Δμgθ) for the bond angle terms but only by a factor of 2. Using the definition of the torsional geometry effect in eq 17, the Δμgω effect is roughly a factor of 5 more important that the torsional charge flux term (Δμqω). This is the expected trend, with charge flux being more important for geometry changes directly affecting the molecular bonding structure, and is in agreement with previous findings.15b
important, while effects due to bond distances are the smallest. As for the atomic charge variations in Table 1, there are compensating effects from each of the different contributions. Table 4 shows the results for the 20 peptide models averaged over conformations, where μ is the reference dipole moment and μ - μq is the difference between the dipole moments according to eqs 1 and 2, i.e. the part that is neglected by fixed charge force fields. The Δμ components are defined in eqs 13−15 and represent a decomposition of the μ - μq term. The 20 peptide models behave very similarly, and only the averages over the 896 conformations will be discussed. The μ - μq difference is on average ∼10% of the reference dipole moment, with most values in the 2−25% range. We note that μq tend to be larger in magnitude than μ, but there are many exceptions where μq is smaller than μ (see the Supporting Information). The average magnitude of the dipole moment is 3.3−5.2 D for all the peptide models, but the total variation is between 0.5 and 10 D. There is for all the systems a significant conformational dependence; the dipole moment for the 22 conformations of the alanine peptide model, for example, varies between 1.2 and 7.4 D, and this is reflected in the Δμgω value of 2.0 D.
5. CONCLUSION We analyze the limitations of fixed charge force fields by decomposing the difference between a quantum mechanical E
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation calculated dipole moment and that calculated by a fixed charge model into geometry and charge flux terms and further decompose these into contributions from changes in bond distances, angles, and torsional angles for ∼900 conformations of amino acid peptide models. The largest component of the variation of the dipole moment with conformation is correctly reproduced by a pure geometry effect involving fixed partial atomic charges, with the missing component accounting for ∼10% of the total effect. The missing effect is in magnitude very similar to the charge flux component due to torsional angles only, with indirect effects due to changes in bond distances and angles being ∼10% and ∼25%, respectively. The summation to more than 100% shows that the individual components to some extent provide compensating effects, but the amount of compensation and which terms that partly cancel vary between conformations and systems. This implies that attempts of designing force fields to account for the geometry dependence of the molecular dipole moment (and implicitly also higher order electric moments) by including charge flux effects must be able to model the dependence on not only torsional degrees of freedom but also indirect effects due to changes in bond angles and, to a lesser extent, also changes in bond distances.
■
ring-conformational equilibria in hexopyranose-based carbohydrates chains. J. Comput. Chem. 2016, 37 (3), 354−365. (5) Jensen, F. Introduction to Computational Chemistry, 3rd ed.; Wiley: Chichester, UK, 2017. (6) 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 (9), 4046−4063. (7) 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 (10), 1048−1057. (8) 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 (12), 5430−5449. (9) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97 (40), 10269−10280. (10) Stone, A. J. Distributed multipole analysis: Stability for large basis sets. J. Chem. Theory Comput. 2005, 1 (6), 1128−1132. (11) Gagliardi, L.; Lindh, R.; Karlstrom, G. Local properties of quantum chemical systems: The LoProp approach. J. Chem. Phys. 2004, 121 (10), 4494−4500. (12) (a) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11 (4), 431−439. (b) Breneman, C. M.; Wiberg, K. B. Determining AtomCentered Monopoles from Molecular Electrostatic Potentials: The Need for High Sampling Density in Formamide ConformationalAnalysis. J. Comput. Chem. 1990, 11 (3), 361−373. (c) Hu, H.; Lu, Z.; Yang, W. Fitting Molecular Electrostatic Potentials from Quantum Mechanical Calculations. J. Chem. Theory Comput. 2007, 3 (3), 1004− 1013. (13) Jakobsen, S.; Jensen, F. Systematic Improvement of PotentialDerived Atomic Multipoles and Redundancy of the Electrostatic Parameter Space. J. Chem. Theory Comput. 2014, 10 (12), 5493−5504. (14) Jakobsen, S.; Jensen, F. Searching the Force Field Electrostatic Multipole Parameter Space. J. Chem. Theory Comput. 2016, 12 (4), 1824−1832 2999. (15) (a) Dinur, U. Force Related Atomic Multipoles in Planar Molecules: Derivation of Atomic Quadrupole and Octupole Moments. J. Comput. Chem. 1991, 12 (1), 91−105. (b) Dinur, U.; Hagler, A. T. Geometry-Dependent Atomic Charges: Methodology and Application to Alkanes, Aldehydes, Ketones, and Amides. J. Comput. Chem. 1995, 16 (2), 154−170. (16) (a) Bader, R. F. W. A Quantum-Theory of Molecular-Structure and Its Applications. Chem. Rev. 1991, 91 (5), 893−928. (b) Popelier, P. L. A. Atoms in Molecules: An Introduction; Pearson Education: Harlow, Great Britain, 1999. (17) 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 (31), 8238−8249. (18) (a) Koch, U.; Popelier, P. L. A. Characterization of C-H-O Hydrogen Bonds on the Basis of the Charge Density. J. Phys. Chem. 1995, 99 (24), 9747−9754. (b) Fletcher, T. L.; Popelier, P. L. A. Toward amino acid typing for proteins in FFLUX. J. Comput. Chem. 2017, 38 (6), 336−345. (19) (a) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular-Dynamics Simulations. J. Phys. Chem. 1991, 95 (8), 3358− 3363. (b) Patel, S.; Brooks, C. L. CHARMM fluctuating charge force field for proteins: I parameterization and application to bulk organic liquid simulations. J. Comput. Chem. 2004, 25 (1), 1−15. (20) Yuan, Y. N.; Mills, M. J. L.; Popelier, P. L. A.; Jensen, F. Comprehensive Analysis of Energy Minima of the 20 Natural Amino Acids. J. Phys. Chem. A 2014, 118 (36), 7876−7891.
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00296. Atomic charges and decomposition according to eq 10 (XLSX) Dipole moments and decomposition according to eq 16 for all 896 conformations (XLSX)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Frank Jensen: 0000-0002-4576-5838 Funding
This work was supported in part by the Danish Natural Science Research Council by Grant No. 4181-00030B and a grant from the Villum foundation. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) (a) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, UK, 1987. (b) Frenkel, D.; Smith, B. Understanding Molecular Simulations; Academic Press: San Diego, USA, 2002. (2) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D., Jr. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31 (4), 671− 690. (3) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem. 2003, 24 (16), 1999−2012. (4) Plazinski, W.; Lonardi, A.; Huenenberger, P. H. Revision of the GROMOS 56A6(CARBO) force field: Improving the description of F
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation (21) (a) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10 (44), 6615−6620. (b) Jensen, F. Unifying General and Segmented Contracted Basis Sets. Segmented Polarization Consistent Basis Sets. J. Chem. Theory Comput. 2014, 10 (3), 1074−1085. (22) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian-09; Gaussian Inc.: 2009.
G
DOI: 10.1021/acs.jctc.7b00296 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
