Theoretically grounded approaches to account for polarization effects in fixed-charge force fields

Molecular polarization describes the process by which a molecule’s electronic structure responds to the presence of an external electric field, which can be (and, in the context of this paper, is) caused by the surrounding molecules in a condensed phase.^1,2 This electric field causes a distortion in the geometry of the molecule, i.e., in the relative position of its atomic nuclei, and also in the electronic cloud from its ground state. These rearrangements increase the dipole moment of the molecule relative to the gas phase value, with the difference between the two being called the induced dipole moment. Such distortions of the molecular wavefunction carry an energy penalty, alternately referred to as the distortion energy, polarization cost, or self-energy of polarization.³ However, the polarized molecule interacts more strongly with the surrounding molecules than if it had remained in the ground state, leading to a favorable stabilization energy.⁴ The latter is larger in magnitude than the distortion energy, leading to an overall lowering of the free energy of the system. While this general picture has been well understood for many decades,^1,5–7 its implications for molecular modeling are still the subject of heated scientific debate.

Polarization is inherently a many-body cooperative effect—each molecule contributes toward the electric field acting on other molecules, and is in turn influenced by that electric field—which makes it challenging to describe in classical molecular models. Several approaches have been proposed over the years to explicitly include polarization effects and develop what are called polarizable force fields; these are not covered here, and the reader is directed to prior reviews on the topic.^8–11 Despite such advances, however, polarizable force fields are still harder to parameterize (e.g., they have considerably more parameters) and are much more computationally expensive than their non-polarizable counterparts, so the latter remain the “weapon of choice” for most current applications of molecular simulations.

It is important to emphasize that the term “non-polarizable” simply means that the electrostatic degrees of freedom of the molecule, normally described by a set of partial charges interacting through Coulomb’s law, remain fixed and do not respond to changes in the surrounding environment. It does not necessarily imply that polarization interactions are completely neglected. In fact, polarization and other many-body effects are most often accounted for indirectly by assigning effective pairwise potential parameters. For example, in the early development of the OPLS-UA force field for alcohols, point charges were selected so that the model dipole moment was ∼25% larger than the corresponding gas phase dipole moment,¹² a practice which, notwithstanding some variations in implementation, has persisted to this day.

Because they provide a good balance between accuracy and computational expense, generic non-polarizable force fields that cover a wide variety of compounds (e.g., CGenFF,¹³ GAFF,¹⁴ OPLS-AA,¹⁵ TraPPE,¹⁶ and OpenFF¹⁷) have found use in a variety of fields and have been at the forefront of many scientific advances over the past 30 years or so. However, as both computer power and data availability increase, there are signs that generic fixed-charge force fields are approaching their limits. For example, large benchmark studies have revealed systematic deviations between experimental data and predictions from non-polarizable force fields for key physical properties such as the dielectric constant and the surface tension.^18–20 In addition, despite their much higher computational expense and predictive potential, molecular models have not yet surpassed more empirical approaches, such as group-contribution or implicit solvent models, in solvation free energy predictions.^21,22 These observations suggest that there is some essential missing physics in non-polarizable models, which arises, at least in part, from an improper consideration of polarization effects.

This Perspective argues for an intermediate approach between explicitly polarizable models and conventional non-polarizable force fields, whereby the computational expedience of the latter is kept, but polarization effects are rigorously accounted for in a theoretically grounded manner. Efforts to develop a comprehensive theoretical framework for describing electronic polarization in non-polarizable simulations began with the ground-breaking work of Leontyev and co-workers on their Molecular Dynamics in Electronic Continuum (MDEC) model,^23–27 which introduced key concepts such as screened Coulomb interactions, the effective dipole moment, and the scaled ion charge, and proposed polarization corrections for the solvation free energy and the dielectric constant. This model laid the foundations for a significant body of work that followed from many research groups, including the development of an alternative approach for estimating the optimal effective dipole moment of non-polarizable models, which we will call the “halfway-charge” approach.^28,29 Based on these ideas, our group has recently proposed the Polarization-Consistent Approach (PolCA) for force field development and applied it to build improved models for alkanes,³⁰ alcohols,³¹ ketones,³² and organosilicates³³ in the liquid phase. The core principles mentioned above are also at the basis of the so-called scaled charge models of ionic liquids and electrolyte solutions,^34–38 as well as recent approaches for force field parameterization, such as IPolQ,²⁹ QUBE,^39,40 and RESP2,⁴¹ albeit with some rather important differences. The remainder of this paper starts with a historical background on approaches to implicitly account for polarization effects, after which the fundamental principles behind PolCA are described and compared with other similar strategies. Then, several important properties/systems are examined in turn, discussing if (and how) polarization effects can be effectively accounted for using this kind of approach. This discussion focuses mainly on force fields for neutral molecules but also includes a brief reflection on ionic liquids and electrolytes. This paper ends with an outlook on the future of this exciting area of research.

It was recognized very early on in the development of non-polarizable force fields for liquids, during the mid-1980s, that point charges assigned to a given molecule should yield dipole moments that are larger than those of an isolated molecule, so as to at least partly account for polarization effects. One option, adopted by the OPLS developers,¹² was to treat the charges as empirical fitting parameters, but subject to the constraint that they would increase the dipole moment of the model by ∼25% over the gas-phase experimental value. Other force fields, e.g., AMBER,⁴² preferred to avoid empirically fitting point charges, instead assigning them on the basis of gas-phase quantum mechanical (QM) calculations. It was found that the Hartree–Fock (HF) method^43,44 with the 6-31G* basis set⁴⁵ systematically over-polarized molecules and was therefore seen as a good choice to approximately account for implicit polarization in the liquid phase.⁴⁶ Recent versions of the most widely employed general force fields still make use of this idea, e.g., OPLS3⁴⁷ uses CM1A-BCC charges,⁴⁸ which are themselves calibrated against HF/6-31G* charges; GAFF¹⁴ uses RESP charges⁴⁶ fitted to HF/6-31G* calculations; and CGenFF¹³ uses MK charges⁴⁹ fit to MP2/6-31G*.

