Evaluating the Effects of Geometry and Charge Flux in Force Field

Apr 26, 2018 - Atomic charges were obtained by fitting to results from density functional theory calculations using the HLY procedure and their geomet...
4 downloads 5 Views 1001KB Size
Download PDF
Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

pubs.acs.org/JPCA

Evaluating the Effects of Geometry and Charge Flux in Force Field Modeling Elaheh Sedghamiz*,†,‡ and Farhad Ghalami†,‡ †

Department of Chemistry, Aarhus University, Langelandsgade 140, Aarhus DK-8000, Denmark Department of Chemistry, University of Isfahan, Isfahan 81746-73441, Iran

‡

S Supporting Information *

ABSTRACT: We apply a model for analyzing the importance of conformational charge flux to 11 molecules with the R−(CH2)n−R structure (R = Cl, F, OH, SH, COOH, CONH2, and NH2 and n = 4−6). Atomic charges were obtained by fitting to results from density functional theory calculations using the HLY procedure, and their geometry dependence is decomposed into contributions from changes in bond lengths, bond angles, and torsional angles. The torsional degrees of freedom are the main contribution to the conformational dependence of atomic charges and molecular dipole moments, but indirect effects due to changes in bond distances and angles account for ∼15% of the variations. While the magnitude of charge flux and geometry effects have been found to be independent of the number of internal degrees of freedom, the nature of the R- group has a moderate influence. The indirect effects are comparable for all of the R-groups and are approximately one-half the magnitude of the corresponding effects in peptide models. However, the magnitudes are different, yet the relative importance of geometry and charge flux effects are completely similar to those of the peptide models, which suggests that modeling the charge flux effects for changes in bond lengths, bond angles, and torsional angles should be considered for developing improved force fields.

1. INTRODUCTION The development of accurate force fields is an important area in computational chemistry, which allows a precise description of large molecular systems. In contrast to electronic structurebased methods, force field methods use computational efficient empirical energy functions to allow one to simulate systems involving 103−106 atoms.1,2 The electrostatics energy is in production type force fields described as a Columbic interaction between fixed atomic point charges, typically obtained by fitting to the electrostatic potential determined from ab initio calculations on model compounds.3 Fixed charge force fields have provided impressive success in many applications, but the neglect of intra- and intermolecular ultimately limits the accuracy. They failed to incorporate electronic polarization at a fundamental level, which is especially important in simulations concerned with polar medium such as water. In addition, they possess low accuracy in modeling the dynamics of biological molecules such as proteins in diverse environments.4−7 The intramolecular polarization can in the lowest order be modeled as geometrydependent atomic charges, while intermolecular polarization is an essential component of solvation. Inclusion of polarization effects in force fields is currently pursued by three main approaches: the fluctuating charge,8,9 the Drude shell (dispersion oscillator),10 and induced atomic dipole11 models. The fluctuating charge model attempts to © XXXX American Chemical Society

describe the geometry dependence of the atomic charges by a second-order Taylor expansion of the energy in terms of atomic charges, but the presence of artificial long-range charge-transfer and unphysical scaling of the polarizability has prevented its widespread use. However, the QTPIE model, which has been developed to solve this problem, is able to provide a realistic description of in-plane polarizabilities.12 The fluctuating charge models have been employed successfully in simulations of DPPC−water monolayers, exploring free energetics underlying ion transport processes in biological ion channels and also transport of charged amino acid analogues across the water− bilayer interface.13−15 The induced dipole model implemented by assigning an isotropic polarizability tensor to atomic sites is a popular choice, and is employed, for example, by the AMOEBA16 force field. The Drude oscillator model can be considered as a simplified version of the induced dipole model, where dipole polarization is described by partitioning the atomic charge into two contributions connected by a spring. Because the usefulness of a model is ultimately decided by the quality of the results and the computational efficiency, it is important to select a consistent set of parameters. In the Dinur classification, the zeroth-order model considers only the atomic Received: December 11, 2017 Revised: April 20, 2018 Published: April 26, 2018 A

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

