PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Theor Chem Acc. Author manuscript; available in PMC 2010 September 1.
Published in final edited form as:
Theor Chem Acc. 2009 September; 124(1-2): 11–28.
doi:  10.1007/s00214-009-0617-x
PMCID: PMC2888514
NIHMSID: NIHMS197600

Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability. Theory and applications

Abstract

A current emphasis in empirical force fields is on the development of potential functions that explicitly treat electronic polarizability. In the present article, the commonly used methodologies for modelling electronic polarization are presented along with an overview of selected application studies. Models presented include induced point-dipoles, classical Drude oscillators, and fluctuating charge methods. The theoretical background of each method is followed by an introduction to extended Langrangian integrators required for computationally tractable molecular dynamics simulations using polarizable force fields. The remainder of the review focuses on application studies using these methods. Emphasis is placed on water models, for which numerous examples exist, with a more thorough discussion presented on the recently published models associated with the Drude-based CHARMM and the AMOEBA force fields. The utility of polarizable models for the study of ion solvation is then presented followed by an overview of studies of small molecules (e.g. CCl4, alkanes, etc) and macromolecule (proteins, nucleic acids and lipid bilayers) application studies. The review is written with the goal of providing a general overview of the current status of the field and to facilitate future application and developments.

1. The need for polarizable force fields

Molecular mechanical (MM) or empirical force fields are widely used in molecular modeling and dynamics studies of biological and materials systems that contain 100,000 or more atoms. This capability is based on the utilization of simplified potential energy functions for determination of the energies and forces acting on large heterogeneous systems. However, such simplified forms of the energy function are also a weakness of MM methods as they limit their inherent accuracy. One of the major limitations with respect to the accuracy of the majority of current MM force fields is the way in which the charge distribution of the molecules is treated. Typically, effective partial fixed charges are assigned to the atoms independent of the configurations, which are adjusted to account for the influence of induced polarization in an average way. The functional form of the Coulomb interaction potentials thus created is not capable of adapting the charge distributions to changes of polarity in the environment. Such force fields, which are commonly termed “additive”, are currently used for most biomolecular simulations [1-5]. These force fields all share the same functional form to determine the potential energy as a function of the geometry, U (r) :

U(r)=Ubond(r)+UvdW(r)+Uelect(r)
(1a)

Uelect(r)=i=1Njiqiqjrij
(1b)

In equation 1 Ubond (r) represents the bonded terms (bonds, valence angles, dihedral or torsion angles, etc.), U vdW (r) represents the van der Waals (vdW) contribution, typically being a Lennard-Jones (LJ) 6-12 term, and U elect (r) is the electrostatic interaction term of the Coulomb form (Eq. 1b), where qi and qj are the partial atomic charges of atoms i and j separated by a distance, rij. While the functional form in equation 1 has been widely used (eg. the CHARMM [3], AMBER [2] and OPLS [6] force field papers in combination have been cited over 8500 times) the inability of the charge distribution to vary and adapt as a function of the local electric field is considered a major limitation of current models, significantly diminishing their ability to accurately treat intermolecular interactions in a variety of environment [7-10]. Accordingly, it has become clear that the inclusion of electronic polarization will play a central role in the next generation of force fields for molecular dynamics (MD) and Monte Carlo (MC) simulations [11,12].

In simple terms it is known that the molecular dipole moment of individual molecules change significantly when they are transferred from the gas to liquid phase. A prime example of the importance of electronic polarization to reproduce these phenomena is the dipole moment of water in different environments. In the gas phase, an isolated water molecule has a dipole moment of 1.85 D [13], but the average molecular dipole is 2.1 D in the water dimer and increases in larger water clusters [14]. In the condensed phase it reaches, a value between 2.4–2.6 D, as suggested from classical MD simulations of the dielectric properties [15-17], and 2.95 D, a value obtained from ab initio MD simulations [18-20] and from analysis of experimental data [21,22]. This bulk value is close to the maximum dipole any charged or polar molecule can induce in a water molecule. For example, the presence of a sodium ion or a dimethyl phosphate anion in bulk water does not significantly add to the induction effect that occurs in pure water [23].

By explicitly including polarization, a force field may be parameterized to reproduce accurate gas-phase quantum mechanical (QM) or experimental data and also perform well in condensed phases, due to the fact that it is able to respond to environmental effects. According to Rayleigh-Schrödinger quantum mechanical perturbation theory, a dipole linearly proportional to the local electric field from the environment is induced in the molecule. Thus, if the electric field is not too large, such that hyperpolarization effects are absent, the induced dipole μ on an atom is the product of the total electric field E and the atomic polarizability α.

μ=αE
(2)

The total electric field, E, is composed of the external electric field, E0, from the permanent atomic charges and the contribution from other induced dipoles. This is the basis of most polarizable force fields currently being developed for biomolecular simulations. Methods for this treatment of polarizability will be discussed in more detail in the following section.

It should be noted that the present review is biased towards the polarizable force fields implemented in the program CHARMM and towards applications to systems of biological interest. To all groups active in developing polarizable force fields and not referenced in this work, the authors wish to express their apologies. Much of the work not covered in the present review has been described in the special issue of the Journal of Chemical Theory and Computation dedicated to polarization [24].

2. Methodologies commonly used for polarizable biomolecular simulations

2.1 Induced dipole model

One method for treating polarizability consists of including both partial atomic charges and inducible dipoles on the atoms comprising the molecular system. In the most common variation currently in use [25-29], point inducible dipoles are added to some or all atomic sites while in the methodology proposed by Allinger and co-workers bond dipoles, in combination with atomic charges, are considered [30]. In the induced dipole approach the dipole moment, μ , induced on a site i is proportional to the electric field at that site, Ei. The proportionality constant is the polarizability tensor, α. The dipole feels an electric field both from the permanent charges of the system and from the other induced dipoles. The expression for μ is

μi=αiEi=αi[Ei0jiTijμj]
(3)

where E0 is the field from the permanent charges and Tij is the dipole field tensor.

A feature of the induced dipole polarizable model, as well as all polarizable models, is that the assignment of the electrostatic parameters is in principle easier than for additive models. Charges can be assigned based on experimental gas phase dipole moments or can be determined using QM ab initio methods and the polarizabilities can be obtained from the literature or QM calculations. For example, seminal work by Thole [31] and Applequist [32] showed how a set of simple polarizabilities based on atom type could reproduce the dipole moments and polarizabilities of a range of molecules. This is in contrast to nonpolarizable models, in which charges are systematically overestimated to have some enhanced permanent dipole moment to reflect the enhanced polarization required to accurately treat condensed phases [2,33,34]. Indeed, determining the degree of charge enhancement is part of the art of constructing additive potentials and constrains the utility of these potentials to a limited range of environments [12,35].

An implementation of the induced dipole method in CHARMM has been reported [36], based on the polarizable intermolecular potential functions (PIPF) model of Gao and co-workers [37,38]. The PIPF potential combined with the CHARMM22 force field has been designated PIPF-CHARMM, although the model has not yet been extended to cover all the amino acids. In this model infinite polarization is avoided by using Thole’s electrostatic damping scheme [31,39]. A method to accelerate the convergence of the induced dipoles for systems employing the PIFF potential functions has been described [40].

An approximation to the induced dipole model was proposed by Ferenczy and Reynolds [41]. This induced charge method involves point charges only, and those depend on the environment. It is based on the idea of representing atomic point dipoles by point charges on neighboring atoms. The method was extended so that both the polarization energy and its derivatives can be determined.

More recently, Ponder and co-workers [42-49] developed the AMOEBA force field based on a modification of the formulation of Applequist and Thole. It uses a modification of Eq. 3 with the static electric field typically treated with permanent atomic charges replaced by permanent multipoles:

μi=αiEi=αi[jiTijαMj+kiTikαβμk]
(4)

where M (Mi = (qi,μi,x,μi,y,μi,z,Qi,xx,Qi,xy,Q,xzQi,zz)T) is the vector of permanent atomic multipole components, up to quadrupole, and T is the interaction matrix. Accordingly, the induced atomic dipoles must respond to the contribution of the multipoles as well as the contribution of the other induced dipoles to the electric field. While such an extension represents an increased computational demand, the inclusion of atomic multipoles more accurately treats the interactions between molecules as a function of orientation, avoiding the inclusion of particle representative of lone pairs as has been incorporated into the CHARMM Drude force field (see below).

2.2 Classical Drude oscillator model

Another method to include polarization consists in modeling the polarizable atomic centers using dipoles of finite length, represented by a pair of point charges. A variety of different models of polarizability have used this approach, but especially noteworthy is the classical Drude oscillator models (also known as the “shell” or “charge on spring” model) frequently used in simulations of solid state ionic materials and recently extended to water and organic compounds. The Drude model can trace its origin to the work of Paul Drude in 1902 and was developed as a simple way to describe the dispersive properties of materials [50]. It represents electronic polarization by introducing a massless charged particle, attached to the atomic center of each polarizable atom by a harmonic spring. The position of these “auxiliary” particles is then adjusted self-consistently to their local energy minima for any given configuration of the atoms in the system, thereby taking into account the permanent electric field due to the fixed charges and the contribution of the induced dipoles to the electric field. A quantum version of the model (including the zero-point vibrations of the oscillator) has been used in early applications to describe the dipole–dipole dispersion interactions [51-55]. A semiclassical version of the model was used more recently to describe molecular interactions [56], and electron binding [57]. The classical version of the model has been quite useful in statistical mechanical studies of condensed systems and in recent decades has seen widespread use in MD or MC simulations. Example applications include ionic crystals [58-63], a range of simple liquids [36,64-70], liquid water [71-77], and the hydration of small ions [78,79]. In recent years, the Drude model was extended to interface with QM approaches in QM/MM methods [80]. A particularly attractive aspect of the Drude oscillator model is that it preserves the simple particle-particle Coulomb electrostatic interaction, such that its implementation in standard biomolecular simulation programs may be performed in a relatively straightforward way.