In a recent systematic study, Zhou et al.⁵⁰ showed that HF/6-31G* charges generally underestimate liquid phase polarization and, more importantly, lead to inconsistent levels of polarization across different types of molecules. The latter inconsistency was confirmed in recent systematic benchmark studies of ab initio methods for the calculation of dipole moments.^51,52 This suggests that the success of force fields based on HF/6-31G* charges (or similar) is largely due to error cancellation, since any shortcomings in the electrostatic interactions can be compensated by empirical fitting of other non-bonded parameters, such as repulsion/dispersion parameters. This increases the interdependency between model parameters and ultimately leads to poor transferability when moving away from the systems and/or conditions used during parameterization.

This is essentially the same as Eq. (3) in Karamertzanis et al.,²⁸ but where only dipole moments are considered and the summations over different sites on each molecule A and B have been omitted for simplicity. The first term on the right-hand side of the above equation represents the ground state electrostatic energy between all molecules in their unpolarized state. The second and third terms added together represent the induction or polarization energy arising from the interaction between polarized molecules—in this case, represented by halfway-polarized dipole moments. The fourth and final term is a residual energy that was neglected by Karamertzanis et al.,²⁸ but which will be revisited later in this paper.

The effectiveness of the halfway-charge approach relies on the accuracy of the linear response approximation. Cerutti et al.²⁹ demonstrated its validity for a series of simple examples composed of interacting point dipoles. In a series of papers, Orozco and co-workers^4,54,55 provided numerical evidence in support of the above approximation by showing, based on quantum mechanical (QM) calculations with a dielectric continuum solvent for a variety of systems, that the magnitude of the stabilization energy is indeed approximately twice that of the distortion energy. Jorge et al. have independently confirmed this for water based on a computationally more accurate QM/MM approach.⁵⁶ It is important to note, however, that Eq. (4) only captures the average polarization energy of the system in a mean-field sense; i.e., while it represents a useful approximation in the context of non-polarizable models, it cannot be expected to accurately describe instantaneous fluctuations in the dipole moments for real liquids.²⁹ 

The halfway charge idea was subsequently adopted by Cerutti et al.²⁹ when developing a point charge assignment scheme for a new version of the AMBER force field. In their scheme, called IPolQ, the model charges for a particular molecule are calculated as the average between two QM calculations at the MP2/cc-pV(T + d)Z level of theory—one for the isolated solute molecule in vacuum and the other for the solute molecule surrounded by an explicit water solvent (i.e., a QM/MM calculation). More precisely, they described the electrostatic field acting on the solute molecule by a set of explicit point charges obtained from MD simulation snapshots, but with their magnitude adjusted to correspond to a “polarized version” of the TIP4P-Ew water model⁵⁷—i.e., they increased the value of the charges so that the dipole moment of the surrounding water molecules was greater than that of normal TIP4P-Ew by the same degree that TIP4P-Ew exceeds the dipole moment of water vapor.²⁹ As will be discussed later, this attempt to represent the polarization state of the surrounding liquid in a realistic way is essential to the success and self-consistency of the halfway charge approach.⁵⁸ 

More recently, several other studies proposed modified versions of the halfway charge approach. For example, the IPolQ-Mod scheme⁵⁹ used the same idea but replaced the computationally intensive QM/MM calculation of the liquid state with a much simpler and faster implicit solvent calculation with a dielectric constant equal to that of liquid water. The underlying assumption is that the dielectric continuum model is able to provide an accurate representation of the real polarization environment of the liquid state, which turns out to be a rather crude approximation.^56,60–63 In fact, it was shown that hydration free energy predictions obtained with IPolQ-Mod charges were inconsistent with those of the reference potential method, which applies a more rigorous explicit representation of the surrounding solvent,⁶⁴ and do not seem to lead to significant improvements over “standard” HF/6-31G*-based charges.^59,65,66

Similarly, the point charge assignment in the QUantum mechanical BEspoke (QUBE) force field developed by Cole and co-workers^39,40 is based on QM calculations of the target molecule in a dielectric continuum model with ε = 4. This particular value led to model dipole moments that were, on average, approximately halfway between those of the gas and the liquid, where the latter was described by a dielectric continuum model with ε = 80 to represent water. This idea is also at the heart of the recently proposed RESP2 method of force field charge assignment,⁴¹ where point charges are computed as linear combinations of gas- and liquid-phase charges, again with the reference liquid state being described by an implicit solvent model representing water. Interestingly, the scaling factor was treated as an adjustable parameter whose optimal value was found to be ∼0.6—i.e., somewhat larger than the ½ value derived by Karamertzanis et al.²⁸ 

A common misconception when applying Eq. (6) is to assume, as Berendsen et al. did,⁶⁷ that the non-polarizable model accurately describes the level of polarization of the real liquid, i.e., μ_L = μ_M. This approach can be represented by a thermodynamic diagram,⁶⁹ as shown in Fig. 1(a), where the distortion energy coming from Eq. (6) is added to the MD prediction for phase change (free) energies before comparison with the corresponding experimental property. This assumption was challenged by Leontyev and Stuchebrukhov in a series of papers in which they put forth their MDEC model^24–27 [also called Electronic Continuum Correction (ECC) scheme]. The MDEC model is based on the realization that classical fixed-charge force fields can only account for the nuclear polarization response, i.e., arising from changes in relative positions and orientations of the atomic nuclei, but cannot describe purely electronic polarization effects, i.e., those that arise from fluctuations of the electronic clouds of the molecules. The electronic response of a given fluid is typically much faster than the nuclear response and is characterized by its infinite-frequency dielectric constant, ε_∞, which can be approximated as the square of the refractive index.²⁴ This means that the interaction energy between the polarized molecules and the surrounding electronic medium in a condensed phase (E_Elec) also needs to be added to predictions from MD simulations, since such interactions take place in experimental systems but not in non-polarizable simulations. This is represented by an additional step in the thermodynamic cycle of Fig. 1(b).