2. COMPUTATIONAL DETAILS Eleven molecular systems with the general structure R− (CH2)n−R, (R = NH2, OH, CONH2, F, Cl, COOH, and SH and n = 4−6) were chosen as models (Figure 1). The model systems hereafter would be identified by the number of carbons in their alkyl chain and their functional groups (CnR). An exhaustive conformational analysis was performed using the Merck Molecular Force Field (MMFF)30 using the ConfGen utility implemented in the Schrodinger program package.31All conformations were energy minimized at the wB97XD/augpcseg-1 level28,29 using the Gaussian 09 program package,32 and duplicate and symmetry related conformations were discarded. Atomic charges were determined by fitting to the electrostatic potential by the HLY procedure33 and have been constrained to reproduce the molecular dipole moment. The HLY procedure corresponds to a particular choice of sampling points for probing the electrostatic potential and differ from other alternatives by employ a larger molecular volume space sampling and a weighting factor to emphasize the sampling points in the medium distance range.34,35 It is important to mention that the ESP is relatively insensitive to the level of theory; therefore, it is unlikely that the general conclusions will change if other DFT or wave functional methods were used. Figure 2 shows the different parts of calculations that have been performed to provide the data for the model. The model is well described in the Theory.

charges, and they are forced to reproduce the dipole moment; the next level would include atomic dipoles and dipole fluxes forced to reproduce the molecular quadrupole moment. From an order-expansion point of view, it would be desirable to include charge polarization before dipole polarization, but an efficient parametrization of geometry dependence charges seems to be lacking in the literature.17 Intramolecular charge polarization effects can be introduced by allowing partial charges to change in response to conformational changes,18 which often is called a “charge flux” effect, and this is also essential for modeling, for example, intensities in infrared spectroscopy.19−23 Therefore, the charge flux could enter the force fields as parameters that determine the change of atomic charges with the change in the molecular internal coordinates. Currently, what is really needed is an understanding of the importance of the charge flux effect in a quantitative manner to implement it in force field calculations. Holt and Karlström suggested an induction correction model for including intramolecular polarization into the force fields. Their model describes the changes in the molecular charge distribution occurring as a response to changes of dihedral angles in the molecule based on a multicenter multipole distribution of the molecular charge distribution. Their findings showed that the induction correction model performs very well as compared to the fixed charged force fields.24−26 Molecular conformations are usually characterized by their torsional angles, but conformations also differ in terms of bond distances and angles. We have in a previous study analyzed ∼900 conformations of peptide models of the 20 neutral amino acids to assess the relative importance of geometry and charge flux effects, and partitioned these into contributions from bond length, bond angle, and torsional angle changes.27 The charge flux arising from the torsional degrees of freedom was found to be the major deficiency of fixed charge force fields, but indirect effects due to changes in bond angles and bond distances were ∼10% and ∼25% of the torsional effects, respectively. These findings were consistent for all of the 20 neutral amino acid peptide models, which perhaps are not surprising, as they all have polar amide groups as the common theme. The present study expands the analysis to a selection of systems where the polarity is more diverse to investigate whether the findings regarding the relative importance of geometry and charge flux effects for different degrees of freedom are generally valid. In the present work, we consider ∼500 conformations of 11 molecules each containing two groups of varying polarity as shown in Figure 1. As in our previous study, we work within a

Figure 2. Illustrative precedure of different sets of calculations performed in this study to provide the data for the model.

The increase in energy by replacing optimized bond lengths by conformationally averaged values is in all cases small (0.2 kJ/ mol on average, standard deviation = 0.19 kJ/mol, maximum = 1.8 kJ/mol), which reflects mean absolute bond length changes (mean absolute deviation (MAD)) of 0.0018 Å (standard deviation = 0.0019 Å, maximum = 0.017 Å). The energy increase by employing average bond angles is 3.7 kJ/mol on average with standard deviation = 4.1 kJ/mol and maximum = 41.2 kJ/mol. These values reflect mean absolute bond angle changes of 0.84° (standard deviation = 0.82°, maximum = 7.7°).

Figure 1. Model molecular systems considered in this study with the chemical formula CnH2nR2 (R shows different functional groups).

model where the molecular dipole moment is reproduced by atomic charges only.17 The analyses are done using a density functional theory method with an augmented polarized doubleζ type basis set,28,29 but the conclusions are likely not dependent on this specific choice.