In the Drude oscillator model polarization is determined by a pair of point charges separated by a variable distance d. For a given atom with charge q assigned to the atomic center a mobile Drude particle (or Drude oscillator) carrying a charge qD is introduced. The charge on the atom is replaced by qqD in order to preserve the net charge of the atom-Drude oscillator pair. The Drude particle is harmonically bound to the atomic particle with a force constant kD. The mathematical formulation of the Drude model is, in fact, an empirical method of representing the dipolar polarization of the atomic center on which it is introduced. A related method to introduce polarization in FF simulations was developed by Sprik and Klein {Sprik, 1988 #97}, where polarization is represented by closely spaced point charges. In their water model, four charges are placed around the oxygen in addition to the three permanent charges on the oxygen and hydrogen atoms. The introduction of many rigid point charges makes this approach computationally more expensive though it may possibly be more stable since there are no moving particles.

In the absence of an electric field during an MD simulation, the Drude particle oscillates around the position of the atom, r, and the atom appears on average as a point charge of magnitude q. In the presence of a uniform field E, the Drude particle oscillates around a displaced position r + d. The Drude separation d is related to kD, E and qD :

d=qDEkD
(5)

and the formula for the induced atomic dipole, μ, as a function of E is

μ=qD2EkD
(6)

which results in a simple expression for the isotropic atomic polarizability, α,:

α=qD2kD
(7)

Therefore, in the Drude polarizable model the only relevant parameter is the combination qD2kD that is responsible for the atomic polarizability. It is of utility to reiterate that the electrostatic interaction in the Drude model is implemented using only the Coulombic term (Eq 1b) already present in MM simulation codes. No new interaction types, such as the dipole field tensor Tij of Eq. 3 are required. The great practical advantage of not having to compute the dipole-dipole interactions is balanced by the extra charge-charge calculations. However, significant computational savings may be gained in the Drude model by only attaching Drude particles to the non-hydrogen atoms that dominate the molecular polarizability, thereby increasing the total number of interactions pairs by a factor much smaller than 2 [70,76].

The polarizable Drude model in CHARMM under development in our laboratories is geared towards proteins, lipids, and nucleic acids [68-70,76-79,81-84]. Important progress has been made thus far. The algorithm has been generically described in reference [81] including the formulation for MD simulations based on an extended Lagrangian as required for computational efficiency. A description of extended Lagrangian integrators for MD simulations is done below. Two water models, that are a generalization of the TIP4P model, have been developed and are referred to as the SWM4-DP (Simple Water Model, 4 points with Drude polarizability) [76] and SWM4-NDP (Simple Water Model, 4 points with negative Drude polarizability) [77]; the later model will act as the basis of the full biomolecular force field. Protocols to determine the partial charges, atomic polarizabilities [82] and atom-based Thole damping factors [83] have been presented and parametrization of a number model compounds representative of biomolecules have been published, including alcohols [69], alkanes [85], aromatic [70] and heteroaromatic compounds [10], ethers [86] and amides [84]. Studies of ions in aqueous solution were also performed [78,79]. Notable extensions of the model include the inclusion of lone pairs and anistropic atomic polarizabilities on N, O and S hydrogen bond acceptors [83].

The parametrization protocol developed for the polarizable Drude model in CHARMM is well defined. A procedure for determining atomic center and Drude charges [82,83] and atom-based Thole damping factors [84] has been developed and is analogous to work by Friesner and co-workers [25,35,87-89], A series of maps of the electrostatic potential (ESP) that surrounds the model compound monomer are evaluated using QM density functional theory on a set of specified grid points, differing by the presence and/or location of a perturbing ion in the environment surrounding the molecule. Electrostatic parameters are then fitted to minimize the difference between the QM and Drude ESP maps [82,83]. Optimization of the components of the parameters not dependent of the Drude oscillator positions, namely the bonding and Lennard-Jones terms, are adjusted as described previously for the additive CHARMM force field [34 ,90].

2.3 Fluctuating charges model

Polarizability can also be introduced into standard energy functions (Eq. 1) by allowing the values of the partial charges to respond to the electric field of their environment, thereby altering the molecular polarizability. This may be performed by coupling the charges to their environment using electronegativity equalization (EE) or chemical potential equalization (CPE) schemes. This method for treating polarizability has been called the “fluctuating charge” method [27,91], the “chemical potential (electronegativity) equalization” method [92-107], or the “charge equilibration” method [108-114] and has been applied to a variety of systems. Examples include application to liquid water [91], vapor-liquid equilibrium [115], studies of ions in aqueous solution [116,117], studies of peptides [88], aqueous solvation of amides [118] and water and cation-water clusters [119]. A practical advantage of this approach is that it introduces polarizability without introducing new interactions. Compared to the Drude model, this can be done using the same number charge-charge interactions as would be present in a nonpolarizable simulation. However, while the fluctuating charge model does not introduce any additional terms or particles as compared to additive force fields it does require a significantly shorter integration time step for stable MD simulations [8], leading to a significant increase in computational costs over additive models.

In the fluctuating charge (FC) method [108], as it is commonly referred to in liquid state and biomolecular force field studies [27,91], variable discrete charges are located on atomic sites within the molecule. Their value is computed, for a given molecular geometry, by minimization of the electrostatic energy. In a multi-molecular system with Nmolec molecules and each molecule consisting of Natom atoms and Nsite charged sites, the electrostatic energy, Uelect(r,q) is as follows:

Uelect(r,q)=i=1Nmolecα=1Nsite[χiα0qiα+12Jiαiα0qiα2+i=1Nmolecα=1Nsiteβ>αNsiteJiαiβ(riαiβ)qiαqiβ+i=1Nmolecj>1Nmolecα=1Nsiteβ>αNsiteJiαjβ(riαjβ)qiαqjβ]
(8)

The energy given by Eq. 8 replaces the Coulomb energy qiqj/rij in Eq. 1. In Eq. 8 χα0 is the “Mulliken electronegativity” and Jαα0 is the “absolute hardness;” [120] these terms represent the electrostatic parameters in the fluctuating charge model and are optimized to reproduce molecular dipoles, interactions with water and the molecular polarization response, typically determined from QM calculations [27]. The charges qi are thus treated as independent variables, and the polarization response is determined by variations in the charges. These charges depend on the interactions with other molecules as well as other charge sites on the same molecule, and will change for every time step or configuration sampled during a simulation.

In most cases, charge is taken to be conserved for each molecule, so there is no charge transfer between molecules. However, in quantum mechanics charge transfer is an important part of the interaction energy, so there are reasons to remove this constraint [121-125]. Unfortunately, this procedure often leads to large overestimation of the polarizability as the molecular size increases as charge can now flow along covalent bonds at a small energetic cost. Thus, this method is suitable for small molecules but generally not applicable to macromolecules.

A solution to the over-polarization problem was developed in terms based on the concept of atom–atom charge transfer (AACT) [126]. In this approach the energy is Taylor expanded in terms of charges transferred between atomic pairs within the molecule, rather than in terms of the atomic charges themselves. Similar in spirit is the bond-charge increment (BCI) model [87,88], which allows for charge to only flow between two atoms that are directly bonded to each other, the method guarantees that the total charge of each set of bonded atoms is conserved. A related approach is the Atom-Bond Electronegativity Equalization Method (ABEEM) [127-131] which has been developed based on concepts from density functional theory. In this model, the total electronic energy of a molecule in the ground state is a complex function of different quantities: (a) valence-state chemical potential of atom a, bond ab and lone-pair electrons, lp, (b) valence-state hardness of atom a, bond ab, and lp, (c) partial charges of atom a, bond ab, and lp, and (d) distances between the different atoms, bonds and lone pairs. ABEEM has been successfully incorporated into the intermolecular electrostatic interaction term in MM models of water [132,133].

The interpretation of the FC model in the framework of QM theory is well defined. It is possible to derive the FC terms from density functional theory, from which the concept of electronegativity equalization arises naturally [134] or, as developed by Field, in terms of semiempirical MO theory [135]. A major difference between the FC and semiempirical MO methods is the arbitrariness in defining the atomic charges. In the FC model the charge on each atom can take any value with the restriction that the sum of atomic charges equals the total molecular charge. In the semiempirical MO methods the atomic charges are limited by the occupation and form of the orbitals. Despite this difference, the connection of the FC and MO methods can be explored to reduce the arbitrariness in the introduction of polarization inherent to the induced dipole or Drude models.

There are other approaches to include polarization in MD simulations. For example, a pure QM based method, conceptually related to the fluctuating charge model, was developed by Gao [136,137]. Each molecule is treated with a QM method, for example AM1 was used in ref [137], and the remaining molecules are represented by a Hartree product of the individual monomer wavefunctions. In this approximation exchange interactions are neglected and a LJ term is included to compensate.

2.4 Implementation of extended Langrangian integrators for MD simulations

An essential feature of MM methods for the treatment of biomolecular systems is their computational efficiency. The inclusion of polarizability increases the computational demand due to the addition of dipoles or additional charges centers and, in the context of MD simulations, the requirement for shorter integration timesteps. In addition, for every energy or force evaluation it is necessary to solve for all the polarizable degrees of freedom in a self-consistent manner. Traditionally, this is performed via a self-consistent field (SCF) calculation in which the induced polarization is solved iteratively until a satisfactory level of convergence is achieved [74,138,139]; the SCF equation can also be solved in a single step with matrix inversion [140-142]. While these methods have been widely used they involve significant computational costs (converging the SCF procedure requires about 15 iterations and the inversion of a large matrix is slow), making their use in MD simulations problematic. To overcome this MD integrators have been developed in which the polarizable degrees of freedom are included as dynamic variables. These approaches, whose origins go back to the Carr-Parrinello approach for quantum mechanical simulations [143], are referred to an extended Langrangian methods [144-146]. Such approaches have been developed for the majority of polarization methods discussed above. These include implementations for induced dipoles [147], Drude oscillators [81] and fluctuating charge methods [91,148]. In the case of the Drude oscillator model, it was explicitly demonstrated that the dynamics of an extended Lagrangian system, in which a small mass is attributed to the auxiliary Drude particles and the amplitude of their oscillations away from the local energy minimum is controlled separately with a low-temperature thermostat, provides a close approximation to the SCF regime [81]. The availability of extended Lagrangian methods are central to the future success of MD simulations that include electronic polarization, allowing the methods to approach computational speeds approaching that of the more approximate additive force fields.

3. Application of polarizable force fields

3.1 Water simulations

Water is an essential component in the chemistry of life, and a high quality water force field is essential for meaningful simulation studies of biological systems. Furthermore, any effort to develop a force field for biomolecular systems must start with a model for water. To meet the need of a computationally tractable yet appropriately accurate model of water an extremely large number of polarizable potentials for water have been developed. These include models based on induced point dipoles [140,149-167], Drude oscillators [72,74-77], fluctuating charge models [15,25,91,101,115,148,168,169] and hybrid induced dipole/fluctuating charge methods [35]. The reader is referred to a previous review on a number of additive and polarizable water models for additional information [16].

Polarizable water models generally perform a good job in reproducing both water dimer interactions energies and enthalpies of vaporization of liquid water, a combination that may be considered a minimum requirement for a model that will be appropriate for a range of environments. Most polarizable water models have dimer interaction energies close to −5.0 kcal/mol, the accepted ab initio value [170]. Examples include −4.69 kcal/mol for the POL3 model [164], −4.51 kcal/mol for the TIP4P-FQ model [91], −5.33 kcal/mol for the model of Burnham et al. [171], −5.00 kcal/mol for the MCDHO model [73], −5.0 kcal/mol for the POL5/TZ model [35] and −5.2 kcal/mol for the SWM4-DP and SWM4-NDP models [76,77]. For comparison the value of −6.5 kcal/mol for the additive TIP3P is significantly overestimated, as required to obtain accurate pure solvent properties [1].

The ability of polarizable models to accurately treat both the gas phase water dimer energy and bulk water properties allows them to perform better in reproducing the molecular dipoles in environments where the hydrogen bonds network is perturbed [14,171]. Explicit account of polarization seems essential to accommodate the local disruption of the hydrogen bond network created by anions such as chloride [116,150,172,173] or fluoride [174-178] or to reproduce the polarization effects of small multivalent cations on the first hydration shell [179,180]. Although explicit polarizability does not appear to have any significant effect on the reorganization of water molecules at liquid–hydrophobic [181] or liquid–vapor [164,182] interfaces, it may play a decisive role for the specific water–water interactions near small nonpolar moieties [183,184]. In short, polarizability is essential to obtain accurate energetics in the vicinity of highly polar moieties (such as carbonyl groups), small ions (such as sodium or chloride), and also in anisotropic nonpolar environments.

Liquid phase properties of several polarizable force fields have been studied in great detail, including structural properties (eg. radial distribution functions, RDF), dielectric constants, and dynamical properties (eg. diffusion constant and NMR relaxation times). The structure of liquid water is characterized by a short range order and a long-range disorder. This is reflected by the radial distribution function g(r), which can be derived from neutron [185-191] and X-ray [192-196] scattering experiments. Sorenson et al. [197] and Head-Gordon and Hura [198] provide a summary of experimental and simulated atom-atom RDF results. Over the years Soper et al. [185,190,199,200] have reported different results emphasizing the difficulty in unambiguously determining RDFs from experimental scattering data. However, two groups, including results from Soper, have now reported almost identical experimental RDFs based on independent analysis of neutron scattering experiments [190] or X-ray scattering experiments [196], representing the best RDF estimates currently available. It should be noted that the latest structural data reported by Soper in 2000 [190] is a revised analysis of the experimental data obtained in 1986 [185].

Head-Gordon and Hura [198] compared partial correlation functions of empirical water models such as the non-polarizable TIP5P model [201], the rigid polarizable model based on fluctuating charges, TIP4P-pol-1 [115], and a flexible polarizable model based on induced dipoles NCC-vib [202]. The TIP5P five-site additive model is particularly noteworthy given its excellent prediction of the gOO(r) data. The polarizable models also show very good agreement with gOO(r) and, without fitting, often predict the temperature of maximum density of water. None of the models does an outstanding job in reproducing the reanalyzed neutron scattering gOH(r) and gHH(r) data. Sorensen et al. [197] compared experimental radial distribution functions with predictions made for both polarizable and non-polarizable water models. They found that the calculated RDFs for the polarizable models (CC [163], TIP4P-FQ [91]) were generally in better agreement with experiment than those for the additive TIP3P [1] and SPC [203] three-site models. The RDFs for the polarizable force fields, however, were not significantly better than those for the four-site TIP4P model and not as good as those for the additive five-site TIP5P model [201].

Additional polarizable water models were analyzed by Sorenson et al. [197]. They include the fluctuating charge version of TIP4P, TIP4P-FQ [25], an extension of TIP4P that introduces an additional coupling between the Lennard-Jones interaction parameters for the oxygen site and their partial charges [115], TIP4P-Pol-1, an extension of the MCY water model to include flexible bonds and angles, as well as many-body effects, NCC-vib [202], the polarizable point charge model [139], PCC, and a simple polarizable model developed by Chialvo and Cummings to reproduce water properties over a wide range of conditions, CC [163]. The CC model shifts all RDF peaks to larger r, and has a very large peak and shows a loss of density at the first minimum in the RDF. While its structure does least well among the polarizable water models, its reproduces non-ambient states better [163]. The NCC-vib model [202] also overemphasizes the loss of density at the first minimum, but is otherwise good in reproducing the experimentally determined gOO(r). For the TIP4P-FQ model and the PPC models [139], it is evident that the overall agreement is excellent, although the position of the experimental first peak is not as well-reproduced as the nonpolarizable TIP5P model. The TIP4P-pol-1 water model shows improvements in first peak positions relative to TIP4P-FQ, but overemphasizes the loss and gain of density at the first minimum and second maximum, respectively. While the polarizable models perform well overall, many of these models are not as optimal performers at ambient temperature as their nonpolarizable partners. This and their inability to reproduce the RDF as well as the TIP5P model suggests that improved polarizable water models are accessible. Such models may need to include some representation of anisotropic charge distributions (e.g. off-center charges or ‘lone pairs’), as well as polarizability.

The ability of models to describe various properties of water in a broad range of thermodynamic states has been accounted to different degrees by the different force fields. The inclusion of the polarizability improves the performance of the water models in various respects. This includes reproduction of the temperature of maximum density [139,204,205], and describing the elongation of the hydrogen bonds with increasing temperature [206]. Interestingly polarizable water models proved to be less successful than some of the additive models in reproducing, for instance, the thermodynamic properties of water around the critical point (Table 1) [206]. It has been suggested that better polarizable models may require a different form for the vdW repulsive energy or a more realistic description of polarizability than is possible with a single polarizable atomic center [206]. Another possibility is that molecular flexibility may need to be included in these models. This conjecture is based on the fact that studies have shown that introduction of flexible bond angles (and bond lengths) affects the thermodynamic properties of nonpolarizable models considerably [132,133,207-209]. Perhaps, in taking one step forward (inclusion of polarizability) it is necessary to take another (inclusion of flexibility) to yield a more balanced representation of water.

(Table 1)
Estimates of critical temperature and density of polarizable and nonpolarizable water models.

More recent polarizable water models developed for biomolecular simulations are the SWM4-DP and SWM4-NDP model of Lamoureux et al. [76,77] and the AMOEBA model of Ren and Ponder [45,46]. The SWM4-DP model is a 4-point rigid model, analogous to the TIP4P model, with polarizability described by a Drude oscillator on the oxygen; the SWM4-NDP model is similar to the original SWM4-DP but carries a negative charge on the Drude particle on the oxygen to mimic the electron charge and has different LJ parameters. The RDFs for these models are characterized by the narrow shape of the first peak in the gOO(r) radial distribution and the first intramolecular peak of the gOH(r) distribution is slightly outward (rOH(1)=1.85 instead of 1.78 Å). The height of the first peak in the oxygen-oxygen radial distribution function is 3.07, which is somewhat high but almost within the experimental error from the neutron diffraction results (gOO(1)(r)=2.7±0.3 [216]). Notably, the SWM4 models accurately reproduce the diffusion coefficient and are in satisfactory agreement with experiment for the Debye and NMR relaxation times, indicating their accurate treatment of dynamic properties. These models also accurately reproduce the dielectric constant of water, ε, which is not surprising given that ε was included as target data during optimization of those models. Interestingly, that effort lead to the observation that the gas phase polarizability of water, as well as other molecules, may not be appropriate for the condensed phase. This is consistent with observations based on quantum calculations, as discussed in detail below. While the question of polarizability scaling is still being addressed (see below), it should be emphasized that a proper treatment of the dielectric behavior of water as well as other molecules is important for accurate treatment of solvation energies in different environments and, accordingly, its accurate reproduction by a model may be considered an essential feature.

The AMOEBA water model of Ren and Ponder [45] is fully flexible and was compared with experimental and QM data. Studies of single isolated molecules, molecular clusters, liquid water and ice were performed. AMOEBA calculated dipole and quadrupole moments and polarizability of an isolated water molecule were shown to be in good agreement with experimental and QM results. Tests of the water dimer were also conducted and it was found to be in good agreement with recent theoretical results [217,218]. Bonding energies and geometries of small water clusters from the trimer to the the hexamer were also found to be in good agreement with QM results [219,220]. Simulations of liquid water were performed and thermodynamic, transport and structural results were compared with experimental data. Of the several quantities computed, density at room temperature and heat of vaporization were in excellent agreement with experimental values, the dielectric constant was slightly higher than the experimental value, the self-diffusion coefficient was lower than the experiment and the viscosity was higher as expected. The structure of liquid water was characterized by computation of O…O, O…H and H…H RDFs sampled from NPT simulations and compared with Soper’s 2000 results [190]. It is noteworthy that the experimental curves gOO reported by Soper in 2000 (neutron-diffraction) are almost identical to the ones of Sorensen et al. (X-ray experiments) [197]. The position of the first peak in the gOO (r) radial distribution from the AMOEBA simulations is 0.08 Å longer and its height is higher than that of the Soper 2000 RDF. The first peaks of the gOH (r) and gHH (r) RDFs are also higher that the Soper 2000 data and the positions of the two peaks are shifted to larger distances. In general, the model provides a credible description of the structural properties of bulk liquid water at room temperature. Studies of two ice forms, ice Ih and XI, were also reported through energy minimization of atomic positions and crystal lattices and MD simulations. The computed results are in very good agreement with the experimental data. The same authors published another study where the temperature and pressure dependence of the AMOEBA water model are analyzed [46].

3.2 Ion Solvation

One of the most critical needs for a biological polarizable force field is the treatment of both atomic and molecular ions. Ion solvation is important in chemistry, including surface chemistry, environmental chemistry, and the study of molecules such as surfactants, colloids, and polyelectrolytes. Biologically, ions are critical to the structure and function of nucleic acids, proteins, and lipid membranes and ion transport in and out of the cell plays a central role in numerous physiological processes [221-225]. The structure of nucleic acids is affected by nonspecific counterion condensation [226] as well as specific interactions [227,228]. Ion binding to specific protein sites occurs for purposes of stabilization as well as playing central roles in enzyme catalysis [229-231]. Ion permeation across the cell membrane is tightly controlled by specialized proteins called ion channels [232-234]. In physical chemistry, ion solvation is also important in several processes such as chemical purification [235] and chromatographic systems [236] and ion-specific chelators [237].

Simulations of aqueous ionic solutions using nonpolarizable additive force fields have shown that consideration of non-additive effects is important to accurately reproduce the atomic details of ion hydration [43,116,150,173,238-240]. In principle, accurate potential functions for computer simulations can be developed and validated by comparing to experimental data (gas phase and bulk) and to the results of high level QM ab initio computations performed on ion-water dimers and small ion-solvent clusters [43,79]. Experimental target data available to parametrize and validate MM models of ion-water systems include gas phase energies of small hydrated clusters [241], bulk hydration free energies [242,243], structural properties (radial distribution function, coordination numbers, etc…), and transport coefficients (diffusion constant, mobility, conductivity). Åqvist was the first to develop additive ion-water interaction potentials for the most common ions in biology using calculations of the absolute hydration free energy in bulk water [244]. The models were constructed for the nonpolarizable SPC water model, but have also been translated for the TIP3P and TIP4P model, giving rise to a number of unanticipated issues (see recent article by Cheatham [245]). A similar route was taken to develop an independent set of ions for the TIP3P model [246,247].

Correct interpretation of the experimental hydration free energies of ions data and its use in constructing an accurate computational model is, in fact, not as straightforward as one would wish [78,79]. One difficulty arises because experimental thermodynamic or electrochemical measurements involve neutral macroscopic systems and have access only to irreducible, conventional hydration free energies that are either the sum of the absolute free energies of an ion and a counterion, ΔGhyd(M+)+ ΔGhyd(X), or the difference of the free energies of two ionic species of the same valence, ΔGhyd(M1+)+ΔGhyd(M2+). In fact, the absolute hydration free energy of a single ion cannot be resolved from calorimetric or electrochemical experiments alone. So, while the solvation free energy for a neutral salt can be measured, it is impossible to separate it experimentally into contributions from the cation and anion [242,248,249]. An additional extrathermodynamic assumption is required to perform this dissection [250]. A detailed explanation of the pros and cons of the different methods can be found in reference [43] and in [78]. In the same references alternative methods were proposed to deal with this assumption. Lamoureux and Roux examined the sensitivity of the bulk hydration free energy of individual ions to the gas phase monohydrate energy and showed that the absolute scale varies only within a narrow range. The implication is that the gas phase monohydrate energies puts, by itself, a tight constraint on the absolute scale of hydration free energies of ions. This analysis was used to set the absolute scale of hydration free energy and develop a consistent parametrization for the complete alkali-halide series. This type of analysis, relating monohydrate and bulk solvation properties, goes back to the pioneering study of Åqvist [244].

Several studies of the solvation of ions and salts have been published using the different polarizable models implemented in CHARMM. Lamoureux and Roux developed polarizable potential functions for the hydration of alkali and halide ions using the SWM4-DP water model. Patel and co-workers published several studies of solvation of ions and salts in water using the polarizable TIP4PFQ water model implemented in CHARMM [251-253]. In another study from that laboratory [251] the TIP4P-FQ model was compared with the additive TIP4P and the Drude water models. Roux and co-workers have presented a study of aqueous solvation of K+ and compared ab initio, polarizable (SWM4-NDP Drude water model of CHARMM) and additive force field methods [79]. All computational methods yielded hydration numbers between 5.9 (Car-Parrinello PW91/pw) and 6.8 (Drude model) in good agreement with experimental data (6-7).

Other authors have presented studies relevant to understand the performance of electronic polarization in interactions of molecules with ions. Masia and co-workers [254] studied the interaction of a molecule with a cation, via induced dipoles and Drude oscillators. The dimer electric dipole moments as a function of the ion-molecule distance for selected molecular orientations was compared with high-level ab initio calculations for water or carbon tetrachloride close to Li+, Na+, Mg2+, and Ca2+. It was shown that the simple polarization methods are able to satisfactorily reproduce the induced dipole moment of the cation-molecule dimer. In a previous study, the same authors studied the interactions between molecules and point charges [255].

A potential model for Li+-water clusters was presented [256] and the same authors performed a detailed study of the monovalent ions, Li+, Na+, K+, F, Cl, and Br in aqueous solution and the small water clusters M+(H2O)n and M(H2O)n. These studies were based on the ABEEM/MM method. Analysis of results from these studies included solvation structures, charge distributions, binding energies, dynamic properties (diffusion coefficients of ions) and free energies of hydration [257]. The computed quantities were found to be in good agreement with experimental results.

A study of solvation dynamics of divalent cations in water was performed by Piquemal et al. [258] using a modified AMOEBA force field. The model consisted of a cation specific parametrization based on ab initio polarization energies computed by a constrained space orbital variation (CSOV) energy decomposition method [259]. Excellent agreement between computed and experimental condensed phase properties was found despite the use of parameters derived from gas phase ab initio calculations.

A number of studies of the influence of ions on the air-water interface using polarizable models have been published. A useful review of these studies is that by Jungwirth and Tobias [260]. Examples include the study of Salvador et al. of the aqueous solvation of NO3 in interfacial environments with a Car-Parrinello MD simulation of a cluster and classical MD of an extended slab system with bulk interfaces using a polarizable force field based on the Atoms in Molecules analysis [261]. Both in aqueous clusters and in systems with extended interfaces the nitrate anion clearly prefers interfacial over bulk solvation. Archontis and co-workers studied the distribution of iodine at the air-water interface using the Drude based SWM4-DP water model [262,263] and Jungwirth and co-workers performed similar studies using the AMOEBA force field [264]. In all cases it was found that iodine tends to remain closer to the water/vapor and that tendency to stay at the interface increases in the order Cl < Br < I. More recently, studies of several ions and salts at the air-water interface have been published. Tobias and co-workers have presented a mixed X-Ray photoemission spectroscopy/MD study using polarizable potentials of aqueous potassium fluoride solutions [265]. Wang et al. studied NaCl using a Drude model [266] and Warren and Patel compared several polarizable ion models [252,253]. MD studies of aqueous solutions of molecular ions have also been published. Picalek et al. [267] studied the interfacial structure of aqueous solutions of 1-butyl-3-methylimidazolium tetrafluoroborate and 1-butyl-3-methylimidazolium hexafluorophosphate using both non-polarizable and polarizable force fields.

Specialized uses of combined QM and MD methods have been presented. Roos and co-workers [268] studied the coordination environment of the uranyl ion in water. Pair potentials were initially calculated using multiconfigurational wave function calculations and the quantum chemically determined energies were used to fit parameters in a polarizable force field with an additional charge transfer term. Classical MD simulations were performed for the uranyl ion and up to 400 water molecules. The results show a uranyl ion with five water molecules coordinated in the equatorial plane. The U-O(H2O) distance is 2.40 Å and a second coordination shell starts at about 4.7 Å from the uranium atom. Interestingly, no hydrogen bonding is found between the uranyl oxygens and water.

In summary, polarizable force fields have been shown to be very promising for simulating the properties of ionic systems. Particularly satisfying, though not surprising, is their ability to successfully describe interfacial systems more accurately than additive models. This further emphasizes the importance of the ability of polarizable models to accurately treat environments of varying polaritiy to produce a more accurate representation of the experimental regimen.

3.3 Application of Polarizable Models to Small Molecules

A number of small molecules have been studied using polarizable force fields. Examples include the major organic functional groups, for example, pure alkanes, alcohols, thiols, aromatic compounds, aldehydes, ketones, ethers, amines and amides; chlorinated compounds, including CCl4 and CH2Cl2, and the guanidinium ion. Many studies included the compounds in aqueous solution.

Early studies of small organic molecules were limited to determination of electrostatic properties using polarizable methods. No et al. applied an electronegativity equalization method to determine net atomic charges of 25 small organic molecules including alcohols, ethers, esters, aldehydes, ketones, thiols, thioethers, secondary amines and alkanes [269], and to ionic and aromatic compounds [270].

An early condensed phase study of polarizable alkanes was presented by Rick and Berne [184]. The manuscript reported the free energy of methane association in water using a polarizable fluctuating charge model. Two previous studies only included polarizability on the water molecules [183,271]. The hydrophobic interaction was more recently analyzed using polarizable models by Rick [272], who calculated the heat capacity change for methane pair aggregation. Chelli and co-workers [126] applied the fluctuating point charge model (FQ) and the atom-atom charge transfer model (AACT), fitted to the polarizability of small alkanes and polyenes, to larger homologues. The AACT scheme was found to perform better on alkanes of any length and conformation. The AACT scheme also satisfactorily reproduced the polarization response for highly conjugated systems.

A number of additional studies of alkanes using polarizable models have been reported. Bret et al. developed force field parameters for methane in the framework of the chemical potential equalization model [97]. Studies of methane clathrate hydrates were performed by English and MacElroy [273] using flexible and rigid polarizable and nonpolarizable water and flexible and rigid methane models. Parametrization and testing of the ABEEM/MM fluctuating charge force field for alkanes was described by Zhang and Yang [274]. Borodin and Smith [275] reported the development of many-body polarizable force fields for ether, alkane and carbonate-based solvents. Alkane parameters were also developed for the polarizable methods included in CHARMM. MacKerell and co-workers presented a systematic study of a Drude oscillator-based model of alkanes [276], calculating bulk thermodynamic, structural, dielectric, and aqueous solvation properties. Patel and Brooks [277] presented a study on a polarizable model of hexane in the framework of the fluctuating charge method, focusing on bulk liquid phase properties and analysis of the hexane–water interface. Recently, this work was extended to include longer alkanes by Davis et al. [278]. Development of a polarizable intermolecular potential function (PIPF) for liquid amides and alkanes has been reported [36]. Another application reported by Jalkanen and Zerbetto [279] studied the adsorption of organics on a silver surface using an Embedded Atom Model for the metal, a standard bonded potentials for the organics, and a combination of the charge equilibration model and the Morse potential for their electrostatic and nonbonding interactions.

Alcohols were the subject of several studies using polarizable force fields. One of the earliest studies of a nonadditive MD simulation of a pure alcohol was reported by Caldwell and Kollman, who studied the structure and properties of pure methanol [280]. This was followed by calculation of the aqueous solution free energy of methanol [281]. Chelli and co-workers applied the chemical potential equalization method to calculate the optical spectra of liquid methanol [96,98] and investigated the polarization response of methanol by polarizable force field and density functional theory calculations [282]. The method to model polarization developed by Ferenczy and Reynolds [283,284] was also applied to methanol complexes [285]. Methanol was also used in the development of the induced dipole method of Berne, Friesner and co-workers [286]. A polarizable model for simulation of liquid methanol was developed using the Charge-on-Spring (COS) technique and is compatible with the COS/G2 water model. The model was used to study the thermodynamic, dynamic, structural, and dielectric properties of liquid methanol and of a methanol-water mixture [287]. MC simulations of liquid methanol have also been reported using a potential including polarizability, nonadditivity, and intramolecular relaxation [288]. The classical Drude oscillator model implemented in CHARMM was used to study water-ethanol mixtures by Noskov et al. [68]. Interestingly, although the water and ethanol models were parametrized separately to reproduce their respective vaporization enthalpies, static dielectric constants, and self-diffusion constants of the pure liquids, the model was able to reproduce the energetic and dynamical properties of the mixtures accurately. Furthermore, the calculated dielectric constant for the various water-alcohol mixtures is in excellent agreement with experimental data. A revised Drude model for primary and secondary alcohols has been presented by Anisimov et al. [69]. That work indicated significant differences in alcohol-water RDFs as compared to the CHARMM additive force field, suggesting that the inclusion of polarizability alters atomic details of the interactions between these classes of molecule. Parameters for ethanol and methanol have also been developed for the FQ implementation in CHARMM [289,290]. Recently, thermodynamic and structural properties of methanol-water solutions were published using that model [291]. Polarizable force fields have also been developed for thiols and other sulfur containing compounds like thioethers and disulfides. Noteworthy is the work of Kaminski et al. [286,292] and sulfur parameters based on the Drude oscillator model are in progress (X. Zhu and A.D. MacKerell, Jr., Work in progress).

Other classes of small organic molecules studied with polarizable force fields are aromatic and heteroaromatic compounds. Stern et al. parameterized electrostatic parameters of substituted benzenes based on fluctuating charge, induced dipole, and a combined model and applied the resulting parameters to compute conformational energies of the alanine, serine and phenylalanine dipeptides [87]. The effort of Berne, Friesner and co-workers to develop a polarizable force field for small organic molecules within a fluctuating charge approach also included benzene and phenol [292]. The development of the DRF90 force field of Swart and van Duijnen relied on comparison of computed interaction energies and geometries of benzene dimers with ab initio results [29]. Lopes et al. published a study of aromatic compounds using the classical Drude formalism implemented in CHARMM. Benzene dimer interaction energies and geometries were considered and thermodynamic and transport properties in condensed phase were computed and compared with experimental values [70]. Soteras et al. developed models of distributed atomic polarizabilities for the treatment of induction effects in MM simulations within the framework of the induced dipole model. Molecular polarizabilities were computed for benzene, pyridine, imidazole, indole, aniline, benzonitrile, phenol and halogenated benzenes [293]. Mayer and Astrand developed a charge-dipole model for the static polarizability of nanostructures that include aliphatic, olephinic and aromatic systems [294]. MacKerell and co-workers have recently published force field parameters for pyridine, pyrimidine, imidazole, pyrrole, indole and purine [10], an effort that will lay the ground work for the development of a nucleic acids force field. That effort relied heavily on the reproduction of a variety of experimental condensed phase properties including pure solvents, crystals and aqueous solvation. Use of multiple types of condensed phase data is important is it increases the number of types of molecules that can be optimized using experimental thermodynamic data and the types of environments that can be considered during the force field optimization.

Ethers, ketones and aldehydes are among the most studied molecules using polarizable force field methods. Shirts and Stolworthy [295] analysis of a crown ether (18-crown-6) showed that the electrostatic term is the largest contributor to the conformational energy and discussed the desirability of using a polarizable method, such as the charge equilibration algorithms, to include these effects in MM and MD calculations. During development of the FQ method by Berne and co-workers investigations of the aqueous solvation and reoganization energy of other molecules, notably formaldehyde, were performed [296,297]. The development of new schemes of including polarization in classical force fields often used water-formaldehyde complexes as a source of target data for the parametrization. The methods to model polarization developed by Ferenczy and Reynolds [283,284] and Krimm and co-workers included formaldehyde-water complexes [285,298]. Borodin and Smith developed classical polarizable force fields for several molecules including polyethers, ketones, and linear and cyclic carbonates on the basis of QM dimer energies of model compounds and empirical thermodynamic liquid-state properties [275,299-305].

The first studies of amines and amides using polarizable force fields were the hydration calculations by Krogh-Jespersen, Levy and co-workers [306], and of amine hydration by Kollman and co-workers [307]. Despite the very different parametrizations, inclusion of polarizability substantially improved the reproduction of the experimental free energies of aqueous solvation in both studies. Amides, in particular N-methylacetamide (NMA), are extremely important in the development of polarizable force fields for proteins since they constitute the smallest unit representative of the poly-peptide backbone. Recently, MacKerell, Roux and co-workers parametrization of NMA within the context of the classical Drude polarizable method was the first force field to reproduce the large dielectric constant of liquid NMA [84]. Polarizable models of amide compounds that have two (acetamide) and zero (N,N-dimethyl acetamide) polar hydrogen-bond donor atoms were also investigated. In those studies it was shown that a proper representation of both the magnitude and direction of the molecular polarizability tensor, made possible by the use of atom-based Thole damping factors, was essential to obtain the correct dielectric response.

The ability of polarizable force fields to accurately reproduce dielectric constants deserves additional discussion. Although the dielectric is a macroscopic property of bulk system, it has a critical impact on microscopic interactions within the system. This can be qualitatively illustrated by considering the familiar Born model of solvation, which shows that the solvation free energy of an ion of charge q and radius R is given by q2/(2R) (1/ε–1). This expression shows that correct estimation of ε is essential to obtain the proper solvation thermodynamics. The treatment of induced electronic polarization becomes of particular importance in the case of the low-dielectric alkanes. The correct value of ε is approximately 2, which is only possible to attain in polarizable models as additive models with fixed partial charges yield values that are approximately equal to 1 [85]. Obviously, going from a value of 1 to 2 drastically impacts solvation energies, as shown in the context of the Born approximation, such that the ability of force fields to model relative solvation in complex simulations, such as lipid bilayers, will be drastically effected. For example, the neglect of polarization of the hydrocarbon core of lipid membranes has been shown to have great practical consequences in computational studies of ion channels [308,309]. However, attaining the correct dielectric behavior appears to not be trival. In our own hands, the assumption that the gas phase polarizabilities would be applicable to the condensed phase was shown to be incorrect on the first molecule studied, water, as discussed above. This lead to the development of an approach whereby the polarizability of a molecule is considered a free parameter during optimization, with the primary target data being reproduction of the dielectric constant of the corresponding pure solvent. To date, our efforts indicate that the appropriate polarizability depends on the class of molecule under study. As quantified in terms of scaling of the gas phase polarizability, we have empirically determined scaling factors ranging from 0.7 for water, alcohols and sulfur containing species, 0.85 for aromatics, N-containing heteroaromatics and ethers (C. Baker and A.D. MacKerell, Jr. Work in progress) and 1.0 for amides. Given the role of electrostatics in a variety of complex phenomena involving biological systems (eg. pKas, reduction potentials) and the contribution of the dielectric to proper treatment of electrostatics, careful consideration of this important term is central to successful development of polarizable force fields for biological molecules.

QM studies also indicate that polarizability scaling for the condensed phase may be necessary. Based on studies of water clusters, it was suggested that in the condensed phase, polarization is lower than in gas phase because of the energetic cost arising from Pauli’s exclusion principle due to the overlap of neighbouring electronic charge distributions [310]. A recent study by Schropp and Tavan [311] on the polarization of a single QM water molecule within a MM described bulk phase also concluded that the gas phase experimental polarizability cannot be used in molecular simulations but must be reduced to an effective polarizability. It was argued that in the liquid phase the electric field, [left angle bracket]E[right angle bracket], in the excluded volume of each water molecule is strongly inhomogeneous such that the electric field at the position of the oxygen (or hydrogens), E(rO), is not appropriate for calculation of the molecular polarizability. Since [left angle bracket]E[right angle bracket] is smaller than E(rO) and it is necessary to use E(rO) in MD simulations because of computational efficiency, results that E(rO) needs to be scaled to match [left angle bracket]E[right angle bracket]. Remarkably, the scaling factor proposed by Schropp and Tavan (0.68) is close to the empirical value of approximately 0.7 proposed by Lamoureux, MacKerell and Roux [76,77]. Another recent study using semiempirical methods on model compounds representative of phospholipids also indicated that the polarizability of the head group decreased in the presence of water, suggesting the effect is due to making “electrons in hydrogen bonds to be more bound” [312]. However, the AMEOBA water model, which has been developed with inducible point-dipoles on the oxygen as well as on the hydrogen atoms without any scaling, only slightly overestimates the dielectric of bulk water under ambient conditions (see above), indicating that the extent of scaling may also be dependent on the method used to treat polarizability. Thus, while both QM as well as empirical approaches based on reproduction of condensed phase properties indicate the need for polarizability scaling for some classes of molecules, the cause of the effect is still a matter of debate.

In general, the simulation studies of small molecules using polarizable models have shown that the extension of empirical force fields to include polarization is indeed feasible. In many, but not all cases the polarizable models have lead to improved agreement with experiment. In addition, the atomic details of the interactions between components in condensed phase simulations have been observed to differ in polarizable models as compared to additive models [70,85,313]. Thus, it appears that polarizable models will lead to a more accurate picture of the atomic details of condensed phases; however, it should emphasized that careful optimization methods, including optimization of the LJ parameters, followed by rigorous validation of the models is essential to assure that those models are yielding atomic pictures representative of the experimental regimen.

3.4 Application of Polarizable Force Fields to Proteins, Nucleic Acids and Lipid Bilayers

The ultimate goal of polarizable force fields with biological relevance is the development of fully usable, high quality force fields applicable to simulations of large biomolecules: proteins, DNA/RNA, lipid bilayers and carbohydrates. While the development of such force fields that have been fully optimized is still in its infancy, very early studies that applied polarizable models to proteins in MM calculations should be noted. These include a study of lysozyme by Warshel and Levitt [314] who simulated the electrostatic environment by a polarizable force field based on induced dipoles and represented the effect of the surrounding solvent by a microscopic dielectric model. Similar approaches were used for other systems [138,315,316]. While these studies only involved single point calculations (ie. calculation of the polarization response on a single protein conformation), they emphasize that early workers were well aware of the importance of this term in theoretical studies of macromolecules. And given that it has taken over 25 years since those seminal works to start to systematically apply polarizable models to macromolecules, it is clear that the technical hurdles to the implementation and development of polarizable models have and will continue to be large. Only recently has there was a surge of publications of large molecules indicating that many of the polarizable force fields being developed in the past 10 years are nearing completion [317]. At the time of writing, many studies are still focused on validating the various force fields that have been developed over the years. However, fully featured studies that address specific research have already been published.

One of the first modern applications of a polarizable force field to a protein was on crambin using the fluctuating charge method interfaced with the UFF and AMBER force fields [111]. The polarizable charges were found to give more realistic charge redistribution between amino acids in the protein. Berne and co-workers used a combination of permanent and inducible point dipoles with fluctuating and fixed charges to simulate bovine pancreatic trypsin inhibitor (BPTI) in water with two commonly used water models TlP4P-FQ and RPOL. The simulated structures remain within 1 Å of the experimental crystal structure for the 2 ns duration of the simulations [318,319]. The extent of deviation of the structure was similar to that obtain with the OPLS all-atom additive force field. More recently, Liang and Walsh studied aqueous solvation of carboxylate groups present in the glycine zwitterion and the dipeptide aspartylalanine using the AMOEBA force field. Results were compared with Car-Parrinello MD data and additive force fields. The polarizable force field yields carboxylate solvation properties in very good agreement with CPMD results, agreement that was significantly closer than that obtained from traditional force fields [320].

Llinas et al. performed structural studies of human alkaline phosphatase using the TCPEp (topological and classical polarization effects for proteins) force field [321]. The enzyme possesses 4 metal binding sites, two for Zn2+, one for Mg2+ and one Ca2+. In this study Ca2+ was replaced by Sr2+, both showing similar interaction energies at the calcium-binding site. Only at high doses of strontium, comparable to those found for calcium, can strontium substitute for calcium. Since osteomalacia is observed after ingestion of high doses of strontium, alkaline phosphatase is likely to be one of the targets of strontium, and thus the results support the suggestion that the enzyme may be involved in this disease.

In a step forward towards a polarizable force field for proteins Wang et al. optimized the the AMBER polarizable model parameters adjusting the phi and psi torsion angles of the protein backbone by fitting to the QM energies of the important regions: beta, P-II, alpha(R), and alpha(L) regions [322]. Performance of the force field was analysed by comparison of energies against QM data and by the replica exchange molecular dynamics simulations of short polyalanine peptides in water. The populations in these three regions were found to be in qualitative agreement with the NMR and CD experimental results.

The ABEEM/MM method has been tested in studies of peptides and proteins. The first of those studies is a conformational analysis of a peptide in 2006 [323], followed by studies of trypsin inhibitors [324], BPTI in aqueous solution and conformational studies of alpha-conotoxin GI [325]. A study on the geometry of the heme prosthetic group was also published in 2008 [326].

Extension of polarizable force fields to studies of other classes of biomolecules has also been published. Allen and co-workers compared newly developed coarse-grained models, with the additive CHARMM27 force field and the newly developed Drude force field on the energetics of membrane-arginine interactions [327]. An innovative application of polarizable force fields was presented by Masella and co-workers, who combined a polarizable force field and a coarse-grained polarizable solvent with application to a long simulation of bovine pancreatic trypsin inhibitor [328]. The method was found to be much faster than traditional approaches where solvent is treated atomistically. Specialized applications included application of a polarizable force field to calculation of protein-ligand binding energies [329] and calculation of the infrared spectra of the polypeptide backbone [330].

The first application of polarizability to an MD simulation of DNA using the Drude oscillator model was performed by MacKerell and co-workers [82] in combination with the CHARMM27 all-atom nucleic acid force field [34]. The system included a DNA octamer and a full solvent representation based on the SWM4-DP model with sodium counterions. The simulation system was shown to be stable, ultimately for 5 ns, though the RMS difference with respect to canonical B form DNA continually increased, well beyond that of a simulation using the CHARMM27 additive model. This deviation emphasizes the importance of full optimization of a force field when the electrostatic model is changed. This is due to the electronic term affecting geometries, vibrations and conformational energies as well as interactions with the environment. More recently, Sagui and co-workers performed molecular dynamics simulations of a DNA decamer using both additive atomic point charge and polarizable force fields [331-333]. Results showed the polarizable model to yield properties similar to that of the additive FF, which may be associated with the underpolarization of the induced dipole model used in that study.

An updated overview of current DNA modeling with ab initio (Hartree-Fock, density functional theory, and tight binding approximations) and empirical methods including polarizable methods has been published by Cozmuta and Mehrez [334]. The authors extensively review the literature until 2007. Another interesting review worth mentioning discusses the merits and limitations of modeling methods available for guanine quadruplex (G-DNA) molecules [335,336]. Although not focused on polarizable force fields it is an important source of information for the computational chemist/biologist working in the field. The review discusses the relation of simulation results to experimental techniques and the significance of those relationships. Aspects such as: pair-additive approximation of the empirical force fields, sampling limitations due to limited simulation times, accuracy of description of base stacking, H-bonding, sugar-phosphate backbone and ions by force fields, are addressed.

4. Summary

Empirical force fields that explicitly treat electronic polarizability have been known for over 30 years. However, due to technical challenges associated with their implementation, difficulties in the optimization of sufficiently accurate models, and computational demands it is only in recent years that they have started to be increasingly considered. Over the last 15 years a large number of studies using polarizable models of water, ions and of small molecules have been published. The most notable insights from these efforts are in the area of ion solvation and at interfaces and it is evident that atomistic details of intermolecular interactions differ between polarizable and additive force fields. With respect to macromolecules, a number of works have been published over the last several years with emphasis on the capabilities and accuracy of polarizable models, particularly of proteins. While the application of polarizable models to macromolcular systems has been slow to get out of the gate, it is clear that over the next decade highly optimized polarizable force fields for these systems will become available and widely used. We eagerly look forward to the novel insights what will be obtained from these studies.

Acknowledgements

Financial support from the NIH (GM 51501 to A.D.M. and GM 072558 to B.R. and A.D.M.) is acknowledged.

BIBLIOGRAPHY

1. Jorgensen WL, Chandrasekhar J, Madura JD, et al. J Chem Phys. 1983;79:926.
2. Cornell WD, Cieplak P, Bayly CI, et al. J Am Chem Soc. 1995;117:5179.
3. MacKerell AD, Bashford D, Bellott M, et al. J Phys Chem B. 1998;102:3586. [PubMed]
4. MacKerell AD, Jr., Brooks B, Brooks CL, III, et al. In: Encyclopedia of Computational Chemistry. Schleyer PvR, et al., editors. John Wiley & Sons; Chichester: 1998.
5. Scott WRP, Hunenberger PH, Tironi IG, et al. J Phys Chem A. 1999;103:3596.
6. Jorgensen WL, Tiradorives J. J Am Chem Soc. 1988;110:1657.
7. Halgren TA, Damm W. Curr Opin Struct Biol. 2001;11:236. [PubMed]
8. Rick SW, Stuart SJ. In: Reviews in Computational Chemistry. Lipkowitz KB, Boyd DB, editors. Wiley-VCH; Hoboken, NJ: 2002.
9. Warshel A, Kato M, Pisliakov AV. J Chem Theor Comput. 2007;3:2034. [PubMed]
10. Lopes PEM, Lamoureux G, MacKerell AD. J Comp Chem. 2009 In Press.
11. Bayly CI, Cieplak P, Cornell WD, et al. J Phys Chem. 1993;97:10269.
12. Mackerell AD. J Comp Chem. 2004;25:1584. [PubMed]
13. Dyke TR, Muenter JS. J Chem Phys. 1973;59:3125.
14. Gregory JK, Clary DC, Liu K, et al. Science. 1997;275:814. [PubMed]
15. Sprik M. J Phys Chem. 1991;95:2283.
16. Wallqvist A, Mountain R. In: Reviews in Computational Chemistry. Lipkowitz KB, Boyd DB, editors. Wiley-VCH; New York: 1999.
17. Soetens JC, Costa M, Millot C. Mol Phys. 1998;94:577.
18. Silvestrelli PL, Parrinello M. J Chem Phys. 1999;111:3572.
19. Silvestrelli PL, Parrinello M. Phys Rev Lett. 1999;82:3308.
20. Silvestrelli PL, Parrinello M. Phys Rev Lett. 1999;82:5415.
21. Badyal YS, Saboungi ML, Price DL, et al. J Chem Phys. 2000;112:9206.
22. Gubskaya AV, Kusalik PG. J Chem Phys. 2002;117:5290.
23. Kuo IF, Tobias DJ. J Phys Chem B. 2001;105:5827.
24. Jorgensen WL. J Chem Theory Comp. 2007;3:1877. [PubMed]
25. Liu YP, Kim K, Berne BJ, et al. J Chem Phys. 1998;108:4739.
26. Maple JR, Cao YX, Damm WG, et al. J Chem Theor Comput. 2005;1:694. [PubMed]
27. Patel S, Brooks CL. J Comp Chem. 2004;25:1. [PubMed]
28. Patel S, Mackerell AD, Brooks CL. J Comp Chem. 2004;25:1504. [PubMed]
29. Swart M, van Duijnen PT. Mol Simul. 2006;32:471.
30. Ma BY, Lii JH, Allinger NL. J Comp Chem. 2000;21:813.
31. Thole BT. Chem Phys. 1981;59:341.
32. Applequist J, Carl JR, Fung K-K. J Am Chem Soc. 1972;94:2952.
33. Kuyper LF, Hunter RN, Ashton D, et al. J Phys Chem. 1991;95:6661.
34. Foloppe N, MacKerell AD. J Comp Chem. 2000;21:86.
35. Stern HA, Rittner F, Berne BJ, et al. J Chem Phys. 2001;115:2237.
36. Xie WS, Pu JZ, MacKerell AD, Gao J. J Chem Theor Comput. 2007;3:1878.
37. Gao JL, Habibollazadeh D, Shao L. J Phys Chem. 1995;99:16460.
38. Gao JL, Pavelites JJ, Habibollazadeh D. J Phys Chem. 1996;100:2689.
39. van Duijnen PT, Swart M. J Phys Chem A. 1998;102:2399.
40. Xie W, Pu J, Gao J. J Phys Chem A. 2009;113:2109. [PMC free article] [PubMed]
41. Ferenczy GG, Reynolds CA. J Phys Chem A. 2001;105:11470.
42. Ren PY, Ponder JW. J Comp Chem. 2002;23:1497. [PubMed]
43. Grossfield A, Ren PY, Ponder JW. J Am Chem Soc. 2003;125:15671. [PubMed]
44. Grossfield A, Ren PY, Ponder JW. Biophys J. 2003;84:94A.
45. Ren PY, Ponder JW. J Phys Chem B. 2003;107:5933.
46. Ren PY, Ponder JW. J Phys Chem B. 2004;108:13427.
47. Grossfield A. J Chem Phys. 2005;122 [PubMed]
48. Jiao D, King C, Grossfield A, et al. J Phys Chem B. 2006;110:18553. [PubMed]
49. Rasmussen TD, Ren PY, Ponder JW, et al. Int J Quantum Chem. 2007;107:1390.
50. Drude P. The Theory of Optics. Kessinger Publishing Company; 1902. 2008.
51. London F. Trans Faraday Soc. 1937;33:8b.
52. Bade WL. J Chem Phys. 1957;27:1280.
53. Bade WL, Kirkwood JG. J Chem Phys. 1957;27:1284.
54. Bade WL. J Chem Phys. 1958;28:282.
55. Whitfield TW, Martyna GJ. Chem Phys Lett. 2006;424:409.
56. Amos AT. Int J Quantum Chem. 1996;60:67.
57. Wang F, Jordan KD. J Chem Phys. 2002;116:6973.
58. Dick BG, Overhauser AW. Phys Rev. 1958;112:90.
59. Hanlon JE, Lawson AW. Phys Rev. 1959;113:472.
60. Jacucci G, McDonald IR, Singer K. Phys Lett A. 1974;50:141.
61. Lindan PJD, Gillan MJ. J Phys: Condens Matter. 1993;5:1019.
62. Mitchell PJ, Fincham D. J Phys: Condens Matter. 1993;5:1031.
63. Lindan PJD. Mol Simul. 1995;14:303.
64. Hoye JS, Stell G. J Chem Phys. 1980;73:461.
65. Pratt LR. Mol Phys. 1980;40:347.
66. Lado F. J Chem Phys. 1997;106:4707.
67. Cao J, Berne BJ. J Chem Phys. 1993;99:2213.
68. Noskov SY, Lamoureux G, Roux B. J Phys Chem B. 2005;109:6705. [PubMed]
69. Anisimov VM, Vorobyov IV, Roux B, et al. J Chem Theor Comput. 2007;3:1927. [PMC free article] [PubMed]
70. Lopes PEM, Lamoureux G, Roux B, et al. J Phys Chem B. 2007;111:2873. [PMC free article] [PubMed]
71. Saint-Martin H, Medina-Llanos C, Ortega-Blake I. J Chem Phys. 1990;93:6448.
72. de Leeuw NH, Parker SC. Phys Rev B. 1998;58:13901.
73. Saint-Martin H, Hernandez-Cobos J, Bernal-Uruchurtu MI, et al. J Chem Phys. 2000;113:10899.
74. van Maaren PJ, van der Spoel D. J Phys Chem B. 2001;105:2618.
75. Yu HB, Hansson T, van Gunsteren WF. J Chem Phys. 2003;118:221.
76. Lamoureux G, MacKerell AD, Roux B. J Chem Phys. 2003;119:5185.
77. Lamoureux G, Harder E, Vorobyov IV, et al. Chem Phys Lett. 2006;418:245.
78. Lamoureux G, Roux B. J Phys Chem B. 2006;110:3308. [PubMed]
79. Whitfield TW, Varma S, Harder E, et al. J Chem Theor Comput. 2007;3:2068. [PMC free article] [PubMed]
80. Lu ZY, Zhang YK. J Chem Theor Comput. 2008;4:1237.
81. Lamoureux G, Roux B. J Chem Phys. 2003;119:3025.
82. Anisimov VM, Lamoureux G, Vorobyov IV, et al. J Chem Theor Comput. 2005;1:153. [PubMed]
83. Harder E, Anisimov VM, Vorobyov IV, et al. J Chem Theor Comput. 2006;2:1587. [PubMed]
84. Harder E, Anisimov VM, Whitfield TW, et al. J Phys Chem B. 2008;112:3509. [PMC free article] [PubMed]
85. Vorobyov IV, Anisimov VM, MacKerell AD. J Phys Chem B. 2005;109:18988. [PubMed]
86. Vorobyov I, Anisimov VM, Greene S, et al. J Chem Theor Comput. 2007;3:1120. [PubMed]
87. Stern HA, Kaminski GA, Banks JL, et al. J Phys Chem B. 1999;103:4730.
88. Banks JL, Kaminski GA, Zhou RH, et al. J Chem Phys. 1999;110:741.
89. Kaminski GA, Stern HA, Berne BJ, et al. J Comp Chem. 2002;23:1515. [PMC free article] [PubMed]
90. Yin DX, Mackerell AD. J Comp Chem. 1998;19:334.
91. Rick SW, Stuart SJ, Berne BJ. J Chem Phys. 1994;101:6141.
92. Nalewajski RF. Int J Quantum Chem. 1991;40:265.
93. Chelli R, Procacci P. J Chem Phys. 2002;117:9175.
94. Itskowitz P, Berkowitz ML. J Phys Chem A. 1997;101:5687.
95. Nalewajski RF. Pol J Chem. 1998;72:1763.
96. Chelli R, Ciabatti S, Cardini G, et al. J Chem Phys. 1999;111:4218.
97. Bret C, Field MJ, Hemmingsen L. Mol Phys. 2000;98:751.
98. Chelli R, Ciabatti S, Cardini G, et al. J Chem Phys. 2000;112:5515.
99. Nalewajski RF. Int J Quantum Chem. 2000;78:168.
100. Llanta E, Ando K, Rey R. J Phys Chem B. 2001;105:7783.
101. York DM, Yang WT. J Chem Phys. 1996;104:159.
102. Smith PE. J Phys Chem B. 2004;108:16271.
103. Chelli R, Barducci A, Bellucci L, et al. J Chem Phys. 2005;123 [PubMed]
104. Medeiros M. Theor Chem Acc. 2005;113:178.
105. Piquemal JP, Chelli R, Procacci P, et al. J Phys Chem A. 2007;111:8170. [PubMed]
106. Warren GL, Davis JE, Patel S. J Chem Phys. 2008;128 [PubMed]
107. Zhang Y, Lin H. J Chem Theor Comput. 2008;4:414. [PubMed]
108. Rappé AK, Goddard WA. J Phys Chem. 1991;95:3358.
109. Kitao O, Ogawa T. Mol Phys. 2003;101:3.
110. Nistor RA, Polihronov JG, Muser MH, et al. J Chem Phys. 2006;125 [PubMed]
111. Ogawa T, Kurita N, Sekino H, et al. Chem Phys Lett. 2004;397:382.
112. Sefcik J, Demiralp E, Cagin T, et al. J Comp Chem. 2002;23:1507. [PubMed]
113. Tanaka M, Siehl HU. Chem Phys Lett. 2008;457:263.
114. Brodersen S, Wilke S, Leusen FJJ, et al. Cryst Growth Des. 2005;5:925.
115. Chen B, Xing JH, Siepmann JI. J Phys Chem B. 2000;104:2391.
116. Stuart SJ, Berne BJ. J Phys Chem. 1996;100:11934.
117. Stuart SJ, Berne BJ. J Phys Chem A. 1999;103:10300.
118. Rick SW, Berne BJ. J Am Chem Soc. 1996;118:672.
119. Toufar H, Baekelandt BG, Janssens GOA, et al. J Phys Chem. 1995;99:13876.
120. Parr RG, Pearson RG. J Am Chem Soc. 1983;105:7512.
121. Kitaura K, Morokuma K. Int J Quantum Chem. 1976;10:325.
122. Weinhold F. J Mol Struct: THEOCHEM. 1997;398-399:181.
123. van der Vaart A, Merz KM. J Am Chem Soc. 1999;121:9182.
124. Korchowiec J, Uchimaru T. J Chem Phys. 2000;112:1623.
125. Jeziorski B, Moszynski R, Szalewicz K. Chem Rev. 1994;94:1887.
126. Chelli R, Procacci P, Righini R, et al. J Chem Phys. 1999;111:8569.
127. Yang ZZ, Wang CS. J Phys Chem A. 1997;101:6315.
128. Wang CS, Li SM, Yang ZZ. J Mol Struct: THEOCHEM. 1998;430:191.
129. Wang CS, Yang ZZ. J Chem Phys. 1999;110:6189.
130. Cong Y, Yang ZZ. Chem Phys Lett. 2000;316:324.
131. Yang ZZ, Wang CS. J Theor Comput Chem. 2003;2:273.
132. Yang ZZ, Wu Y, Zhao DX. J Chem Phys. 2004;120:2541. [PubMed]
133. Wu Y, Yang ZZ. J Phys Chem A. 2004;108:7563.
134. Parr RG, Yang W. Density-Functional Theory of Atoms and Molecules. Oxford University Press; Oxford: 1989.
135. Field MJ. Mol Phys. 1997;91:835.
136. Gao JL. J Phys Chem B. 1997;101:657.
137. Gao JL. J Chem Phys. 1998;109:2346.
138. van Belle D, Couplet I, Prevost M, et al. J Mol Biol. 1987;198:721. [PubMed]
139. Svishchev IM, Kusalik PG, Wang J, et al. J Chem Phys. 1996;105:4742.
140. Bernardo DN, Ding YB, Kroghjespersen K, et al. J Phys Chem. 1994;98:4180.
141. Medeiros M, Costas ME. J Chem Phys. 1997;107:2012.
142. Campbell T, Kalia RK, Nakano A, et al. Phys Rev Lett. 1999;82:4866.
143. Car R, Parrinello M. Phys Rev Lett. 1985;55:2471. [PubMed]
144. van Belle D, Wodak SJ. Comput Phys Comm. 1995;91:253.
145. Tuckerman ME, Martyna GJ. J Phys Chem B. 2000;104:159.
146. Tuckerman ME, Martyna GJ. J Phys Chem B. 2001;105:7598.
147. van Belle D, Froeyen M, Lippens G, et al. Mol Phys. 1992;77:239.
148. Sprik M, Klein ML. J Chem Phys. 1988;89:7556.
149. Stillinger FH, David CW. J Chem Phys. 1978;69:1473.
150. Lybrand TP, Kollman PA. J Chem Phys. 1985;83:2923.
151. Barnes P, Finney JL, Nicholas JD, et al. Nature. 1979;282:459.
152. Rullmann JAC, van Duijnen PT. Mol Phys. 1988;63:451.
153. Ahlström P, Wallqvist A, Engström S, et al. Mol Phys. 1989;68:563.
154. Cieplak P, Kollman P, Lybrand T. J Chem Phys. 1990;92:6755.
155. Niesar U, Corongiu G, Clementi E, et al. J Phys Chem. 1990;94:7949.
156. Kuwajima S, Warshel A. J Phys Chem. 1990;94:460.
157. Caldwell J, Dang LX, Kollman PA. J Am Chem Soc. 1990;112:9144.
158. Dang LX. J Chem Phys. 1992;97:2659.
159. Kozack RE, Jordan PC. J Chem Phys. 1992;96:3120.
160. Wallqvist A, Berne BJ. J Phys Chem. 1993;97:13841.
161. Brodholt J, Sampoli M, Vallauri R. Mol Phys. 1995;85:81.
162. Brodholt J, Sampoli M, Vallauri R. Mol Phys. 1995;86:149.
163. Chialvo AA, Cummings PT. J Chem Phys. 1996;105:8274.
164. Dang LX, Chang TM. J Chem Phys. 1997;106:8149.
165. Wallqvist A, Ahlström P, Karlström G. J Phys Chem. 1990;94:1649.
166. Halley JW, Rustad JR, Rahman A. J Chem Phys. 1993;98:4110.
167. Corongiu G. Int J Quantum Chem. 1992;42:1209.
168. Zhu SB, Singh S, Robinson GW. J Chem Phys. 1991;95:2791.
169. Zhu SB, Yao S, Zhu JB, et al. J Phys Chem. 1991;95:6211.
170. Feyereisen MW, Feller D, Dixon DA. J Phys Chem. 1996;100:2993.
171. Burnham CJ, Li JC, Xantheas SS, et al. J Chem Phys. 1999;110:4566.
172. Dang LX, Smith DE. J Chem Phys. 1993;99:6950.
173. Sprik M, Klein ML, Watanabe K. J Phys Chem. 1990;94:6483.
174. Xantheas SS, Dang LX. J Phys Chem. 1996;100:3989.
175. Bryce RA, Buesnel R, Hillier IH, et al. Chem Phys Lett. 1997;279:367.
176. Cabarcos OM, Weinheimer CJ, Lisy JM. J Chem Phys. 1999;110:8429.
177. Baik J, Kim J, Majumdar D, et al. J Chem Phys. 1999;110:9116.
178. Robertson WH, Diken EG, Price EA, et al. Science. 2003;299:1367. [PubMed]
179. Lightstone FC, Schwegler E, Hood RQ, et al. Chem Phys Lett. 2001;343:549.
180. Bako I, Hutter J, Palinkas G. J Chem Phys. 2002;117:9838.
181. Wallqvist A. Chem Phys Lett. 1990;165:437.
182. Motakabbir KA, Berkowitz ML. Chem Phys Lett. 1991;176:61.
183. New MH, Berne BJ. J Am Chem Soc. 1995;117:7172.
184. Rick SW, Berne BJ. J Phys Chem B. 1997;101:10488.
185. Soper AK, Phillips MG. Chem Phys. 1986;107:47.
186. Dore JC. J Mol Struct. 1991;250:193.
187. Dore JC, Blakey DM. J Mol Liq. 1995;65-6:85.
188. Soper AK. J Phys: Condens Matter. 1997;9:2717.
189. Jedlovszky P, Brodholt JP, Bruni F, et al. J Chem Phys. 1998;108:8528.
190. Soper AK. Chem Phys. 2000;258:121.
191. Nakahara M, Matubayasi N, Wakai C, et al. J Mol Liq. 2001;90:75.
192. Narten AH, Levy HA. In: Water A Comprehensive Treatise. Franks F, editor. Plenum Press; New York: 1972.
193. Page DI. In: Water A Comprehensive Treatise. Franks F, editor. Plenum Press; New York: 1972.
194. Narten AH, Thiessen WE, Blum L. Science. 1982;217:1033. [PubMed]
195. Yamanaka K, Yamaguchi T, Wakita H. J Chem Phys. 1994;101:9830.
196. Hura G, Sorenson JM, Glaeser RM, et al. J Chem Phys. 2000;113:9140.
197. Sorenson JM, Hura G, Glaeser RM, et al. J Chem Phys. 2000;113:9149.
198. Head-Gordon T, Hura G. Chem Rev. 2002;102:2651. [PubMed]
199. Soper AK, Turner J. Int J Mod Phys B. 1993;7:3049.
200. Soper AK, Bruni F, Ricci MA. J Chem Phys. 1997;106:247.
201. Mahoney MW, Jorgensen WL. J Chem Phys. 2000;112:8910.
202. Corongiu G, Clementi E. J Chem Phys. 1992;97:2030.
203. Berendsen H, Postma JPM, van Gusteren WF, et al. In: Intermolecular forces. Pullman B, editor. Reidel Publishing Company; Dordrecht, The Netherlands: 1981.
204. Jedlovszky P, Vallauri R. Mol Phys. 1999;97:1157.
205. Rick SW. J Chem Phys. 2001;114:2276.
206. Jedlovszky P, Richardi J. J Chem Phys. 1999;110:8019.
207. Mahoney MW, Jorgensen WL. J Chem Phys. 2001;115:10758.
208. Saint-Martin H, Hernandez-Cobos J, Ortega-Blake I. J Chem Phys. 2005;122
209. Lopez-Lemus J, Chapela GA, Alejandre J. J Chem Phys. 2008;128 [PubMed]
210. Yezdimer EM, Cummings PT. Mol Phys. 1999;97:993.
211. Kiyohara K, Gubbins KE, Panagiotopoulos AZ. Mol Phys. 1998;94:803.
212. Svishchev IM, Hayward TM. J Chem Phys. 1999;111:9034.
213. Berendsen HJC, Grigera JR, Straatsma TP. J Phys Chem. 1987;91:6269.
214. Panagiotopoulos AZ. Mol Phys. 1987;61:813.
215. Hill PG. J Phys Chem Ref Data. 1990;19:1233.
216. Markovich G, Perera L, Berkowitz ML, et al. J Chem Phys. 1996;105:2675.
217. Klopper W, van Duijneveldt-van de Rijdt J, van Duijneveldt FB. Phys Chem Chem Phys. 2000;2:2227.
218. Tschumper GS, Leininger ML, Hoffman BC, et al. J Chem Phys. 2002;116:690.
219. Lee HM, Suh SB, Lee JY, et al. J Chem Phys. 2000;112:9759.
220. Xantheas SS, Burnham CJ, Harrison RJ. J Chem Phys. 2002;116:1493.
221. Kaplan JH. Annu Rev Biochem. 2002;71:511. [PubMed]
222. Zhang CX, Lippard SJ. Curr Opin Chem Biol. 2003;7:481. [PubMed]
223. Aleman EA, Lamichhane R, Rueda D. Curr Opin Chem Biol. 2008;12:647. [PubMed]
224. Al-Hashimi HM, Walter NG. Curr Opin Struct Biol. 2008;18:321. [PMC free article] [PubMed]
225. Chen PR, He C. Curr Opin Chem Biol. 2008;12:214. [PubMed]
226. Manning GS. Quart Rev Biophys. 1978;11:179. [PubMed]
227. Scott WG. Curr Opin Chem Biol. 1999;3:705. [PubMed]
228. Draper DE. Rna-a Publication of the Rna Society. 2004;10:335.
229. Bartlett GJ, Borkakoti N, Thornton JM. J Mol Biol. 2003;331:829. [PubMed]
230. Harris ME, Christian EL. Curr Opin Struct Biol. 2003;13:325. [PubMed]
231. Sigel RKO, Pyle AM. Chem Rev. 2007;107:97. [PubMed]
232. Doyle DA, Cabral JM, Pfuetzner RA, et al. Science. 1998;280:69. [PubMed]
233. Dutzler R, Campbell EB, Cadene M, et al. Nature. 2002;415:287. [PubMed]
234. Jiang YX, Lee A, Chen JY, et al. Nature. 2003;423:33. [PubMed]
235. Okada T. J Chem Soc, Faraday Trans. 1991;87:3027.
236. Okada T. Anal Sci. 1998;14:469.
237. Takeda Y, Kanazawa M, Katsuta S. Anal Sci. 2000;16:929.
238. Dang LX, Rice JE, Caldwell J, et al. J Am Chem Soc. 1991;113:2481.
239. Roux B. Chem Phys Lett. 1993;212:231.
240. Roux B, Karplus M. J Comp Chem. 1995;16:690.
241. Dzidic I, Kebarle P. J Phys Chem. 1970;74:1466.
242. Marcus Y. J Chem Soc, Faraday Trans I. 1986;83:233.
243. Tissandier MD, Cowen KA, Feng WY, et al. J Phys Chem A. 1998;102:7787.
244. Åqvist J. J Phys Chem. 1990;94:8021.
245. Cheatham TE., III . In: Annual Reports in Computational Chemistry. Spellmeyer DC, editor. Elsevier; Amsterdam: 2005.
246. Beglov D, Roux B. J Chem Phys. 1994;100:9050.
247. Marchand S, Roux B. Protein Struct Funct Genet. 1998;33:265. [PubMed]
248. Rosseinsky DR. Chem Rev. 1965;65:467.
249. Kalidas C, Hefter G, Marcus Y. Chem Rev. 2000;100:819. [PubMed]
250. Grunwald E, Baughman G, Kohnstam G. J Am Chem Soc. 1960;82:5801.
251. Warren GL, Patel S. J Chem Phys. 2007;127 [PubMed]
252. Warren GL, Patel S. J Phys Chem C. 2008;112:7455.
253. Warren GL, Patel S. J Phys Chem B. 2008;112:11679. [PMC free article] [PubMed]
254. Masia M, Probst M, Rey R. J Chem Phys. 2005;123:164505. [PubMed]
255. Masia M, Probst M, Rey R. J Chem Phys. 2004;121:7362. [PubMed]
256. Li X, Yang ZZ. J Phys Chem A. 2005;109:4102. [PubMed]
257. Yang ZZ, Li X. J Phys Chem A. 2005;109:3517. [PubMed]
258. Piquemal JP, Perera L, Cisneros GA, et al. J Chem Phys. 2006;125 [PubMed]
259. Bagus PS, Illas F. J Chem Phys. 1992;96:8962.
260. Jungwirth P, Tobias DJ. Chem Rev. 2006;106:1259. [PubMed]
261. Salvador P, Curtis JE, Tobias DJ, et al. Phys Chem Chem Phys. 2003;5:3752.
262. Archontis G, Leontidis E, Andreou G. J Phys Chem B. 2005;109:17957. [PubMed]
263. Archontis G, Leontidis E. Chem Phys Lett. 2006;420:199.
264. Tuma L, Jenicek D, Jungwirth P. Chem Phys Lett. 2005;411:70.
265. Brown MA, D’Auria R, Kuo IFW, et al. Phys Chem Chem Phys. 2008;10:4778. [PubMed]
266. Wang XW, Watanabe H, Fuji M, et al. Chem Phys Lett. 2008;458:235.
267. Picalek J, Minofar B, Kolafa J, et al. Phys Chem Chem Phys. 2008;10:5765. [PubMed]
268. Hagberg D, Karlstrom G, Roos BO, et al. J Am Chem Soc. 2005;127:14250. [PubMed]
269. No KT, Grant JA, Scheraga HA. J Phys Chem. 1990;94:4732.
270. No KT, Grant JA, Jhon MS, et al. J Phys Chem. 1990;94:4740.
271. van Belle D, Wodak SJ. J Am Chem Soc. 1993;115:647.
272. Rick SW. J Phys Chem B. 2003;107:9853.
273. English NJ, MacElroy JMD. J Comp Chem. 2003;24:1569. [PubMed]
274. Zhang Q, Yang ZZ. Chem Phys Lett. 2005;403:242.
275. Borodin O, Smith GD. J Phys Chem B. 2006;110:6279. [PubMed]
276. Vorobyov IV, Anisimov VM, MacKerell ADJ. J Phys Chem B. 2005;109:18988. [PubMed]
277. Patel SA, Brooks CL. J Chem Phys. 2006;124:204706. [PubMed]
278. Davis JE, Warren GL, Patel S. J Phys Chem B. 2008;112:8298. [PubMed]
279. Jalkanen JP, Zerbetto F. J Phys Chem B. 2006;110:5595. [PubMed]
280. Caldwell JW, Kollman PA. J Phys Chem. 1995;99:6208.
281. Cieplak P, Caldwell J, Kollman P. J Comp Chem. 2001;22:1048.
282. Chelli R, Pagliai M, Procacci P, et al. J Chem Phys. 2005;122 [PubMed]
283. Ferenczy GG. J Comp Chem. 1991;12:913.
284. Winn PJ, Ferenczy GG, Reynolds CA. J Phys Chem A. 1997;101:5437.
285. Winn PJ, Ferenczy GG, Reynolds CA. J Comp Chem. 1999;20:704.
286. Kaminski GA, Friesner RA, Zhou RH. J Comp Chem. 2003;24:267. [PubMed]
287. Yu HB, Geerke DP, Liu HY, et al. J Comp Chem. 2006;27:1494. [PubMed]
288. Valdez-Gonzales M, Sanit-Martin H, Hernandez-Cobos J, et al. J Chem Phys. 2007;127
289. Patel S, Brooks CL. J Chem Phys. 2005;123:164502. [PubMed]
290. Patel S, Brooks CL. J Chem Phys. 2005;122:024508. [PubMed]
291. Zhong Y, Warren GL, Patel S. J Comp Chem. 2008;29:1142. [PMC free article] [PubMed]
292. Kaminski GA, Stern HA, Berne BJ, et al. J Phys Chem A. 2004;108:621.
293. Soteras I, Curutchet C, Bidon-Chanal A, et al. J Chem Theor Comput. 2007;3:1901.
294. Mayer A, Astrand PO. J Phys Chem A. 2008;112:1277. [PubMed]
295. Shirts RB, Stolworthy LD. J Inclusion Phenom Mol Recognit Chem. 1994;20:297.
296. Rick SW, Stuart SJ, Bader JS, et al. J Mol Liq. 1995;65-6:31.
297. Bader JS, Cortis CM, Berne BJ. J Chem Phys. 1997;106:2372.
298. Mannfors B, Palmo K, Krimm S. J Mol Struct. 2000;556:1.
299. Borodin O, Smith GD, Douglas R. J Phys Chem B. 2003;107:6824.
300. Borodin O, Smith GD. J Phys Chem B. 2006;110:6293. [PubMed]
301. Borodin O, Smith GD. J Phys Chem B. 2006;110:4971. [PubMed]
302. Borodin O, Smith GD. Macromolecules. 2006;39:1620.
303. Borodin O, Smith GD, Fan P. J Phys Chem B. 2006;110:22773. [PubMed]
304. Borodin O, Smith GD. Macromolecules. 2007;40:1252.
305. Borodin O, Smith GD, Sewell TD, et al. J Phys Chem B. 2008;112:734. [PubMed]
306. Ding YB, Bernardo DN, Kroghjespersen K, et al. J Phys Chem. 1995;99:11575.
307. Meng EC, Caldwell JW, Kollman PA. J Phys Chem. 1996;100:2367.
308. Allen TW, Andersen OS, Roux B. Proc Natl Acad Sci U S A. 2004;101:117. [PubMed]
309. Allen TW, Andersen OS, Roux B. Biophys J. 2006;90:3447. [PubMed]
310. Morita A, Kato S. J Chem Phys. 1999;110:11987.
311. Schropp B, Tavan P. J Phys Chem B. 2008;112:6233. [PubMed]
312. Botek E, Giribet C, de Azua MR, et al. J Phys Chem A. 2008;112:6992. [PubMed]
313. Anisimov VM, Roux B, MacKerell AD. Biophys J. 2007:152A.
314. Warshel A, Levitt M. J Mol Biol. 1976;103:227. [PubMed]
315. Warshel A. J Phys Chem. 1979;83:1640.
316. Warshel A. Biochemistry. 1981;20:3167. [PubMed]
317. Friesner RA. Advances In Protein Chemistry. 2006;72:79. [PubMed]
318. Harder E, Kim BC, Friesner RA, et al. J Chem Theor Comput. 2005;1:169. [PubMed]
319. Kim BC, Young T, Harder E, et al. J Phys Chem B. 2005;109:16529. [PMC free article] [PubMed]
320. Liang T, Walsh TR. Phys Chem Chem Phys. 2006;8:4410. [PubMed]
321. Llinas P, Masella M, Stigbrand T, et al. Protein Sci. 2006;15:1691. [PubMed]
322. Wang ZX, Zhang W, Wu C, et al. J Comp Chem. 2006;27:781. [PMC free article] [PubMed]
323. Yang ZZ, Zhang Q. J Comp Chem. 2006;27:1. [PubMed]
324. Guan QM, Yang ZZ. J Theor Comput Chem. 2007;6:731.
325. Jiang N, Ma J. J Phys Chem A. 2008;112:9854. [PubMed]
326. Cui BQ, Guan QM, Gong LD, et al. Chem J Chin Univ-Chinese. 2008;29:585.
327. Vorobyov I, Li LB, Allen TW. J Phys Chem B. 2008;112:9588. [PubMed]
328. Masella M, Borgis D, Cuniasse P. J Comp Chem. 2008;29:1707. [PubMed]
329. Jiao D, Golubkov PA, Darden TA, et al. Proc Natl Acad Sci U S A. 2008;105:6290. [PubMed]
330. Schultheis V, Reichold R, Schropp B, et al. J Phys Chem B. 2008;112:12217. [PubMed]
331. Baucom J, Transue T, Fuentes-Cabrera M, et al. J Chem Phys. 2004;121:6998. [PubMed]
332. Babin V, Baucom J, Darden TA, et al. Int J Quantum Chem. 2006;106:3260.
333. Babin V, Baucom J, Darden TA, et al. J Phys Chem B. 2006;110:11571. [PubMed]
334. Cozmuta I, Mehrez H. J Comput Theor Nanosci. 2007;4:349.
335. šponer J, špačková N. Methods. 2007;43:278. [PMC free article] [PubMed]
336. Stefl R, Cheatham TE, špačková N, et al. Biophys J. 2003;85:1787. [PubMed]