An even more important implication of the MDEC theory is that the effective dipole moment of a non-polarizable model is scaled down from the real liquid dipole moment by a factor of $1/ε∞$ due to dielectric screening by the electronic degrees of freedom of the solvent/liquid. Since the value of ε_∞ is close to 2 for most liquids, this corresponds to a scaling factor of ∼0.7. In other words, an MD simulation with fully polarized molecules surrounded by a dielectric continuum of ε_∞ is equivalent to a simulation carried out in vacuum but with a scaled-down dipole moment or charge. This equivalence is represented on the right-hand side of Fig. 1(b). Note that while the vast majority of simulations with non-polarizable models adopt the approach depicted in the right-most arrow in Fig. 1(b) (i.e., simulations are run in vacuum), some authors have attempted to account for electronic polarization effects by running simulations with fully polarized molecules within an explicit dielectric continuum.⁷⁰ 

This expression, used in the PolCA framework, provides a way to estimate the dipole moment (and, indirectly, the point charges) of a non-polarizable model that best incorporates the energetic contributions of polarization in an implicit way.³² It has recently been demonstrated that this approach indeed yields model charges/dipoles that provide an optimal description of the potential energy surface of the liquid phase.⁵⁸ It predicts a model dipole moment of 2.262 D for water, which is very close to that of the best performing fixed-charge water models (e.g., TIP4P/2005 has a dipole moment of 2.305 D⁷⁷), and 2.108 D for methanol, which matches very closely the empirically fitted value of 2.07 D.³¹ 

A second point in common between MDEC theory and the halfway charge theory is the prediction that the use of scaled-charge models leads to net polarization corrections to the energy that are very small or zero. For both approaches, this can be represented by thermodynamic cycles, as depicted in Fig. 1. The cycle for the MDEC theory, shown in Fig. 1(b), has been discussed in Sec. II B. By applying this theory to TIP4P water, Leontyev and Stuchebrukhov have shown that indeed the electronic and distortion energy contributions, estimated from Eqs. (6) and (7), nearly canceled each other out.²⁵ More recently, Milne and Jorge demonstrated that the difference between both components, and hence the net polarization correction, was less than 1 kJ/mol if a more accurate approach to estimate E_dist and E_pol, based on expansions up to the quadrupole moment, is applied.⁶⁹ Barrera and Jorge, when developing the PolCA model for alcohol molecules,³¹ estimated the net polarization correction for several solute/solvent pairs using Eqs. (6) and (7), examples of which are shown in Table I. For non-polar molecules in polar or non-polar solvents (first four rows), the result is trivial: both contributions are approximately zero as there are no meaningful solute polarization effects in such systems. For polar compounds in polar solvents (next five rows), as was the case for water,²⁵ both contributions are large in magnitude but opposite in sign, nearly canceling out and yielding net polarization energies that are very close to zero. Notably, however, when polar molecules are dissolved in non-polar solvents (last four rows), the electronic contribution dominates and the net polarization correction is non-negligible. This result agrees with the conjectures of Leontyev and Stuchebrukhov²⁴ and has important implications for the prediction of solvation free energies (see Sec. III C).

The thermodynamic cycle for the halfway-charge theory, as proposed by Cole et al.³⁹ based on ideas put forth by Chipot,⁷⁸ is depicted in Fig. 1(c). Here, the first step involves polarizing the target molecule up to an intermediate level of polarization (i.e., μ_M), followed by an MD calculation of the (free) energy using the non-polarizable model with intermediate charges, while the third and final step comprises polarizing the molecule all the way to the full liquid state polarization level (i.e., μ_L). This final step includes two contributions, one being the remainder of the distortion energy and the other being the favorable stabilization energy arising from interaction with the surrounding liquid. As described in Sec. II A, the theory predicts that the energies arising from the first and third steps cancel out when the dipole moment of the model has the correct intermediate degree of polarization.

An important consistency test for this theory has recently been performed by Jorge et al.,⁵⁸ who estimated the favorable and unfavorable contributions to the polarization energy using a procedure adapted from Cole et al.³⁹ They showed that the predicted enthalpy of vaporization of both water and methanol was independent of the choice of model dipole moment when the corresponding polarization corrections were applied—hence, the approach was self-consistent. However, this was only the case when the correct dipole moment of the reference liquid state [i.e., the end state of step 3 of the cycle shown in Fig. 1(c)], calculated using a high-accuracy QM/MM method,⁶⁰ was used. In contrast, when the liquid phase dipole moment was estimated from a dielectric continuum model, as was done in the original QUBE, IPolQ-Mod, and RESP2 methods, the consistency test was not met. This is because, as will be discussed in more detail in Sec. III A, dielectric continuum methods are unable to fully capture the correct degree of polarization for hydrogen-bonding liquids. Incidentally, this underestimation of the liquid dipole moment provides a possible explanation for the fact that the optimal scaling factor for RESP2 point charges was found to be between 0.5 and 0.7, i.e., somewhat larger than the factor 0.45 predicted by Eq. (16).

The recently proposed PolCA force field parameterization strategy takes advantage of the similarities between the MDEC and halfway charge theories to propose a new workflow for parameterizing non-polarizable force fields so as to provide a theoretically consistent consideration of implicit polarization effects. This works as follows:³² 

1. Calculate the liquid phase dipole moment of the molecule of interest. The accuracy of this step is of paramount importance, and a method that is able to describe both local and mean-field polarization contributions is essential, as will be discussed in more detail in Sec. III A.

2. Apply Eq. (16) to calculate the model dipole moment (μ_M) that leads to net zero polarization correction energy for the pure fluid.

3. Calculate initial point charges on the molecule using an appropriate method to partition the electron density into atomic contributions (e.g., the DDEC method⁷⁹ is a suitable choice).

4. Scale the above charges uniformly, so the dipole moment matches μ_M. This is the simplest approximation, although other options are possible (e.g., scaling the charge on each atom proportionally to its electronegativity) as long as the total dipole moment is constrained to be equal to μ_M.

5. Fit the LJ parameters (or a subset thereof) to a selected set of condensed-phase properties (typically density, enthalpy of vaporization, and self-diffusion).

6. Polarization energy corrections are NOT applied to pure liquid phase-change properties (e.g., enthalpy of vaporization and self-solvation free energy).