3. THEORY The following briefly recapitulates the analysis presented in ref 27. The molecular dipole moment μk for a given conformation k is written in terms of atomic charges qik and positions Rik as shown in eq 1. B

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

indicates a contribution due to charge flux at a fixed geometry, and superscript gq indicates a small second-order term. Equation 12 shows the explicit definition of the first two terms on the right-hand side of eq 11, with similar definitions of the bond angle and torsional terms.19

Natom

μk =

∑ qik R ik(rk, θk , ωk )

(1)

i=1

The atomic positions Rik are determined by the bond distances rk, the bond angles θk, and the torsional angles ωk. Conformational-independent atomic charges are taken as the conformationally averaged values, and this provides the fixed charge approximation to the conformational-dependent dipole moments given by eq 2. μkq =

Δμkgr = Δμkqr =

qi R ik(rk , θk , ωk )

(3)

k=1

∑ ( qirθ − qi)R ikrθ( r ̅ , θ ̅ , ωk )

(13)

The major part of the conformational variation of the molecular dipole moment is correctly reproduced by fixed charge models, and this pure geometrical effect can be quantified by eq 14, where the last term is the conformational average length of the dipole moments using conformational average atomic charges.

Natom

∑ qikrθ R ikrθ( r ̅ , θ

Natom i=1

(4)

i=1

μkrθ =

Δμkqqω =

Natom

∑ qikr R ikr( r ̅ , θk , ωk )

(12)

The last term in eq 11 is a residual term arising from the chosen decomposition involving conformational-independent charge differences.

Nconf

∑ qik

∑ Δqikr R ikr( r ̅ , θk , ωk ) i=1

Atomic charges and molecular dipole moments are in addition determined for geometries where all bond distances are constrained to their conformationally averaged values, and where both bond distances and angle are constrained to their conformationally averaged values. μkr =

Natom

(2)

i=1

1 qi = Nconf

∑ (qikr − qi)ΔR ikr(rk, θk , ωk ) i=1

Natom

∑

Natom

Natom

̅, ωk )

Δμkgω

(5)

i=1

(6)

qikr = qikrθ + Δqikθ

(7)

qikrθ = qirθ + Δqikω

(8)

Natom

−

qirθ R ikrθ(r, θ , ωk )

(14)

4. RESULTS AND DISCUSSION Figure 3 shows the atomic charges and Δq values due to the changes in bond lengths, bond angles, and torsional angles for the non-hydrogen atoms in a representative conformation of HO−(CH2)4−OH. The torsional degrees of freedom provide the largest variation, ranging between 7% and 44% of the absolute value for the different atoms. The torsional dependence of the charges is on average 24% of the absolute values, while the variation due to the bond angles is ∼3% and the variations due to the bond distances are almost negligible. Table 1 shows the non-hydrogen atom and conformationally averaged values for the absolute charge differences in eq 6 for the 11 systems in Figure 1. The torsional variation (Δqω) is the largest by 100% (as measured relative to the Δqrθω value), which in magnitude is identical to the total variation. The indirect effects of bond distance and bond angle changes with conformation are 2% and 17%, respectively, while the conformational independent term is 8%. The summation to more than 100% reflects that the charge variations are average over absolute values, and thus implies that some of the changes are in the opposite direction.

(9)

R ik(rk , θk , ωk ) = R ikrθ( r ̅ , θ ̅, ωk ) + ΔR ikθ ( r ̅ , θk , ωk ) (10)

Equations 9 and 10 allow a decomposition of difference in the dipole moments from eqs 1 and 2 as shown in eq 11, and this is the contribution not included in fixed charge force fields. Δμkrθω = μk − μkq = Δμkgr + Δμkqr + Δμkgqr + Δμkgθ + Δμkqθ + Δμkgqθ + Δμkqω + Δμkqqω

qirθ R ikrθ( r ̅ , θ ̅ , ωk )

The decomposition in eq 11 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. As for the atomic charges, we will quantify each contribution in terms of the magnitude, and the values should thus only be taken as an indication of the relative importance. The Δμgqr and Δμgqθ second-order terms are always significantly smaller (few percent) than the first-order terms (Δμgr, Δμqr, Δμgθ, and Δμqθ) and will not be discussed.