7. Apply polarization corrections to solvation free energies in different solvents, estimated, for example, from Eqs. (6) and (7). This is because, as discussed above, the favorable and unfavorable contributions do not generally cancel out in such cases.

8. Apply polarization corrections to the dielectric constant predictions. The precise functional form of this correction will be discussed in Sec. III B.

The PolCA framework has so far been applied to aliphatic hydrocarbons,³⁰ alcohols,³¹ ketones,³² and organosilicates,³³ yielding significant improvements over state-of-the-art force fields. This is particularly the case for predictions of the dielectric constant and solvation free energies in non-polar solvents, where polarization corrections are required. In Sec. III, the impact of rigorously accounting for polarization effects on several properties is critically discussed.

Non-polarizable force fields are typically parameterized to reproduce the potential energy surface of the target compounds. This means that bulk phase properties that depend solely on the PES (e.g., density and self-diffusion) are normally well reproduced by applying the standard parameterization approach—i.e., obtaining effective interaction parameters that account, at least in some approximate way, for bulk polarization effects. However, properties that involve a change of polarization environment (e.g., phase-change energies and interfacial properties) or that depend explicitly on the dipole moment surface (e.g., dipole moment and dielectric constant) require a more theoretically grounded approach to eliminate systematic deviations arising from a neglect of polarization effects. In some cases, the implicit approaches described in Sec. II are able to yield remarkable improvements in predictions, while some cases remain where explicit approaches, such as polarizable models, may be required.

As discussed previously, non-polarizable fixed-charge models are intrinsically unable to describe polarization effects and therefore cannot accurately predict the molecular dipole moment in condensed phases (e.g., liquids and solutions). In contrast, the dipole moment of classical polarizable models is able to respond on-the-fly to the surrounding environment, and thus, at least in principle, those models provide predictions of the liquid dipole moment. However, in most cases, they underpredict the magnitude of this dipole moment.⁵⁶ 

In principle, a much more accurate approach to calculate liquid-phase dipole moments is to use QM methods. However, this poses challenges of its own, as recently discussed.⁵⁶ The simplest approach, widely employed for many decades, is to apply a polarizable continuum method⁸⁰—i.e., the molecule of interest is surrounded by a dielectric continuum that represents the effect of the surrounding liquid. Despite its many successes, continuum methods describe solvent polarization effects in a mean-field way and hence do not account for local interactions such as hydrogen bonds. This assumption leads to a significant underestimation of the degree of polarization in polar solvents and, hence, to liquid-phase dipole moment values that are much lower than those of the real liquid.^56,60–63

At the other end of the complexity/accuracy scale are Ab Initio Molecular Dynamics (AIMD) methods, which share the advantages of classical MD (explicit description of individual molecules, and hence local effects) and QM calculations (explicit description of electronic degrees of freedom). Due to their very high computational expense, AIMD calculations are typically limited to relatively small systems and relatively low levels of QM theory—although such limitations are becoming less significant as computer power increases.⁸¹ More importantly, though, an explicit QM description of all molecules in the liquid introduces ambiguity in the separation of individual molecular contributions, and hence in the calculated molecular dipole moment values.⁸² The current state-of-the-art approach to decouple individual contributions to the dipole moment is to use maximally localized Wannier functions,^72,83 but this approach is not without its challenges.

A good compromise between accuracy and computational speed is provided by QM/MM methods,⁸⁴ where the central molecule is described at a QM level and solvent effects are represented by discrete classical particles. These methods allow for much larger systems and higher levels of QM theory than AIMD, because only one molecule is described at this level, yet they retain the explicit description of electronic effects on the central molecule. In contrast with dielectric continuum methods, they are able to describe both mean-field and local contributions to polarization, due to the explicit description of the surrounding solvent molecules. However, several conditions need to be observed for QM/MM methods to provide accurate values of liquid phase dipole moments,⁵⁶ including the need to self-consistently polarize the surrounding solvent molecules. The recently proposed Self-Consistent Electrostatic Embedding (SCEE) method strikes a good balance between accuracy and computational efficiency for dipole moment calculations, by replacing an expensive polarizable embedding calculation with three low-cost static electrostatic embedding calculations.^56,60

The SCEE method has, so far, been applied to calculate dipole moments of water, alkyl alcohols, and ketones. For water, the most recent value under ambient conditions, 2.76 D,⁶⁰ is in very good agreement with previous AIMD and QM/MM calculations (see Ref. 56 for a detailed comparison) and is well within the uncertainty range of the experimentally reported value of 2.9 ± 0.6 D.⁷¹ For alcohols, the dipole moment values are between 2.6 and 2.7 D,⁶⁰ again in good agreement with recent AIMD calculations of methanol.^85, Table II shows a summary of these previous results (shown as the induced dipole moment, i.e., the difference between the liquid and gas values) and a comparison with equivalent calculations using a dielectric continuum IEFPCM model.⁸⁶ Both water and alcohols are hydrogen-bonding fluids, and this has a strong effect on the induced dipole moment for those compounds, which in both cases is ∼0.9 to 1.0 D. This is much larger than what is predicted by dielectric continuum models (typically ∼0.3 to 0.4 D), due to the neglect of hydrogen bond effects in those models. Meanwhile, for ketones, which do not form hydrogen bonds in the pure liquid state, the SCEE induced dipole moment (∼0.8 D³²) is of the same order, but actually slightly lower than predictions from the IEFPCM method (∼1.1 D), presumably because local effects are negligible in those pure liquids. More extensive studies of molecules with a variety of functional groups are needed to confirm the generality of these trends.

Once the dipole moment for a pure liquid has been determined self-consistently, it can be used as a fully polarized solvent in the MM part of the calculation to determine the dipole moment of different solutes. The inherent assumption here is that the presence of the solute molecule does not significantly change the dipole moment of the solvent relative to that of its pure state. Using this approach, Jorge et al.⁶⁰ calculated the dipole moment of methanol in different solvents and compared it with the results of IEFPCM calculations (Fig. 2). The results show a substantial enhancement of the solute dipole moment when solvated in polar solvents (water and alcohols), significantly above the predictions of the continuum model. In contrast, when methanol is solvated in alkanes, the extent of polarization is small and is well predicted by a simple mean-field approach. Interestingly, solvation in ketones, which are hydrogen-bond acceptors but not donors, leads to a mild but non-negligible increase in the methanol dipole moment relative to the PCM calculations. This further supports the conclusion that the dipole enhancement is intimately related to the propensity for H-bond formation.

Recently, SCEE has been applied to compute the dipole moment of water over a wide range of thermodynamic conditions, spanning gas, liquid, and supercritical states.⁸⁷ It was shown that the induced dipole moment, being very significant at room temperature, decreases with temperature but is not very sensitive to pressure, at least within moderate pressure ranges. In the supercritical regime, a continuous transition was observed, with the dipole moment distributions gradually shifting from gas-like to liquid-like behavior. Interestingly, a linear correlation was observed between the SCEE dipole moment and the average number of hydrogen bonded neighbors of the central (i.e., QM) molecule (Fig. 3), which was valid over all thermodynamic states considered.⁸⁷ As expected, however, when a continuum dielectric model was used, which neglects local effects, such a correlation was not observed. Further evidence for the strong effect of hydrogen bonds on the dipole moment was provided by Bakó et al.,⁸⁸ who showed from AIMD simulations of room temperature water that the dipole moment is mainly determined by the number of hydrogen bonds formed by the central molecule, with only a small effect caused by different types of hydrogen-bonding environment.

One can see that if k assumes the form proposed in the original MDEC model [Eq. (9)], then the precise dielectric scaling rule proposed by Leontyev and Stuchebrukhov [Eq. (8)] is recovered. However, Eq. (18) is more general and allows for the predictions of any model (including polarizable models, although this has not yet been tested) to be reconciled with experimental data. Indeed, Jorge and Lue⁹⁰ applied this approach to the data shown in Fig. 4(a) and were able to practically eliminate the systematic deviations observed [see Fig. 4(b)]. The parameter k used by Jorge and Lue was determined empirically by fitting the dielectric constant of methanol [i.e., a single point in Fig. 4(b)], yielding a value of 1.26. Subsequent calculations of the dipole moment of liquid methanol using SCEE,⁶⁰ however, allow for k to be predicted from first principles by combining Eqs. (9) and (16). This yields k = 1.25, in very close agreement with the original value of Jorge and Lue, thus demonstrating the consistency of the PolCA approach.

Despite the above successes and the fact that the underlying theories of polarization have been known for quite some years, the mismatch between the PES and the DMS is still mostly “off the radar” of force field developers. For example, several recent fixed-charge water models, such as H₂O-DC,⁹² TIP4P-FB,⁹³ OPC,⁹⁴ and TIP4P/ε,⁹⁵ have explicitly included the experimental dielectric constant among target properties for parameter fitting. The authors of the TIP4P/ε model go even further and argue that fitting the dielectric constant is essential in order to achieve an accurate model⁹⁵ and have subsequently extended their approach to a few other compounds.^96–98 As explained above, a natural consequence of trying to fit both the PES and the DMS is that the resulting dipole moment of the water model has to strike a balance between being close to that of the real liquid (so as to describe the DMS) but not very far from the “half-polarized” level (which provides the optimal description of the PES⁵⁸). Indeed, the dipole moments of the above models are 2.48 D for OPC, ∼2.43 D for both TIP4P-FB and TIP4P/ε, and 2.42 D for H2O-DC, in comparison with ∼2.2 to 2.3 D for most widely used models, including TIP4P/2005,⁶⁹ reflecting this shift in balance between the DMS and the PES.

A corollary of the theories of polarization described in Sec. II is that one would expect to see a deterioration of the description of the PES for models fitted to the dielectric constant, due to the required increase in the dipole moment away from the optimal level. Unfortunately, due to the high degree of coupling between different model parameters (most notably point charges and Lennard-Jones parameters) and to the varying parameter fitting targets and protocols used in the past, it is not trivial to conclusively demonstrate this. Nevertheless, some advances in that direction have been recently made. Sedano et al.⁹⁹ attempted to answer this question by carrying out a systematic comparison (the so-called VA-test¹⁰⁰) between TIP4P/2005 and OPC in their ability to predict a wide range of experimental properties. TIP4P/2005, which is almost entirely consistent with the polarization theories of Sec. II, obtained a higher performance score (7.2 vs 6.3) than OPC, which is inconsistent with those theories, providing strong evidence in favor of the above hypothesis. Similar conclusions were reached by Perrone et al.¹⁰¹ in their systematic study of three-site water model optimization. Those authors concluded that adopting a scaling factor approach was able to improve agreement with dielectric constant predictions without significant deterioration in most other properties. Finally, the ECCw2024 model of Jungwirth and co-workers¹⁰² was developed according to the tenets of the MDEC/ECC theory of Leontyev and Stuchebrukhov^24–27 and hence was intentionally designed to reproduce a scaled-down version of the experimental dielectric constant of water, obtained by applying Eq. (8). The ECCw2024 model was shown to provide predictions of a wide range of water properties at least as accurately as state-of-the-art non-polarizable water models.

Despite these advances, however, a “smoking gun” that unequivocally resolves this issue is yet to appear, due to some limitations of the above studies: (i) the TIP4P/2005 model, despite its outstanding performance, is not 100% consistent with the above theories because it included the enthalpy of vaporization in the parameterization process after application of the (incorrect) Berendsen correction (see Sec. II B); (ii) the study of Perrone et al. is limited to three-site models, which are well known to underperform when compared with four-site models¹⁰¹ and therefore likely suffer from larger coupling between charges and LJ parameters; and (iii) the ECCw2024 model uses Eq. (8) instead of the combination of Eqs. (16) and (18), as advocated here as more appropriate for neutral molecules. Further work is therefore required for a complete resolution of this question.