Because charge differences can be either positive or negative, we quantify the relative importance of each of the terms in eq 9 as an average of absolute values. The molecular geometry can similarly be written as a sum of terms involving only torsional, bond angle, and bond distance contributions. + ΔR ikr(rk , θk , ωk )

∑ i=1

Δqikrθω = qik − qi = Δqikr + Δqikθ + Δqikω + Δqiqq Δqiqq = ( qirθ − qi)

∑ i=1

This allows the conformational dependence of the atomic charges to be decomposed into terms arising from changes in bond distance, bond angles, and torsional angles, and a residual conformational-independent term according to eqs 6−8.

qik = qikr + Δqikr

=

(11)

The superscript g indicates a contribution due to a geometric displacement of fixed atomic charges, while superscript q C

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the largest effect and has values very close to those of Δμrθω. This is a term that is already captured by fixed charged force fields and shows their ability in considering the conformational variation of the molecular dipole moments. The μ−μ q difference represents the part that is neglected in fixed charge force fields, and this is on average ∼16% of the reference dipole moment, with most values in the 6−31% range (Table S3). The average magnitude of the dipole moment is 1.3−3.3 D for all of the systems, but there is a large conformational dependence. The dipole moment for the 36 conformations of C4CONH2 functional group, for example, varies between 0.0 and 7.8 D, and this is reflected in the Δμgω value of 2.17 D. The bond distance effect (Δμgr + Δμqr) shown in Table 3 is ∼3% of the total effect (μ−μq), the bond angle effect (Δμgθ + Δμqθ) is ∼10%, and the torsional charge flux term (Δμqω) is of magnitude almost equal to that of the total effect. This again reflects that some compensation due to sign differences takes place. Comparing charge flux and geometry effect for each degree of freedom, the charge flux effect (Δμqr) is 1 order of magnitude larger than the geometry effect (Δμgr) for the bond distance terms, while the charge flux effect (Δμqθ) dominates over the geometry effect (Δμgθ) for the bond angle terms by a factor of 2. For the torsional degrees of freedom, the geometry effect (Δμgω) is about a factor of 4 more important than the torsional charge flux term (Δμqω), which is clear from Table 3. Table 3 (and also Table S2) show that there are no significant changes of different effects by varying the length of the alkyl chain. This most likely reflects that the charge flux and geometry effects do not depend on the number of internal degrees of freedom of a molecular system. From Table 3 and Figure 4, the charge flux and geometry effects, for example, for C4NH2, C5NH2, and C6NH2 molecular systems are very close in magnitude (and also in the percentage of contribution to the total dipole variation (Table S2)) but differ slightly from the corresponding values for C4CONH2. It indicates that the geometry and charge flux effects for different degrees of freedom depend slightly on the nature of the R-group. The importance of the indirect (bond distance and angle) effects relative to the direct (torsional) effect is roughly a factor of 2 smaller for the present systems as compared to the peptide models, and this is also reflected in the magnitudes of the geometry changes (MADbond = 0.0018 Å and MADangle = 0.84° for the present systems as compared to the values for the

Figure 3. Atomic charges for the non-hydrogen atoms in a representative energy conformation for the HO−(CH2)4−OH and charge differences defined in eqs 6−8. The values represented on the atoms are partial charges, and the values in parentheses are Δq values due to the changes in bond lengths, bond angles, and torsional angles, respectively.

There is no significant change in the charge variations with increasing length of the alkyl chain in the model molecules (Table S1). For example, the charge variations in the C6NH2 system, which has the largest set of conformations and thus the largest number of internal degrees of freedom, is similar to those of the C4NH2 system. There is little dependence of the magnitude of the charge variations with the nature of the Rgroup, and substituents such as F and Cl display smaller changes than, for example, SH and CONH2. The averaged values of charge variations in Table 1 are roughly one-half of the values found for the peptide models, but the relative importance is almost identical.27 Table 2 shows the Cartesian dipole components for the decomposition according to eq 8 for the same representative conformation for the HO−(CH2)4−OH model as in Figure 2. The charge flux contribution from the torsional angles (Δμqω) dominates and is slightly larger than the total effect. The charge flux contribution from the bond angles (Δμqθ) is the second most important, while the charge flux effect due to bond distances (Δμqr) is the smallest. The averaged dipole contributions over all conformations are shown in Table 3. The Δμgω term (geometry effects due to torsional angles) has

Table 1. Non-Hydrogen Atom and Conformationally Averaged Values for the Absolute Charge Differences in Equations 6−8 for the Models in Figure 1a

a

molecular system

Nconf

Δqrθω

Δqr

Δqθ

Δqω

Δqqq

C4Cl C4F C4CONH2 C4COOH C4OH C4SH C4NH2 C5Cl C6Cl C5NH2 C6NH2 average %

11 10 36 31 65 23 27 23 74 78 150 528

0.0315 0.0229 0.0312 0.0267 0.0400 0.0701 0.0496 0.0377 0.0419 0.0589 0.0596 0.0427 100

0.0006 0.0007 0.0014 0.0009 0.0009 0.0008 0.0009 0.0007 0.0007 0.0010 0.0009 0.0008 2

0.0048 0.0029 0.0051 0.0046 0.0077 0.0117 0.0096 0.0062 0.0070 0.0102 0.0094 0.0072 17

0.0331 0.0237 0.0309 0.0280 0.0395 0.0698 0.0509 0.0371 0.0404 0.0572 0.0583 0.0426 100

0.0029 0.0018 0.0013 0.0014 0.0034 0.0047 0.0057 0.0035 0.0031 0.0055 0.0049 0.0035 8

Superscripts r, θ, and ω indicate changes due to changes in bond lengths, bond angles, and torsional angles, respectively. D

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 2. Changes in the Cartesian Molecular Dipole Moment (debye) for a Representative Conformation of HO−(CH2)4−OHa component

μ

μ−μq

Δμgr

Δμqr

Δμgθ

Δμqθ

Δμqω

Δμqqω

X Y Z length %

−0.866 0.960 1.732 2.161

−0.004 −0.355 −0.532 0.640 100

0.002 0.001 0.004 0.005 1

−0.008 0.004 −0.005 0.011 2

0.006 0.002 0.012 0.014 2

0.020 0.003 0.012 0.023 4

−0.033 −0.357 −0.547 0.654 102

0.009 −0.007 −0.007 0.014 2

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 12−14. a

Table 3. Changes in the Magnitude of Molecular Dipole Moments (debye) Averaged over All Conformations of the 11 Modelsa molecular system

Nconf

μ

Δμrθω

μ−μq

Δμgω

Δμgr

Δμqr

Δμgθ

Δμqθ

Δμqω

Δμqqω

C4Cl C4F C4CONH2 C4COOH C4OH C4SH C4NH2 C5Cl C6Cl C5NH2 C6NH2 average %

11 10 36 31 65 23 27 23 74 78 150 528

1.965 1.824 3.283 2.209 2.097 1.882 1.277 2.689 2.549 1.649 1.631 2.096

1.368 1.307 2.120 0.880 0.725 0.936 1.002 0.639 1.144 0.368 0.669 1.014

0.130 0.168 0.207 0.145 0.527 0.456 0.405 0.179 0.188 0.434 0.502 0.304 100

1.447 1.361 2.169 0.830 0.888 1.034 1.128 0.704 1.226 0.498 0.919 1.109

0.000 0.000 0.001 0.001 0.001 0.001 0.002 0.001 0.001 0.002 0.002 0.001 0.4

0.005 0.006 0.018 0.007 0.008 0.003 0.006 0.007 0.006 0.007 0.007 0.007 2.4

0.006 0.003 0.006 0.009 0.010 0.011 0.015 0.005 0.008 0.013 0.016 0.009 3.1

0.011 0.007 0.041 0.018 0.024 0.023 0.024 0.011 0.012 0.027 0.026 0.020 6.7

0.104 0.163 0.225 0.147 0.528 0.477 0.388 0.156 0.170 0.431 0.492 0.298 98.3

0.025 0.015 0.018 0.015 0.011 0.022 0.017 0.023 0.025 0.019 0.019 0.019 6.2

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 eq 11.

Figure 4. Contribution of geometry and charge flux effects in percentage of the total variation for different degrees of freedom as a function of the nature of the R-group (the values are reported in Table S2).

peptide models of MADbond = 0.0024 Å and MADangle = 1.2°).27 Amino acids have diverse functional groups in their molecular structure and thus have larger dipole moments as compared to the studied molecular systems. Therefore, they possess more intramolecular interactions resulting in more conformational dependence on bond distance and angles.

in different degrees of freedom for 11 molecular systems R− (CH2)n−R (with R = Cl, F, OH, SH, COOH, CONH2, NH2 and n = 4−6) and compared them to previous results for peptide models. The charge flux effect dominates over geometry effects for changes in bond distances and bond angles, while the geometry effect is the dominant effect for torsional degrees of freedom. The average difference between the dipole moment of different conformations and the dipole moment calculated by fixed atomic charges is ∼16% of the total dipole moment. This effect, which is not represented in fixed

5. CONCLUSION On the basis of density functional calculations, we have evaluated the relative importance of geometry and charge flux E

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A charge force fields, is very similar in magnitude to the charge flux due to torsional angles only, but indirect effects due to changes in bonds and angles are ∼15% of the total variation. This reflects that different degrees of freedom in molecules provide partly compensating effects. Furthermore, it has been shown that the nature of the R group for a molecular system is more important than the number of internal degrees of freedom in determining the magnitude of charge flux and geometry effects. However, the relative importance of charge flux and geometry effects is the same for all studied systems and very similar to the previous findings for peptide models, which suggest that the findings are general. Therefore, it is suggested that the improved force fields should take into account the geometry dependence of the molecular dipole moment, which could be reached by including the charge flux effects that are able to model the dependence on both torsional degrees of freedom and changes in bond angles. The indirect effects of bond distances have a smaller effect, and the importance of their implications could remain a choice that depends on the molecular systems.