As explained in Sec. II C, for non-polarizable models with dipole moments that are intermediate between the gas and liquid state dipoles, polarization corrections to phase change energies of the pure fluid (e.g., enthalpy of vaporization and self-solvation free energy) are zero or negligible. The same could be said of the vapor pressure, since it is related thermodynamically to the self-solvation free energy.¹⁰³ However, for solvation free energies where the solute and solvent are different compounds, energetic polarization corrections do not necessarily cancel out (see Table I). When the solute is non-polar, e.g., alkanes, polarization corrections are still expected to be nearly (or even exactly) zero, simply because the solute practically does not respond to changes in the surrounding environment. However, when the solute is polar and the solvent is less polar than the solute, polarization corrections are likely to be significant. Indeed, it has been shown that, when applied, such polarization corrections are able to improve the description of solvation free-energies of alcohols and ketones in alkane solvents.^31,32 However, these are just two examples, and more wide-ranging studies, involving molecules with functional groups of different polarities, are required to assess the validity of this hypothesis.

In principle, one could make use of the large datasets for simulated solvation free energies available in the literature (e.g., the FreeSolv database¹⁰⁴) to test whether application of polarization corrections improves agreement with experimental solvation free energies. Unfortunately, this is not as straightforward as it would seem, primarily because none of the previously tested force fields (e.g., GAFF and OPLS) has been parameterized in a way that is consistent with the polarization theories described in Sec. II. Whenever the point charges of a given model (whether empirically fitted or obtained from QM calculations) deviate from the optimal level of effective polarization determined by Eq. (16), the LJ parameters of that model, which are almost always empirically fitted, will compensate at least partially for any missing or excessive polarization energy. This parametric coupling implies that when those LJ parameters are combined with charges that better describe the effects of polarization, even if calculated rigorously from Eq. (16), there is no guarantee that the prediction accuracy of that model will improve.

Moreover, previous attempts to directly assess the “improved” versions of generic force fields that effectively account for polarization have used approaches that are not fully accurate. For example, some authors^59,64–66 have compared the performance of point charges obtained from the “standard” RESP method with those calculated using the IPolQ-Mod approach,⁵⁹ which in principle accounts for polarization effects through halfway-polarized charges. However, as discussed in Sec. II A, the IPolQ-Mod method, unlike the original IPolQ formulation,²⁹ is inadequate for this purpose because it uses a dielectric continuum model as the reference state for liquid-phase polarization, thus neglecting local contributions to the polarization process. It has been shown, as explained in Sec. III A, that such an approximation introduces rather serious errors in liquid dipole moment calculations. Other authors^105,106 have instead compared the performance of RESP charges against fully polarized charges obtained in a dielectric continuum model. Once again, this is inadequate because there is no guarantee that the charges are polarized to the correct extent. In fact, based on the SCEE liquid dipole moment calculations discussed in Sec. III A, it seems that polarizable continuum models yield charges that are approximately halfway polarized for some molecules (such as water and methanol), but fully polarized charges for others (e.g., acetone). This means that in all the studies mentioned above, there is no way to make sure that each solute molecule has the correct level of polarization, and this has led to inconclusive results.

In comparison with all the other properties discussed above, the impact of polarization on interfacial properties, such as the interfacial tension, is much less certain. Because of the non-isotropic nature of interfaces, one might expect that the dipole moment of interfacial molecules should be intermediate between that in the bulk liquid and bulk gas phases. Indeed, Kuo and co-workers carried out AIMD studies of the water/vapor¹⁰⁷ and methanol/vapor¹⁰⁸ interfaces and observed that the dipole moment decreased gradually by ∼0.6 to 0.7 D from the bulk to the interface. However, the uncertainty in their calculations is rather large (about ±0.3 D) due to the system size limitations of AIMD. The authors ascribed the dipole moment decrease to the weaker hydrogen bonds in the vicinity of the interface. Similarly monotonic decreases in the dipole moment between the bulk and the interface were reported from simulations with polarizable models.^109,110

One shortcoming of all those studies, however, is that they have sampled the interface using a “global” framework, which ignores the capillary wave-induced corrugations of fluid interfaces. A correct description of corrugated interfaces, and hence an unequivocal connection between the dipole moment of a particular molecule and its surrounding environment and exact location with respect to the interface, is only possible using intrinsic analysis methods.^111,112 To the author’s best knowledge, only a few studies examined interfacial systems described by polarizable models using intrinsic analysis. Kiss et al.¹¹³ applied the ITIM intrinsic analysis method¹¹⁴ to the vapor/liquid interface of water described by the polarizable BK model.¹¹⁵ They observed that the average water dipole moment in the surface layer was lower than in bulk, but all subsequent sub-surface layers exhibited a bulk behavior with regard to that property. Furthermore, the decrease in the dipole moment at the surface was much lower (∼0.1 D) than reported previously through non-intrinsic sampling methods.^107–110 

A more detailed study of the same polarizable water model was carried out by Hantal et al.¹¹⁶ as part of a series of papers on liquid/vapor interfaces of alkali halide aqueous solutions.^116–119 They calculated both intrinsic and non-intrinsic dipole moment profiles for the water/vapor interface, as well as detailed distributions of both total and induced dipole moments in subsequent layers beneath the surface. Their results confirmed that the dipole distribution (and hence also the average dipole moment) in the surface layer was ∼0.1 D lower than in subsequent layers and in the bulk liquid, so only the surface layer was affected by the change in polarization environment. This likely reflects the importance of local effects in the polarization of water, as discussed in detail in Sec. III A. They also observed an oscillatory behavior in the dipole moment orientation that only became apparent by applying intrinsic analysis. Furthermore, a direct comparison between their intrinsic and non-intrinsic profiles suggests that molecules that are located at the “fingers” of the corrugated interface, and therefore correspond to the outermost molecules in the global profile, have a much lower dipole moment than those located elsewhere on the surface layer. However, more in-depth studies are needed to confirm this interpretation.

The change in polarization environment, and consequently of the molecular dipole moment, at interfaces is likely to propagate to other interfacial properties, such as the surface tension. Attempts to assess if polarizable models provide better predictions of the surface tension than non-polarizable models have yielded inconclusive results (see Ploetz et al. for a detailed discussion¹²⁰). Rivera et al.¹²¹ simulated the water/vapor interface using a polarizable model and decomposed the contributions to the surface tension arising from different intermolecular interactions, concluding that polarization was responsible for about one third of the total positive contributions, hence could not be neglected. However, as discussed previously, it is not the case that non-polarizable models completely neglect polarization interactions, just that they include them implicitly into pairwise parameters. In fact, a later study by Iuchi et al.¹²² showed that a non-polarizable water model fitted specifically to match bulk properties of a polarizable model was able to reproduce the surface tension predicted by the latter. In an interesting recent study, Chen and Voth compared three models for the vapor/liquid interface of an ionic liquid.¹²³ They showed that significantly different properties were obtained in the polarizable and non-polarizable full-charge models, while a non-polarizable model with scaled ionic charges (see also Sec. III E) provided a reasonable, albeit not perfect, approximation of the behavior of the polarizable model. However, neither of the models was able to accurately predict the experimental surface tension for that system.

Thus, a more relevant question in the context of the present review is not if polarization contributions to the surface tension are important, but if polarization corrections should be applied to the surface tension predictions of non-polarizable models. Indeed, if one thinks of the surface tension as a surface free energy, then it stands to reason that polarization corrections analogous to those derived in Sec. II B (i.e., E_Dist and E_Elec) would apply at the interface, although to which extent they cancel out in the same way as observed for bulk liquids (see Table I) is as yet unknown.

The present paper focuses primarily on the effects of polarization on force fields for neutral molecules, and the polarization theories described in Sec. II have been discussed mostly from that point of view. However, polarization effects are arguably even more important in ionic systems, such as aqueous electrolyte solutions or ionic liquids, and therefore, it is worth ending this Perspective with some brief considerations about such systems.

For ionic liquids, it was observed that earlier “full-charge” models led to unphysically slow dynamics of the bulk liquid, significantly underpredicting self-diffusion coefficients and overpredicting viscosities.^124,125 However, when non-polarizable force fields were parameterized on the basis of a reduced net charge on individual ions (roughly between 70% and 90% of the formal charge), much better predictions of dynamic properties were obtained without compromising accuracy on structural or thermodynamic properties.^38,126–129 While many attempts to develop scaled-charge ionic liquid models used the net charge as an additional empirical parameter, some studies tried to rationalize this need for a reduced charge, often based on QM calculations on ionic liquid pairs and clusters, and invoking charge transfer, single ion polarization effects, or both.^37,127,130 However, such studies do not represent the true liquid environment experienced by the ions and neglect the important contribution of electronic screening of ion–ion interactions by the surrounding solvent, as described in the MDEC model.^23,24,131 Even for AIMD studies that more realistically represent the condensed phase environment (e.g., Ref. 132), it has not been possible so far to decouple all the different contributions that lead to the reduced net ionic charge for ionic liquid systems. In fact, the scaled-charge approach is far from being consensual, with several authors advocating the use of full-charge models with adjusted Lennard-Jones parameters and/or alternative functional forms.^133,134 Interestingly, however, a recent AIMD study on simple spherical ions demonstrated unequivocally that ion–ion interactions are indeed screened by ∼0.7 to 0.8 even when only electronic polarization is present,¹³⁵ lending strong support to the MDEC theory for these systems.

The uncertainty regarding the physical nature of the reduced charge for non-polarizable ionic liquid models may have prevented the idea of electronic screening from being adopted for parameterization of ions in other environments, such as electrolyte solutions. Nevertheless, the shortcomings of “full-charge” models for electrolyte solutions have been recognized, including excessive ion pairing in solution and too low solubility predictions.^136–138 As observed for ionic liquids, reducing the ion charges, again to something between 70% and 90% of their formal charge, significantly improved the performance of the models.^{34–36,75,76,139,140} However, while the benefits of scaled-charge over full-charge ion models are now well established, the actual degree of charge scaling that should be applied is still a matter of ongoing debate. While developing their MDEC model, Leontyev and Stuchebrukhov¹³¹ showed that interaction energies and potentials of mean force between ions in solution approached those obtained in QM calculations when the ion charges were scaled using Eq. (9). The ECC model of Jungwirth and co-workers adopts precisely this form of scaling^34,75,76 and is hence entirely compatible with MDEC theory. The scaling factor for the ion charge would thus be almost exactly 0.75 for aqueous solutions (since ε_∞ = 1.776 for water) and would be close to 0.7 for most ionic liquids. Meanwhile, Kann and Skinner suggested that the ion charge should be scaled by the ratio of the solvent dielectric constants from simulation and experiment.³⁵ This means that the scaling, as well as being solvent-dependent, would also be model-dependent. For example, for TIP4P/2005, the charges should be scaled by ∼0.85.¹⁴⁰ This idea has been followed by Vega and co-workers when developing their scaled-charge Madrid-2019 force field for electrolyte solutions.³⁶ More recently, however, Blazquez et al.¹⁴⁰ have shown that different scaling factors yield optimal predictions of different types of properties: 0.85 for thermodynamic properties (as in the Madrid-2019 model), 0.75 for transport properties, and 0.92 for interfacial properties. Consequently, they advocated considering the net ion charge as an additional fitting parameter during model development.

It is important to reiterate that the effectiveness of a particular charge scaling factor is likely to be intimately linked to the underlying water model³⁵ (in the case of aqueous electrolyte solutions), and hence, it is hard to decouple the two effects. In this regard, one might argue that the Madrid force field is compatible with TIP4P/2005,¹³⁹ while the ECC force field is compatible with the recently developed ECCw2024 water model,¹⁰² since the latter was specifically designed to yield a dielectric constant obtained from Eq. (8). However, neither TIP4P/2005 nor ECCw2024 are fully compatible with PolCA, as described in Sec. II C: TIP4P/2005 has a dipole that is slightly too high [2.305 D compared with 2.262 D predicted by Eq. (16)] and overestimates the enthalpy of vaporization in the absence of the Berendsen correction; conversely, the dipole moment of ECCw2024 is slightly too low (2.168 D). It would be interesting to develop a water model with a dipole moment of exactly 2.262 D and assess the performance of different ion models in combination with that model.