■

(6) Simonson, T.; Brooks, C. L. Charge Screening and the Dielectric Constant of Proteins: Insights from Molecular Dynamics. J. Am. Chem. Soc. 1996, 118, 8452−8458. (7) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, J. A. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 27, 4983−5013. (8) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358−3363. (9) Patel, S.; Brooks, C. L. Empirical force fields for biological macromolecules: Overview and issues. J. Comput. Chem. 2004, 25, 1− 15. (10) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D. Force Field for Peptides and Proteins based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430− 5449. (11) 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. J. Comput. Chem. 2001, 22, 1048−1057. (12) Chen, J.; Martínez, T. J.; QTPIE. Charge Transfer with Polarization Current Equalization. A Fluctuating Charge Model with Correct Asymptotics. Chem. Phys. Lett. 2007, 20, 315−20. (13) Lee, W. G.; Davis, J. E.; Patel, S. Origin and Control of Superliner Polarizability Scaling in Chemical Potential Equalization Methods. J. Chem. Phys. 2008, 14, 144110. (14) Patel, S.; Davis, J. E.; Bauer, B. A. Exploring Ion Permeation Energetics in Gramicidin A Using Polarizable Charge Equilibration Force Fields. J. Am. Chem. Soc. 2009, 10, 13890−1. (15) Lucas, T. R.; Bauer, B. A.; Patel, S. Charge Equilibration Force Fields for Molecular Dynamics Simulations of Lipids, Bilayers, and Integral Membrane Protein Systems. Biochim. Biophys. Acta, Biomembr. 2012, 1, 318−29. (16) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. The Polarizable Atomic Multipole-based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046−4063. (17) Dinur, U. Atomic Multipoles and Perpendicular Electrostatic Forces in Diatomic and Planar Molecules. J. Comput. Chem. 1991, 12, 91−105. (18) 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. (19) Person, W. B.; Zerbi, G. Vibrational Intensitites in Infrared and Raman Spectrocopy; Elsevier Scientific Publishing: Amsterdam, 1982. (20) Torii, H. Intra- and Intermolecular Charge Fluxes Induced by the OH Stretching Mode of Water and Their Effects on the Infrared Intensities and Intermolecular Vibrational Coupling. J. Phys. Chem. B 2010, 114, 13403−13409. (21) Galimberti, D.; Milani, A.; Castiglioni, C. Charge Mobility in Molecules: Charge fluxes from Second Derivatives of the Molecular Dipole. J. Chem. Phys. 2013, 138, 164115. (22) Milani, A.; Galimberti, D.; Castiglioni, C.; Zerbi, G. Molecular Charge Distribution and Charge Fluxes from Atomic Polar Tensors: The Case of OH Bonds. J. Mol. Struct. 2010, 976, 342−349. (23) Milani, A.; Castiglioni, C. Modeling of Molecular Charge Distribution on the Basis of Experimental Infrared Intensities and First-Principles Calculations: The Case of CH Bonds. J. Phys. Chem. A 2010, 114, 624−632. (24) Holt, A.; Karlstrom, G. An Intramolecular Induction Correction Model of the Molecular Dipole Moment. J. Comput. Chem. 2008, 29, 1084−1091. (25) Holt, A.; Karlstrom, G. Inclusion of the Quadrupole Moment when Describing Polarization. The Effect of the Dipole-Quadrupole Polarizability. J. Comput. Chem. 2008, 29, 2033−2038. (26) Holt, A.; Karlstrom, G. Inclusion of the Quadrupole Moment when Describing Polarization. The Effect of the Dipole-Quadrupole Polarizability. J. Comput. Chem. 2008, 29, 2485−2486.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.7b12198. Averaged atomic charge and dipole moment variations in percentage of the total variations for 11 molecular systems (PDF) Atomic charges and decomposition according to eqs 6−8 for all 528 conformations (XLSX) Dipole moments and decomposition according to eq 11 for all 528 conformations (XLSX)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected], elahe.sedghamiz@gmail. com. ORCID