A final word is in order regarding the predictions of properties of ionic liquids or solutions that depend on the DMS. Just as demonstrated for molecular liquids, properties that depend on the DMS but are calculated using fixed-charge force fields parameterized to describe the PES require corrections before being compared with experimental data. The dielectric constant is obviously one such property (see Sec. III B), but so is the electrical conductivity. Indeed, it has been shown that good agreement with experimental conductivity data can be achieved by reanalyzing MD trajectories sampled with a scaled-charge model but reassigning the full charge to the ions.^141,142

Over the past 10–15 years, significant advances have been made in theoretical approaches to effectively account for polarization in non-polarizable force fields, and it appears as if the picture of how to optimally achieve this is finally emerging. However, a few pieces of the puzzle are still missing, and work to complete them is required. For example, proof that the thermodynamic cycles depicted in Fig. 1(b) (MDEC theory) and Fig. 1(c) (halfway charge theory) are formally equivalent would allow one to establish an unequivocal relation between the different energy terms considered in the two approaches—so far, as discussed in Sec. II C, this equivalence has been demonstrated only on the basis of empirical evidence. In addition, it would be desirable to develop a method to calculate E_Elec directly from QM (or QM/MM) calculations without requiring the presence of a dielectric continuum model, as the energy values calculated using the former approach can exhibit some dependence on the technical details of the continuum model (e.g., solute cavity radius and shape). This would allow for more accurate polarization corrections to be estimated for solvation free energies.

The PolCA workflow attempts to reconcile the two theories and proposes a practical way to parameterize non-polarizable force fields so that polarization effects are accounted for in a rigorous way. The first requirement is that point charges are assigned in an optimal way to approximate the correct level of polarization interactions in the liquid phase—i.e., yielding a dipole moment that is intermediate between that of the gas and the liquid states [more precisely, following Eq. (16)]. This leads to an almost exact cancellation of the positive and negative polarization energies that are not captured by fixed-charge models and hence obviates the need for polarization corrections to pure fluid phase-change energies (e.g., enthalpy of vaporization and self-solvation free energy). While point charges in most existing generic fixed-charge force fields are assigned values that generally lead to a higher dipole moment than in the gas phase, this assignment has so far been based on empirical or theoretically inconsistent approaches (although the original IPolQ method²⁹ is a possible exception). As such, it is hard to know a priori if a given molecule described using a particular force field possesses the right degree of polarization. In this regard, it would be interesting to carry out more detailed tests of generic force fields within the framework presented herein—e.g., calculating SCEE dipole moments for a wide range of compounds and checking how often Eq. (16) is obeyed.

On the basis of the current evidence, it is likely that to take polarization effects into account in a completely rigorous way would require a full reparameterization of generic fixed-charge force fields. Although recent efforts have been made to develop new models that are consistent with PolCA, this has only been done for a small number of compounds. In order to reach a much wider diversity of molecules and functional groups, it would be useful to embed the principles of PolCA within state-of-the-art force field development toolkits, such as those provided by the OpenFF initiative.¹⁷ Another option that would obviate the need for extensive parameter fitting exercises is to embed PolCA into the development of bespoke QM-derived force fields, for example, the QUBE approach of Cole and co-workers.^39,40,143 As discussed previously,⁵⁸ the main effort required would be to carry out SCEE calculations of liquid phase dipole moments for the entire QUBE calibration dataset (∼100 molecules), potentially followed by a re-optimization of a small number of required adjustable parameters—much less than a full reparameterization of all atom type parameters.

In terms of practical applications of the methods presented here, extending the data of Fig. 2 to a much wider variety of solutes and solvents would allow one to more extensively test the hypothesis of a direct relationship between the liquid phase dipole moment and the extent of hydrogen bond formation. In particular, consideration of solutes/solvents with only partial hydrogen bond formation capabilities (e.g., only donors or only acceptors) would be desirable. A more extensive database of SCEE liquid dipoles would also allow for a more comprehensive test of the validity of Eq. (18), not only applied to a variety of compounds but also over a range of thermodynamic conditions. For example, the recently reported data for the dipole moment of water over a wide range of temperatures and pressures⁸⁷ could be used to test dielectric constant predictions in comparison with experimental data.

As discussed at the end of Sec. III B, an unequivocal demonstration of the mismatch between the PES and the DMS for non-polarizable models would require a systematic and rigorous study that parameterizes two models using exactly the same approach but fitted to reproduce either the experimental dielectric constant (i.e., trying to fit both the PES and the DMS simultaneously) or the corrected dielectric constant obtained by applying Eq. (18) (i.e., fitting only the PES and then obtaining the DMS by charge rescaling). Recent studies provide strong evidence in favor of this hypothesis, although an unequivocal proof is still lacking. The same could be said about scaled-charge models of ionic liquids and electrolyte solutions—although empirically it is clear that they outperform full charge models, it would be good to demonstrate this using a systematic comparison whereby the two types of models are parameterized on the same basis.

Finally, the application of polarization corrections to interfacial properties as well as solutions and mixtures is largely underexplored. Much more work is needed along these lines, and this is likely to lead to innovative approaches to force field development and testing. In particular, it would be interesting to systematically test to which extent non-polarizable models can compete with fully polarizable models in describing these heterogeneous systems—in other words, how far can the principles of polarization energy cancellation and post facto polarization corrections be stretched to deal with non-isotropic systems. This would enable modelers to make a more informed choice regarding the type of model that is required for a particular system and hence choose the right balance between accuracy and computational efficiency.

The author is indebted to the multiple Ph.D. students and long-time collaborators who, over the years, contributed to the work discussed here, including (in no particular order) Cecilia Barrera, Andrew Milne, Zoe MacPherson, Chris Ringrose, Nuno Garrido, György Hantal, Daniel Cole, Javier Cardona, Leo Lue, and José Richard Gomes. Thanks are also due to several colleagues, who have made immense contributions to this field, for fruitful discussions, including (again, in no particular order) Bill Swope, Pál Jedlovszky, Marcelo Sega, Natália Cordeiro, Carlos Vega, Pavel Jungwirth, Bill Smith, Jiři Kolafa, and Carlos Nieto-Draghi.

The author has no conflicts to disclose.

Miguel Jorge: Conceptualization (lead); Investigation (lead); Methodology (lead); Writing—original draft (lead).

Data sharing is not applicable to this article, as no new data were created or analyzed in this study.