Elaheh Sedghamiz: 0000-0003-3998-9701 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported in part by the Danish Natural Science Research Council by grant no. 4181-00030B. We would like to thank Prof. Frank Jensen for his guidance.

■

REFERENCES

(1) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, UK, 1987. (2) Frenkel, D.; Smith, B. Understanding Molecular Simulations; Academic Press: San Diego, CA, 2002. (3) 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, 10269−10280. (4) Vierros, S.; Sammalkorpi, M. Phosphatidylcholine Reverse Micelles on the Wrong Track in Molecular Dynamics Simulations of Phospholipids in an Organic Solvent. J. Chem. Phys. 2015, 142, 142. (5) Vorobyov, I.; Allen, T. W. The Electrostatics of Solvent and Membrane Interfaces and the Role of Electronic Polarizability. J. Chem. Phys. 2010, 132, 132. F

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (27) Sedghamiz, E.; Nagy, B.; Jensen, F. Probing the Importance of Charge Flux in Force Field Modeling. J. Chem. Theory Comput. 2017, 13, 3715−3721. (28) Chai, J.-D.; Head-Gordon, M. Long-range Corrected Hybrid Density Functionals with Damped Atom−Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (29) Jensen, F. Unifying General and Segmented Contracted Basis Sets. Segmented Polarization Consistent Basis Sets. J. Chem. Theory Comput. 2014, 10, 1074−1085. (30) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and Performance of MMFF94. J. Comput. Chem. 1996, 17, 490−519. (31) Watts, K. S.; Dalal, P.; Murphy, R. B.; Sherman, W.; Friesner, R. A.; Shelley, J. C. ConfGen: A Conformational Search Method for Efficient Generation of Bioactive Conformers. J. Chem. Inf. Model. 2010, 50, 534−546. (32) 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.; et al. Gaussian 09; Gaussian, Inc.: Pittsburgh, PA, 2009. (33) Hu, H.; Lu, Z.; Yang, W. Fitting Molecular Electrostatic Potentials from Quantum Mechanical Calculations. J. Chem. Theory Comput. 2007, 3, 1004−1013. (34) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431−439. (35) Breneman, C. M.; Wiberg, K. B. Determining Atom-Centered Monopoles from Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11, 361−373.

G

DOI: 10.1021/acs.jpca.7b12198 J. Phys. Chem. A XXXX, XXX, XXX−XXX