atomistic 
 molecular dynamics ( md ) simulations have become an integral 
 component of the tool set used to examine biomolecular systems . since 
 the first md simulation of a protein in 1977 , the time scales and sizes of computationally tractable simulations 
 have grown by orders of magnitude . systems may now be simulated on 
 the microsecond or millisecond time scale and contain over a million atoms due to increased computing power and ever - improving algorithms for 
 parallel and graphical processing unit ( gpu)-based computing . such 
 capabilities also allow for more rigorous testing of the accuracy 
 of the models used for such md simulations . quantum mechanics 
 ( qm ) provides a highly accurate representation 
 of the energetics of molecules . while some progress has been made 
 in the application of qm methods to biological macromolecules , for 
 example , hybrid quantum mechanics / molecular mechanics ( qm / mm ) methods 
 that have been utilized for decades , applying qm to large 
 biomolecular systems remains infeasible in most cases . as a result , 
 md simulations have relied upon empirical potential energy functions 
 ( so - called  force fields  ) to calculate the energies 
 of configurations and , thus , the forces on all constituent atoms required 
 to perform dynamic simulations . thus , the predictive ability of md 
 simulations relies on the accuracy of the underlying force fields 
 that describe the nature of the interactions between all components 
 in the simulation system . modeling molecular processes accurately 
 is made challenging by 
 the wide range of interactions that atoms or molecules can experience 
 in condensed - phase environments . for instance , a protein in a cell 
 may interact with other proteins with polar and nonpolar surfaces , 
 amphipathic lipid membranes , and highly charged species like ions 
 and nucleic acids ( dna and rna ) , all while diffusing through a predominantly 
 polar , aqueous environment . small molecules ( xenobiotics , metabolites , 
 etc . ) may partition across hydrophobic membranes or bind to specific 
 transporters before being delivered to intracellular receptors or 
 enzymes . thus , their dynamics must also be modeled in a manner that 
 accounts for these diverse environments . describing all of these processes 
 using computationally feasible functional forms remains a challenge 
 in the area of force field development . to date , the majority of force 
 fields for molecular simulations have been based on pairwise or additive 
 force fields in which the partial atomic charges are fixed . however , 
 this convention represents a significant approximation that limits 
 the ability of the force field to be accurate in the diverse range 
 of environments present in biological and chemical condensed - phase 
 systems . to overcome this limitation , more complex force fields that 
 include explicit electronic polarization that allow the charge distribution 
 to vary as a function of environment , the development of a force field based on the drude oscillator 
 model will be described , including the 
 theoretical background and parametrization methodology , as well as 
 applications to biomolecular systems in which explicit polarization 
 has provided unique insights into functional dynamics . traditional 
 additive force fields will be described briefly as an introduction 
 to the fundamentals of force field development . the review will conclude 
 with an outlook on the next steps in drude polarizable force field 
 development and systems to which we anticipate the inclusion of polarizability 
 will be critical . we note that the review will focus on the drude-2013 
 polarizable force field , which has its origins in the charmm additive 
 force field , with the development of the force field initiated in 
 2000 with the first polarizable water model published in 2003 . while there are efforts by other groups to develop 
 polarizable models , a comprehensive evaluation of all of these force 
 fields is beyond the scope of the present review . to provide perspective 
 on the field , we will briefly discuss efforts by other investigators 
 using alternate models in section 3.1 . most md simulations to date have used fixed - charge ( additive ) force 
 fields . a force field has two components , a functional form for calculating 
 energies , and a parameter set that populates the variables in the 
 functional form . these force fields are characterized by point charges 
 on each of the atoms , centered on the atomic nucleus . the sum of these 
 pairwise interactions gives the total electrostatic energy , hence 
 the name  additive . there are many popular atomistic 
 force fields for biomolecular simulations , the most widely used of 
 which are amber , charmm , gromos , and opls . as an example , the functional form of the charmm 
 force field is as follows:1 the calculated 
 potential energy of a configuration , u(r ) , is thus a sum of internal ( bonded ) 
 and noncovalent ( nonbonded ) interactions . bonds ( b ) , valence angles (  ) , urey  bradley ( s ) , dihedral angles ( or torsions ) (  ) , cmap backbone torsion 
 corrections ( , ) , and improper 
 dihedrals (  ) comprise the internal terms . most of these terms 
 are described by simple harmonic functions with a force constant , k , with respect to an equilibrium value ( subscript 0 ) . the 
 dihedral term describing bond rotations is not harmonic ; instead , 
 it is a periodic term characterized by a force constant , multiplicity 
 ( n ) , and phase shift (  ) . bradley 
 cross - terms describe interactions between atoms i and k ( so - called  13 interactions  ) 
 in a valence angle formed by atoms i  j  k . improper dihedrals describe out - of - plane vibrations of planar groups . the van der waals component is described by a lennard - jones 
 ( lj or 612 potential  ) form , and electrostatic 
 interactions are calculated using coulomb s law with point 
 charges qi and qj on atoms i and j , respectively . this potential energy function 
 shares many common features with those of the other force fields . while additive force fields have enjoyed many years of development 
 and success in describing macromolecular dynamics , there are inherent problems and limitations with these 
 models . additive force fields are generally parametrized using a mean - field 
 approximation , such that they respond in an average way to the different 
 environments that a molecule might experience during the course of 
 a simulation , a consideration that is especially important in terms 
 of the electrostatic parameters . a common strategy for implicitly 
 representing polarization response in additive force fields is to 
 overestimate the gas - phase dipole moments from qm calculations , typically 
 on the order of 20% or more . while this 
 approach is tractable from the perspective of parametrization and 
 computational cost , this underlying approximation has important consequences . many biological processes involve the interaction of molecules that 
 pass through or interact with the predominantly aqueous milieu as 
 well as less polar environments , such as substrate binding to a soluble 
 enzyme ( which requires at least partial desolvation prior to entering 
 the binding site ) or passage of small molecules through membranes 
 ( passing into and through the hydrophobic core of the membrane ) . the 
 electronic distribution of these molecules will change in response 
 to variations in the local electric field , but assigning fixed charges 
 prevents such a response and is inherently a deficient description 
 of the interactions in the system that limits accuracy . that is to say , the inaccuracies in additive 
 models are not necessarily due to inadequate parameters ; rather , they 
 reflect a fundamental limitation in the functional form . it is for 
 this reason that an extended description of the physical model , one 
 that includes explicit polarization , is desirable . though the present review focuses on the classical 
 drude oscillator model , it is important to briefly summarize other 
 approaches that have been developed before moving on to a more in - depth 
 discussion of the drude-2013 polarizable force field . more comprehensive 
 reviews of these polarization methods have been presented previously . given that empirical force fields are approximations of qm 
 representations , modeling induced polarization via qm as an extension 
 to the empirical energy function is an intuitive approach . the x - pol 
 potential is one such approach that is 
 related to qm / mm simulations . unlike traditional qm / mm , in which part 
 of the simulation is treated via qm , in the x - pol potential , the addition 
 of a qm layer to the mm system accounts for polarization effects . the application of a semiempirical qm method reduces the computational 
 cost relative to ab initio approaches , but the cost of such simulations 
 remains very high . as such , it is attractive to represent many - body 
 effects classically , as a direct extension or modification of the 
 mm potential . the three primary methods for treating electronic polarization 
 classically are induced point dipole , fluctuating charge , and the drude oscillator ( also called  shell  or  charge - on - a - spring  ) model . an outgrowth 
 of the drude oscillator model is the qmpff general force field , which treats polarization via the attachment of electron clouds 
 to atomic cores . rather than using a discrete particle to represent 
 electronic degrees of freedom , polarization is modeled on the basis 
 of the distance between the cloud center and the core atom . in the induced point dipole representation , each atomic site is 
 assigned a dipole that is calculated using an iterative self - consistent 
 field ( scf ) procedure ( see section 3.2.2 below ) . the most notable example of 
 an induced dipole model is the amoeba force field , which also includes multipoles up to quadrupoles in the treatment 
 of the electrostatic interactions . other induced dipole models include 
 those presented by berne , friesner , and co - workers , as well as the possim model , among 
 others . the fluctuating charge model treats the point charges on the 
 atoms 
 as dynamical variables , such that the topology can vary throughout 
 the simulation . charges on individual molecules can redistribute according 
 to the electronegativity of each atom , though the overall charge on 
 the molecule is maintained . the cheq force field for proteins , lipids , and carbohydrates is one such model that uses fluctuating charges , 
 and other fluctuating charge models have been reported . ligand 
 binding , ion solvation , and diffusion of arginine across a model lipid membrane . the main limitation with fluctuating charge 
 models is that they can not account for out - of - plane polarization in 
 conjugated systems , though , in principle , the model could be extended 
 to account for this phenomenon , e.g. , by using a semiempirical qm / mm 
 approach or in a purely mm approximation 
 by either combination with the drude oscillator model or the use of off - center virtual sites to which charge 
 can be distributed . the drude oscillator model is a direct extension 
 of additive models 
 and retains many of the pairwise features of its functional form . it simply consists of attaching an auxiliary particle carrying a charge 
 to an atom via a harmonic spring to simulate induced polarization 
 via its displacement under the influence of an electric field . in 
 early studies of ion solvation with the fluctuating charge model for 
 water , the drude formalism was used for the ions , since for a monatomic ion there is no other way to redistribute 
 charge in the fluctuating charge method . for these reasons , its implementation 
 in md codes to simulate chemical and biological systems is fairly 
 straightforward , and the approach has been used in a number of laboratories in addition 
 to our own efforts . with all polarizable models , the requirement 
 to solve the magnitude 
 of all the induced dipoles self - consistently is computationally demanding 
 and is typically a limiting factor in the efficiency of such simulations . a number of approaches have been contemplated to overcome the limitation 
 of the scf procedure . the most common method has been to attribute 
 an effective mass and kinetic energy to the degrees of freedom associated 
 with the induced dipoles via an extended lagrangian framework . this 
 concept goes back to the so - called car  parrinello md method 
 to propagate a system in density functional theory . an alternative route has been adopted by the possim force 
 field and the iamoeba model for water , whereby only a single evaluation of mutual 
 polarization is carried out . this idea of explicitly treating first - order 
 induction explicitly while accounting for higher - order contributions 
 approximately via an empirical modification of force field parameters 
 was first introduced by straatsma and mccammon . in this manner , the induced dipoles are solved in the electric 
 field of the fixed charges and multipoles in the system , and the costly 
 iterative scf procedure is avoided . this approach is computationally 
 more efficient than a full scf procedure , but it is an approximation 
 that neglects the explicit many - body effects of the other induced 
 dipoles in the system . consequently , the parametrization of the model 
 accounts for an average effect of induced dipoles , an approximation 
 analogous to traditional assumptions in additive force fields . though 
 the iamoeba approach is likely adequate for addressing polarization 
 effects in homogeneous systems such as water , it remains to be seen 
 if such an approximation is valid for heterogeneous environments or 
 in the presence of highly charged particles . for example , the coordination 
 of an ion is controlled by a delicate balance involving the interactions 
 between the induced dipoles of the ligands , an effect that is neglected 
 in a first - order approximation . a recent alternate formulation of 
 the induced dipole approach treats the polarization response using 
 perturbation theory rather than a variational approach to achieve 
 the scf condition , yielding an improvement in computational efficiency . the remainder of this review will focus 
 on the drude oscillator 
 model of induced polarization . it will cover the functional form of 
 the drude model , parametrization strategy , and development history 
 of the drude-2013 force field and its applications to biological systems 
 and will conclude with a prospective view of future studies that will 
 benefit from this polarizable model . while focusing on drude-2013 , 
 it is important to note other recent drude - based polarizable ff studies 
 include the use of multipole expansion to simulate ionic liquids , with the polarization effects yielding important 
 improvements in the translational and rotational behavior of the molecular 
 ions , and similar behavior was observed in a drude model . explicit treatment of polarizability in these 
 systems led to important observations about contributions to the dielectric 
 constant that are not possible with nonpolarizable force fields , which , 
 by virtue of their fixed charges , often retard dynamics . polarizability led to enhanced fluidity and 
 reduced long - range structure using both induced point - dipole and drude 
 models , each with comparable accuracy . the central 
 feature of the atomic force field based on the classical drude model 
 that we have been developing for molecular simulations is the inclusion 
 of electronic polarizability on all non - hydrogen atoms , as described 
 in detail in the next section . in addition , the model includes virtual 
 particles representative of lone pairs , typically located on hydrogen - bond - acceptor 
 atoms . the remainder of the functional form in the drude-2013 polarizable 
 force field shares many common elements with that of the additive 
 model shown in eq 1 . the isotropic polarizability , 
  , of a given atom is described 
 by the distribution of the total charge on the atom , q , between the core atom and its drude oscillator , which is attached 
 to the core atom via a harmonic spring , as shown schematically in figure 1a . the charges assigned 
 to the drude oscillator ( qd ) and to the 
 core atom ( qa ) are such that q 
 = qa+ qd . the atomic 
 polarizability is defined as2where kd is the 
 force constant of the drude  atom harmonic bond . in principle , qd , kd , or both can 
 be tuned to achieve an appropriate polarization response . in the present 
 drude-2013 force field , the value of kd is set to 1000 kcal mol  for all drude  atom pairs . this convention reduces the complexity 
 of parametrization by reducing the number of parameters to be fit . it follows from this description of polarization that there exists 
 some displacement , d , between the drude oscillator and 
 the core atom in response to an electric field , e , according 
 to3from which the induced atomic dipole can be 
 calculated as4 the electrostatic 
 component of the 
 system s potential energy includes a harmonic self - polarization 
 term , uself , which is calculated on the 
 basis of the harmonic bond energy between the drude oscillators and 
 their core atoms . self - polarization is thus expressed as5with kd and d defined above . as implied by eq 3 , the displacement vector is the result of 
 many - body effects induced by the electric field e resulting 
 from the configuration of the charged particles in the system . the polarizability of a given atom need not be isotropic ; that 
 is , the drude oscillator can deform differentially in all three spatial 
 dimensions . anisotropic polarization was introduced in the drude-2013 
 force field to more accurately describe functional groups acting as 
 hydrogen - bond acceptors . the anisotropic 
 form is achieved by modifying uself , which 
 is written as6 in which d1 , d2 , and d3 represent 
 the components of the drude  atom displacement vector d in a coordinate frame based on the local atomic geometry 
 ( figure 1a ) . a unit 
 vector ni , j is 
 defined between atoms i and j from 
 which each component of the displacement vector is calculated , e.g. , d1 = ni , jd along the first axis . the isotropic force 
 constant kd is expanded to a tensor , k , with off - diagonal elements set 
 to zero . consistent with the isotropic kd , the anisotropic polarization tensor is chosen such that tr[k]/3 = 1000 kcal mol  . as mentioned above , 
 atoms described by anisotropic polarization are also generally augmented 
 by  lone pairs ( virtual sites carrying negative charge ) 
 to further improve the description of electronic distribution around 
 the polarizable heavy atom ( figure 1a ) . the combination of lone pairs and anisotropic polarization 
 leads to an improved description of hydrogen bonding in polar compounds 
 and interactions with ions as a function of orientation . ( a ) the addition 
 of drude 
 oscillators to carbon and oxygen atoms via harmonic springs and the 
 subsequent distribution of charge between the atoms ( qc and qo ) and their respective 
 drude oscillators ( qdc and qdo ) . values of 11 , 22 , and 
 33 are the tensor components of the anisotropic 
 atomic polarizability along the indicated axes . dipole interactions are indicated by blue 
 and red double - headed arrows , respectively . these interactions are 
 not present in additive force fields but are treated explicitly via 
 thole screening in the drude-2013 force field , contributing to anisotropic 
 molecular polarizability . ( c ) the directional response of the drude 
 oscillators to external electric fields ( e ) is shown 
 for a simple diatomic where , due to the 12 dipole dipole 
 interactions when the electric field is perpendicular to the bond , 
 the 12 dipoles damp each other , leading to a smaller polarization 
 response compared to the electric field being parallel to the bond 
 leading to the 12 dipoles enhancing each other . one significant difference between the polarizable 
 and additive 
 functional forms is the explicit inclusion of dipole normally , in additive force fields , nonbonded 
 interactions between atoms connected by one or two bonds ( so - called 
 12 and 13 pairs ) are excluded and their effects are 
 implicitly represented in the parametrization of bond - stretching and 
 angle - bending terms . in the drude polarizable model , interactions 
 between induced dipoles ( but not charge  dipole interactions ) 
 are explicitly included for 12 and 13 atom pairs , 
 with those electrostatic interactions screened by7where rij is the distance 
 between atoms i and j , i and j are respective atomic polarizabilities , 
 and a is a damping constant . equation 7 is equivalent to a smeared charge distribution 
 described by8as originally proposed by thole . a value of 2.6 for a , similar 
 to the 2.089 suggested by thole , was found to reproduce the polarizability 
 anisotropy ratio of benzene and was 
 used initially as the default value in the drude-2013 force field . subsequently , atom - specific values were implemented , as they may be 
 assigned to fine - tune the molecular polarizability tensor9where ai and aj are respective atomistic thole factors . interactions between 
 atoms and drude oscillators separated by three or more bonds are calculated 
 according to the normal ( unscreened ) coulombic potential . it is worth pointing out that without such 12 and 13 
 dipole  dipole interactions , the polarizability tensor of molecular 
 groups is trivially diagonal by default ( i.e. , the polarizability 
 of the molecule in figure 1c would be the same , whether the field is parallel or orthogonal 
 to the chemical bond ) . in the fitting procedure used to optimize 
 the parameters , the starting 
 values for are taken from molecular polarizabilities from 
 the miller atomic hybrid polarizability ( ahp ) model . given that the drude-2013 force field does not treat hydrogen 
 atoms as polarizable , the  values of any hydrogen atoms bonded 
 to a heavy atom are added to the  value of that heavy atom . initial guesses for atomic charges and bonded parameters originated 
 from the additive charmm force field . from these initial guesses , various fitting and refinement procedures 
 have been used to obtain agreement with a wide range of target data 
 ( see section 3.3 ) . a 
 practical consideration in carrying out polarizable md simulations 
 is the efficiency of the integration scheme . there are two principal 
 methods by which simulations using the drude oscillator model can 
 be carried out . the first is the scf regime , which is an intuitive 
 extension of the born  oppenheimer approximation . that is , the 
 electronic degrees of freedom relax to their ground state immediately 
 upon any change in nuclear configuration . in an md simulation with 
 the scf approach , the positions of the drude oscillators are solved 
 via energy minimization at each integration step such that10where the index i runs over 
 all the drude this approach 
 is computationally demanding , as multiple force evaluations are required 
 at each md time step , and convergence is not guaranteed with limited 
 precision . numerical errors could also be accumulated depending on 
 the convergence criterion , initial guess , and minimization algorithms . as an alternative to the scf regime , 
 an extended lagrangian approach was developed by lamoureux 
 and roux for the drude model , in which 
 each drude oscillator is ascribed a small mass ( md ) that is subtracted from the parent atom such that the 
 total mass of the drude atom pair remained equal to the atomic 
 mass . in practice , this mass is set to 0.4 amu . since the drude oscillators 
 are not massless , their positions can be integrated on equal dynamical 
 footing with the real atoms of the system . a velocity verlet 
 algorithm is used to numerically integrate the equations of motion , 
 which use the absolute positions of the atoms and their drude oscillators 
 ( ri and rd , i , respectively ) to compute forces on the 
 drude  atom centers of mass ( ri ) and drude  atom displacements ( di ) , which are1112 to approximate 
 the scf energy surface , dual nos hoover 
 thermostats are used to regulate the temperature , t , of the real atoms ( the  physical  thermostat ) , 
 and the relative velocities of the drude oscillators are regulated 
 at some low temperature , t * , such that t *  t. the equations 
 of motion in the canonical ( nvt ) ensemble are thus 
 expressed as13141516where fr , i and fd , i are given by eqs 11 and 12 . the reduced mass of each drude  atom pair , mi , is calculated as mi = ( 1  md / mi ) . nf and nf * are the number of degrees of freedom in the physical 
 thermostat and the low - temperature drude oscillator thermostat , respectively . the forces acting on each thermostat are dynamical variables that 
 are used to calculate the friction coefficients  and 
  * on each , per eqs 15 and 16 . these friction 
 coefficients are used to scale center - of - mass and relative velocities 
 of the drude  atom pairs . isobaric ( npt ) ensemble 
 commonly used in biomolecular md simulations is also supported . the dual - thermostat extended lagrangian 
 integration scheme has 
 been implemented in charmm , namd ( as 
 langevin dynamics rather than velocity  verlet ) , chemshell qm / mm , openmm , and gromacs and is more efficient than the scf approach . the use of 
 an extended lagrangian allows for highly scalable md simulations that are only about 2 times slower than additive simulations with 
 a comparable number of atoms , assuming the same integration time step , 
 however , drude simulations are typically performed with a 1 fs time 
 step versus 2 fs for additive simulations such that the overall difference 
 between the additive and drude models is approximately a factor of 
 4 . the stability of the integration is enhanced by the use of a  hard - wall 
 constraint  to prevent polarization 
 catastrophe , which can occur when the forces acting on the drude oscillator 
 cause large displacement from its core atom . the constraint reflects 
 the drude oscillator back toward its core atom , along the drude  atom 
 bond , if it reaches a predefined limit , typically set to 0.2  . this approach is a simpler form than what was proposed by donchev 
 et al . in qmpff , which uses a distance - dependent force constant to 
 prevent polarization catastrophe . importantly , 
 the efficiency of the drude extended lagrangian integration algorithms 
 allows for simulations of biologically relevant systems on the microsecond 
 scale . the general 
 parametrization scheme for the drude-2013 polarizable force field 
 is summarized in figure 2 . the protocol relies on a strong connection between target data 
 based on qm calculations for initial parametrization and experimental , 
 condensed - phase observables for parameter validation and refinement . all qm target data are generated using structures optimized at the 
 mp2/6 - 31g(d ) level of theory for neutral molecules and mp2/6 - 31+g(d , p ) 
 for molecular ions . the fitting of electrostatic parameters has been described in detail 
 previously , but a brief discussion is appropriate 
 here . in contrast to the electrostatic component of additive force 
 fields that requires only partial charges to be fit against target 
 data , the drude polarizable model requires that the total charge on 
 the drude  atom pair , atomic polarizability (  , eq 2 ) , and thole screening 
 factors ( a , eq 7 ) all be parametrized . internal ( bonded ) terms are parametrized 
 according to crystal survey data and qm vibrational calculations ( for 
 bonds , angles , selected dihedrals , and improper dihedrals ) and 1-d 
 dihedral potential energy scans ( in the case of proper dihedrals ) 
 to obtain correct conformational energetics . potential energy surfaces 
 as a function of backbone ( , ) dihedrals are obtained 
 for model dipeptides to parametrize cmap terms ( see below ) . partial atomic charges 
 in additive force fields are often assigned 
 by targeting the qm electrostatic potential ( esp ) around the molecule 
 of interest . in the parametrization of the drude-2013 force field , 
 in addition to this unperturbed esp map , the polarization response 
 is calculated by generating multiple esp maps in the presence of individual 
 perturbing point charges with a magnitude of + 0.5 e. the qm calculations of all esp are performed on the mp2-optimized 
 geometries using the b3lyp hybrid functional with cc - pvdz or aug - cc - pvdz basis sets . b3lyp / cc - pvdz systematically underestimates dipole moments relative 
 to experimental values , whereas b3lyp / aug - cc - pvdz yields better agreement 
 and thus is preferred . for large molecules , using the augmented basis 
 set may be too computationally demanding , and for this reason the 
 b3lyp / cc - pvdz values can be used in conjunction with a well - established 
 value ( 1/0.83 ) for scaling of the initial atomic polarizability values ; 
 calculations done with the aug - cc - pvdz basis set require no initial 
 scaling . following fitting of electrostatic 
 parameters , a further scaling factor of 0.85 is generally applied 
 to the fitted atomic polarizabilities . this value represents an approximate 
 mean obtained from the empirical optimization of the scale factor 
 targeting dielectric constants for a range of small molecules ( see section 3.3.3 below ) 
 and is introduced to account for overlap of electron clouds in the 
 condensed phase that opposes induction and/or as a compensation for an inhomogeneous electric field resulting 
 from the excluded volume of water around the solute . the postfitting scaling factors are not universal , as will 
 be described in the following sections , but remain an important guideline 
 and consideration in the development of the drude polarizable force 
 field . scaling factors for specific types of functional groups are 
 discussed in section 3.3.3 below . the perturbed qm esp maps become target data 
 for a fitting function 
 that is described in detail elsewhere . during fitting , charge penalties are introduced to prevent the fitting 
 function from yielding charges with poor chemical significance , thereby 
 favoring  chemically intuitive results ; the restraint 
 is also necessary because the charge - fitting problem is underdetermined . this approach is known as restrained electrostatic potential ( resp ) 
 fitting . for esp calculations , rather 
 than constructing a cubic grid around the molecule , which is computationally 
 expensive and suffers from nonunique axes , or a random grid , which 
 impairs reproducibility , a connolly surface is generated around the molecule being parametrized onto which grid 
 points are placed . the advantages to this approach are that there 
 is shape complementarity to the molecule of interest and the total 
 number of grid points is reduced relative to a cubic grid . since the 
 surface is generated as a function of the atomic radii , multiple nonoverlapping 
 surfaces can be generated at increasing distance from the atomic centers 
 by multiplying the radii by size factors , f ( figure 3 ) . electrostatic 
 parameter fitting in the drude-2013 force field relies on these multiple 
 surfaces , which are composed of either perturbing charges or esp grid 
 points . in practice , five such surfaces are generated ; multiple perturbing 
 charge layers are needed to account for the orientation dependence 
 of the polarizability , and multiple esp surfaces are needed to calculate 
 the polarization response and the molecular dipole moment . the innermost 
 layer ( f = 2.2 ) consists of perturbing charges at 
 typical hydrogen - bonding distance and geometry , with additional charges 
 interspersed at vertices between these points . the next layer ( f = 3.0 ) is an esp surface to calculate the response to 
 these charges . to account for the orientational dependence of the 
 molecular polarizability , an additional layer of charges ( f = 4.0 ) and another esp surface ( f = 5.0 ) 
 are subsequently added . the final layer is placed at the longest distance 
 ( f = 6.0 ) from molecule , to resolve the molecular 
 dipole moment . following initial generation of the charges , polarizabilities , 
 and thole factors , the anisotropic polarizabilities are optimized , 
 typically only for hydrogen - bond acceptors . the anisotropic polarization 
 tensor is fit to esp calculated at the b3lyp / aug - cc - pvdz level of 
 theory in the presence of perturbing charges of + 0.5 e along in- and out - of - plane arcs typically around the acceptor atom . these arcs are typically 2  from the 
 acceptor atom , spanning 90180 with grid points 
 every 15. after optimization of the electrostatic parameters , 
 the internal 
 ( bonded ) parameters are optimized by targeting both qm and experimental 
 data . qm vibrational analysis is performed to determine force constants 
 for internal parameters , supplemented with dihedral potential energy 
 scans carried out for flexible torsions being parametrized . equilibrium 
 bond lengths and valence angles target available crystal survey data 
 for model compounds . in addition , for protein ( , ) conformational 
 properties , a high - level qm 2-d ( , ) potential energy 
 surface is obtained for calculation of the cmap potential . qm potential energy surfaces for dipeptides 
 are obtained at the mp2/6 - 311g(d , p)//rimp2/cbs model chemistry and serve as additional target data for protein 
 force field parametrization . the cmap corrections to ( , ) 
 dihedral cross - terms are a grid - based correction map that are an extension 
 of the force field to obtain better agreement in the energy landscape . 
 in the drude-2013 protein force field , the alanine dipeptide cmap 
 is applied to all amino acids , except glycine and proline , for which 
 specific cmaps are determined on the basis of the respective dipeptides . these simulations generally target crystal data ( molecular volume , vm ; lattice parameters ; and heat of sublimation , 
 hsub ) and properties of neat liquids 
 ( vm and heat of vaporization , hvap ) . additional target data include qm water 
 interaction energies , as well as free energies of 
 hydration . optimization of the lj parameters typically requires multiple 
 iterations and represents one of the most difficult aspects of force 
 field parametrization . the final step toward the finished model 
 consists of condensed - phase 
 simulations of macromolecules . as will be shown below , combining the 
 model compound parameters to construct polymers is nontrivial in the 
 case of the drude polarizable force field . additional effects ( particularly 
 dipole  dipole interactions ) may present themselves only in 
 macromolecules , such that multiple electrostatic models must be considered , 
 and additional dihedral and cmap refinements are often required . the 
 entire refinement protocol is iterative , as implied by figure 2 , such that adjustments at 
 any stage of the protocol require validation in all other steps , until 
 sufficient agreement is achieved for all target data . as most systems 
 of biological interest are immersed in water , the 
 natural starting point for the creation of any force field is the 
 development of a suitable water model . a central challenge in developing 
 an accurate water model is balancing structural , thermodynamic , and 
 dynamic properties . several polarizable water models have been 
 developed for use with the drude-2013 force field , the first of which 
 was a four - site model called swm4-dp , featuring a rigid configuration of hoh atoms , a negatively charged 
 virtual site (  m site  ) , and a drude oscillator attached 
 to the o atom . the m site is placed such that the gas - phase dipole 
 moment and quadrupole moments are correctly reproduced . the drude 
 oscillator and its parent o atom have charges that are equal in magnitude 
 but opposite in sign such that the total charge of the oxygen is zero . 
 in swm4-dp (  simple water model with drude polarization following 
 the development of swm4-dp , a newer model was introduced , called swm4-ndp 
 ( figure 4 ) , with the n denoting a negatively charged drude 
 oscillator ( hence  simple water model with negative drude polarization  ) . 
 given that the drude oscillator is intended to represent electronic 
 degrees of freedom , assigning it a negative charge is more intuitive 
 and may be better suited to reproducing higher - order polarization 
 effects . as such , this convention of assigning a negative charge to 
 all drude oscillators has been used throughout the development of 
 the drude-2013 force field . swm4-ndp dimer , extracted from a simulation 
 of bulk water , with 
 intermolecular and intramolecular geometric properties labeled . atoms 
 are colored by element ( oxygen = red , hydrogen = white ) , the m sites 
 are in cyan , and the drude oscillators are green . the oxygen drude 
 displacements are 0.09  , demonstrating the polarization response 
 along the hydrogen bond . both swm4-dp and swm4-ndp were parametrized to reproduce 
 a variety 
 of gas- and condensed - phase properties . unlike fixed - charge water 
 models , which must inherently make compromises between properties 
 in gas and condensed phases , the explicit inclusion of polarizability 
 allows for better response of the polarizable models to changes in 
 environment and therefore better reproduction of thermodynamic and 
 transport properties . the dipole moment of swm4-dp and swm4-ndp serves 
 as an example of this phenomenon . whereas an isolated water molecule 
 reproduces the experimental dipole moment of 1.85 d , transfer to the 
 liquid state induces an increase in the dipole moment to a value of 
 2.46 d , in good agreement with experimental results for water in a 
 hydrogen - bonded environment . the swm4-ndp 
 dipole moment is also sensitive to the presence of charged groups 
 ( figure 5 ) . in contrast , 
 the commonly used tip3p water model has 
 a fixed dipole moment of 2.347 d , irrespective of environment . other 
 properties , such as the change in internal energy upon liquefaction , 
 static dielectric constant , and molecular volume , are also more accurately 
 reproduced using the swm4-dp and swm4-ndp models relative to tip3p . importantly , the polarizable models also reproduce 
 the self - diffusion constant and viscosity of water much more accurately 
 than fixed - charge models , indicating that the time scales of biomolecular 
 phenomena may be better represented using polarizable models . we note 
 that the water models yield self - diffusion constants that are slightly 
 overestimated when the yeh and hummer correction for periodicity is taken into account due to the models being 
 developed prior to publication of that correction . though swm4-dp 
 and swm4-ndp were parametrized using similar methods , and produce 
 similar physical properties , the swm4-ndp model is preferred because 
 of its intuitive connection between the negatively charged drude oscillators 
 and electronic degrees of freedom , and it is the water model that 
 has been used for the development of the remainder of the drude polarizable 
 force field . ( a ) 
 distributions of water dipole moments in bulk water ( || 
 = 2.46 d ) and in a 1 m solution of mgcl2 ( || 
 = 2.53 d ) . dipole 
 moments as a function of distance from ( b ) glutamate-64 and ( c ) lysine-11 
 in ubiquitin . error bars in panels b and c are the root - mean - square 
 fluctuations over the simulation . it is also worth noting a more recent polarizable water model 
 in 
 this series , known as swm6 . like swm4-ndp , swm6 includes a negatively charged drude oscillator and an m site , 
 but also adds two lone pairs out of the h  o  h plane . the addition of these two lone pairs and subsequent reoptimization 
 of the electrostatic parameters of the existing swm4-ndp model were 
 an effort to minimally perturb the properties of swm4-ndp while achieving 
 better agreement with structural properties of water , which would better capture both microscopic 
 and bulk behavior . the swm6 model better reproduces the structure 
 and dipole moment of a water dimer , while preserving or improving 
 upon the dynamic and thermodynamic properties of swm4-ndp . both swm4-ndp 
 and swm6 capture the polarization effects in water and yield better 
 agreement with a wide variety of physical properties than additive 
 water models commonly used in most md simulations , like tip3p , tip4p , 
 and tip5p , with an exception being the 
 temperature - dependence of the density of water . the 
 development of any force field is generally approached in a progressive 
 manner ; that is , parameters are initially optimized by targeting small 
 model compounds ( see examples in figure 6 ) , with those parameters acting as the foundation 
 for the development of parameters for larger entities , such as amino 
 acids , nucleotides , and ultimately proteins and nucleic acids . in 
 this section , the development of drude polarizable models for organic 
 compounds that are the basis for the protein , nucleic acid , lipid , 
 and carbohydrate aspects of the force field will be discussed in detail , 
 with specific emphasis on target data used for optimization of the 
 nonbonded parameters and the role of explicit polarization in the 
 accuracy of the physical model . in general , target data include qm 
 ( gas - phase ) properties , such as molecular dipole and quadrupole moments , 
 water interactions ( energies and distances ) , and dihedral energy scans , 
 as well as condensed - phase properties with experimental target data , 
 such as diffusion coefficients , heats of vaporization ( hvap ) for liquids , heats of sublimation ( hsub ) for crystals , molecular volumes ( vm ) for liquids and crystals , free energies of 
 hydration ( ghydr ) , and static dielectric 
 constants (  ) of neat liquids . the dielectric constant is a 
 crucial consideration when evaluating parameter sets , as it largely 
 dictates the extent to which gas - phase polarizabilities are scaled , 
 as mentioned above . in practice , the scaling factor becomes a freely 
 tunable parameter during optimization , as will be discussed below . examples 
 of model compounds used in the development of the drude 
 polarizable force field : ( a ) ethanol , ( b ) n - methylacetamide , 
 and ( c ) benzene . atoms are colored by element ( o = red , h = white , 
 c = gray , n = blue ) . lone pairs attached to atoms 
 treated with anisotropic polarization in panels a and b are cyan . drude oscillators in panels a and b are shown in green . in the case 
 of benzene in panel c , a perturbing charge of + 0.5 e ( dark blue sphere ) is placed below the center - of - mass of the ring 
 to illustrate the response of the drude oscillators ( colored green 
 in the unperturbed system and yellow when perturbed ) . the drude oscillators 
 move toward the center of the benzene ring , and also out - of - plane , 
 toward the perturbing charge . aliphatic groups are 
 a major component of biological macromolecules , including hydrophobic 
 lipid tails , amino acid side chains , and the backbone of carbohydrates . alkanes are a fundamental component of any force field , but additive 
 models suffer from an inability to respond to a high - frequency oscillating 
 field that gives rise to the optical dielectric constant , inf , according to the clausius  mossotti equation . as a consequence , additive models lead to a 
 systematic underestimation of alkane dielectric constants , as low 
 as 1 . while the dielectric constants of alkanes are very low , this 
 underestimation can have important consequences for free energies 
 of solvation , which scale with ( 1  1/ ) according to 
 the born approximation . given that electronic 
 polarization is central to the dielectric response in low - polarity 
 media like alkanes , its inclusion in the functional form is likely 
 to improve these physical properties . for this reason , the next logical 
 step toward a fully functional polarizable force field was the development 
 of a polarizable model for alkanes . these alkane parameters were optimized according to the method outlined 
 in section 3.3.1 , with two minor modifications . during parameter fitting , charges 
 tended to accumulate asymmetrically on carbon atoms . while this outcome 
 is not without physical basis , it represents a conformation - dependent 
 solution to the fitting algorithm and reduces transferability between 
 compounds . given that computing esp maps for a large number of rotamers 
 of long alkanes was not computationally tractable , a restraint was 
 imposed on the fitting algorithm such that the charge of any carbon 
 atom ( qc ) was restricted to qc = xqh , where x is the number of hydrogen atoms ( each carrying a charge qh ) bonded to the carbon in the chx group , which by convention is electroneutral . final 
 parameters for chx groups were obtained 
 by averaging over several alkanes in multiple rotameric states . the 
 second modification of the fitting protocol concerns the scaling of 
 atomic polarizability values . the general protocol described above 
 prescribes a scaling factor of 0.85 to be applied if the esp maps 
 are calculated at the b3lyp / aug - cc - pvdz level of theory . in the case 
 of alkanes , applying this scaling led to poor agreement between modeled 
 and experimental observables ( i.e. , with the alkanes a scale factor 
 of 1 was applied ) , indicating that a single scaling factor for all 
 classes of molecules is not appropriate ( see table 1 ) . having developed a transferable 
 set of parameters for ethane , propane , 
 butane , isobutene , and pentane , the force field was assessed in terms 
 of thermodynamic , structural , and dynamic properties . specifically , 
 the condensed - phase behavior of these alkanes was evaluated in terms 
 of hvap , vm , static dielectric constants , and ghydr . notably , the polarizable model improved upon hvap values substantially , achieving agreement 
 within 2% of experimental data , whereas the additive models underestimated 
 hvap by up to 18% . the enhanced 
 interaction due to electronic polarizability resulted in a slight 
 underestimation of diffusion constants , but the magnitude of error 
 was no worse than the charmm additive model , and led to better agreement 
 with nmr relaxation times , which arise from rotational correlation 
 times of c  h bond vectors . most importantly , the drude polarizable 
 alkane model produced significantly better agreement with experimental 
 static dielectric constants . whereas the additive model produced nearly 
 uniform values of 1.0 , irrespective of alkyl chain length ( in the 
 series of ethane , propane , butane , isobutene , heptane , and decane ) , 
 the drude polarizable model produced values in excellent agreement 
 with experimental data , ranging from 1.71 to 2.13 . thus , the polarizable 
 alkane model represented a significant improvement over the existing 
 additive model , and explicit polarization was critical in obtaining 
 good agreement with a wide range of condensed - phase behavior , while 
 still retaining good agreement with qm conformational energetics . these results have implications for the accurate 
 modeling of hydrophobic environments , such as lipid bilayers and the 
 interiors of folded proteins ( see section 4.1.2 ) . following the development 
 of the drude alkane model , parameters for aromatic compounds transferability of parameters was evaluated by combining benzene 
 and methyl parameters from the drude alkane force field to create a model for toluene . these compounds 
 form the basis for aromatic amino acids such as phenylalanine and 
 tyrosine . while many of the structural and dynamic features of the 
 benzene and toluene liquids were comparable between the additive and 
 polarizable force fields , the drude model produced more accurate quadrupole 
 moments and dielectric constants . the additive models of benzene and 
 toluene showed similar behavior to the alkanes described above , producing 
 nearly uniform  values of 1.0 , considerably lower than the 
 experimental values of 2.3 ( benzene ) and 2.4 ( toluene ) ; the drude 
 polarizable model ( applying an atomic polarizability scale factor 
 of 0.724 ) yielded better agreement . of 
 particular importance to biomolecular applications of these aromatic 
 parameters is the interaction of water molecules with aromatic groups . 
 given the electronic structure of aromatic rings , the  electron 
 system gives rise to an accumulation of negative charge above and 
 below the ring with partial positive charge accumulating at the edges . the drude polarizable model produced a better - defined 
 population of water molecules in the first peak of the radial distribution 
 function around the ring at 3.5  , indicative of more structured 
 interactions that reflect this hydrogen bonding . subtle differences 
 in hydration shells around benzene were evident , suggestive of fundamental 
 differences in the atomic interactions in additive versus polarizable 
 systems . a proper description of these subtle effects is an important 
 component of the polarizable model , as it produces a more realistic 
 representation of hydration around these biologically important moieties . 
 in this regard , the drude oscillator model is particularly useful , 
 as other methods , like fluctuating charge , will not capture out - of - plane 
 polarization , as charge can only be redistributed among atoms in the 
 planar ring system . more recently , the benzene model was extended 
 to yield improved agreement with respect to qm data for cation 
 interactions by inclusion of a virtual particle in the center of the 
 benzene ring that uses additional lj interactions to optimize the 
 interaction orientation and energy . the 
 development of polarizable force fields for carbohydrates and nucleic 
 acids rely , in part , upon the development of accurate parameters for 
 cyclic ethers , which form the backbone of ribose and deoxyribose rings , 
 as well as the backbones of the wider classes of furanoses and pyranoses . 
 to this end , ethers , 
 while generally of low polarity , have the ability to participate in 
 hydrogen bonding and ion coordination via their oxygen atoms , but 
 their alkyl substituents are nonpolar . thus , a central challenge in 
 developing an accurate force field for ethers is the balance between 
 electrostatic , dispersion , and repulsive forces . the original 
 electrostatic parametrization method of anisimov et al . was expanded to include a total of eight grids 
 around each of the ethers with perturbing charges and an esp map at f = 1.3 , as well as another esp map at f = 2.2 in addition to the five surfaces described above . these additional 
 surfaces correspond to hydrogen - bonding geometries ( hydrogen and heavy 
 atoms , respectively ) and were necessary to properly account for polarization 
 response upon hydrogen - bond formation . parametrization initially focused 
 on two cyclic ethers , tetrahydrofuran and tetrahydropyran , which form 
 the backbone of furanose and hexopyranose ring systems found in carbohydrates . in addition , cyclopentane and cyclohexane parameters were developed 
 to extend the linear alkane series to these cyclic species . following 
 the goal of transferability , parameters for all chx groups not directly bonded to the oxygen of the ether were 
 transferred from the alkane force field , and the parameters developed for cyclic ethers were subsequently 
 transferred to a series of linear molecules . polarizability scaling 
 was applied to all species developed in this series of parameters , 
 including cyclic alkanes . whereas such scaling was not necessary in 
 the development of the linear alkane series , the scaling here reflects 
 the complex electronic properties , likely due to correlation effects , 
 which lead to cyclic species behaving more like polar compounds . a key outcome in the development of the 
 drude polarizable force 
 field for ethers is the simultaneous agreement in relative energies 
 of various conformations and their dipole moments . several rotamers 
 of linear species , as well as different conformers for the cyclic 
 species , were assessed . while the additive model reasonably reproduced 
 the relative conformational energies , the dipole moments of these 
 conformers were generally too large . this outcome reflects 
 the ability of the polarizable model to accurately model the complex 
 electronic structure of ring systems as well as the response of linear 
 species to rotation about various covalent bonds . the implications 
 for biomolecular force fields are clear ; the conformational sampling 
 of furanose and hexopyranose sugars is important to the stability 
 and dynamics of nucleic acids and glycosylated proteins . other condensed - phase properties were reasonably reproduced by the 
 additive model , indicating that the addition of polarizability in 
 the context of ethers did not necessarily improve certain aspects 
 of the physical model . this limitation 
 is likely due to the desire for transferability of parameters between 
 molecules ; the cyclic molecules represented a distinct improvement 
 as a consequence of including explicit polarization , but linear species 
 did not improve . as a result , the ether force field was subsequently 
 revised in 2010 , using a parametrization 
 scheme that allowed for independently tuning lj parameters on linear 
 and cyclic ether oxygen atoms , as well as carbon atoms bonded to ether 
 oxygen atoms . previous parametrization efforts assumed that the value of a in eq 7 should 
 be held constant at a value of 2.6 . in the refinement of the ethers , the parametrization protocol tuned the value of individual atomic 
 thole screening factors , an approach that was introduced on the basis 
 of the need to improve the molecular polarizability tensor of the 
 amides , as discussed below . the application of a 0.7 scaling factor 
 for ether gas - phase polarizabilities led to a systematic underestimation 
 of liquid dielectric constants ; the updated ether parameter set applied 
 a scaling factor of 0.85 , which achieved significantly better agreement . this updated ether parameter set , with rescaled 
 atomic polarizabilities and atom - specific thole factors , forms the 
 basis of the carbohydrate force field described below . given the ubiquity of 
 alcohol groups in proteins ( amino acids like serine , threonine , and 
 tyrosine ) , nucleic acids ( 2- and 3-hydroxyl groups ) , 
 carbohydrates , and lipids , these functional groups are a major component 
 of any force field . properties of alcohols are interesting not only 
 due to this prevalence , but because their hydration is dictated by 
 competition between polar and nonpolar interactions . thus , a proper 
 description of alcohol interactions with various media will certainly 
 require electronic polarization . the first application study 
 of a polarizable alcohol using the drude oscillator model was conducted 
 by noskov et al . results from that study 
 produced excellent agreement with solution properties ( such as dielectric 
 constant and diffusion of water and ethanol ) and provided important 
 insight into ethanol hydration properties . supporting the hypothesis 
 that water structure is affected by ethanol , the authors found that 
 water water hydrogen bonding was depleted in the vicinity of 
 the ethanol hydroxyl group but that the dominant effect in water perturbation 
 was due to restructuring of the water molecules in the second hydration 
 shell around ethanol . , a 
 more generalized parameter set for polarizable alcohols was subsequently 
 produced to encompass a larger series of primary and secondary alcohols 
 two years later , including a more refined 
 model of ethanol that included lone pairs on the hydroxyl oxygen atom 
 ( figure 5a ) . parametrization 
 of electrostatic terms was carried out in the same manner as for the 
 ethers described above . the parameters that were developed were further 
 evaluated against qm water interaction energies and geometries . the 
 inclusion of explicit polarization led to significantly better agreement 
 with qm water interactions . in the final parametrization of the model , 
 atomic polarizabilities were selectively scaled ; to preserve the quality 
 of the hydroxyl water interactions ( as well as other properties 
 discussed below ) , the polarizability of the hydroxyl oxygen atom was 
 not scaled , but  values of carbon atoms in the alcohols were 
 scaled by 0.7 . off - diagonal terms ( pair - specific lj parameters that 
 override normal lorentz  berthelot combination rules , called 
 nbfix in charmm ) were introduced for the interaction of primary and 
 secondary alcohol oxygen atoms with water oxygen atoms . while this 
 extension of the normal combination rules may limit transferability 
 between molecules , it was necessary to simultaneously achieve good 
 agreement with experimental properties and qm target data . the final 
 model was evaluated in training set molecules as well as longer - chain 
 alcohols that were not part of the training set ( 2-butanol , 1-propanol , 
 1-butanol ) . the polarizable model was found to be a significant improvement 
 over the additive charmm force field for all quantities . most notably , 
 the polarizable alcohol series improved hvap and ghydr , among quantities , 
 and reduced the average error in static dielectric constants from 
 35.9% in the additive model to 2.3% in the drude model . 
 given the similarity in internal terms between these two force fields , 
 the polarization response is clearly responsible for the improvement 
 in agreement of these values , an outcome that reflects the improved 
 treatment of nonbonded terms . while the additive force field was intrinsically 
 overpolarized to mimic condensed - phase conditions , this approximation 
 is ultimately inadequate , and the drude model showed another critical 
 feature , the response of the molecular dipole moment to environment . 
 whereas the additive charmm force field for alcohols produced dipole 
 moments that were largely invariant between gas , neat alcohol , and 
 aqueous systems , the drude model responded dramatically to these different 
 environments , shifting from low dipole moments in the gas phase ( values 
 in agreement with qm calculations ) and , when solvated by benzene , 
 to much higher values in aqueous solution , in agreement with previous 
 theoretical calculations . as a final 
 point , the drude alcohol parameter set influenced the dipole moment 
 of water molecules hydrating the alcohols . small variations in water 
 dipole moment were noted as a function of distance around ethanol 
 as well as hydrogen - bonding geometry ( with water serving as either 
 a donor or acceptor ) , suggesting that mutual polarization between 
 ethanol and water is an important factor in dictating the interactions 
 in aqueous systems . the ability to model these subtle effects in the 
 polarizable alcohol series reflects a distinct advantage over additive 
 force fields , in which rigid water models and nonpolarizable alcohols 
 are incapable of responding to such effects , which may have implications 
 for the dynamics of the many biomolecules that contain hydroxyl groups . the drude model of methanol was recently examined in the context 
 of kirkwood  buff integrals of methanol in aqueous solution using the swm4-ndp model of water . by assessing solutions of varying methanol concentration 
 ( from pure methanol to pure water ) , buff theory 
 relates the microscopic details of the solution structure with thermodynamic 
 ( macroscopic ) properties . while the drude model of methanol resulted 
 in better density in the neat liquid and increasing water self - association 
 in aqueous methanol solution , ultimately there was little difference 
 between the additive charmm and drude polarizable models of methanol 
 with respect to experimental kirkwood  buff integrals , excess 
 coordination numbers , activity derivatives , and excess molar gibbs 
 energy . this outcome suggests that additional improvements can be 
 made to the drude methanol model ( and perhaps alcohols in general ) 
 to even better reproduce the balance of forces between alcohols and 
 water . future work will also assess the quality of alcohol interactions 
 with the swm6 water model , which , like 
 the alcohols , includes lone pairs on the oxygen atom and consequently 
 may achieve a better balance of intermolecular interactions . given that protein backbones 
 are polyamides and two amino acid ( asparagine and glutamine ) side 
 chains contain amides , an accurate model for amide groups is an essential 
 part of a polarizable force field for proteins . the balance of forces 
 between amides and water is an important determinant of protein conformational 
 dynamics . to this end the parametrization was carried out as described above , 
 and it was during the optimization of the amides that the use of atom - based 
 thole screening factor was introduced in the parameter optimization , 
 as was used in the refinement of the ether parameters . using atom - specific thole screening factors 
 is justified due to the electronic resonance in chemical groups such 
 as amides , for which electronic response along different axes may 
 be different , and inclusion of this term leads to significant improvements 
 in the treatment of molecular polarizability tensors . notably , the 
 final , optimized amide parameter set directly used the gas - phase atomic 
 polarizabilities to reproduce the high dielectric of neat n - methylacetamide ( nma , 100 at 373 k ) . additive models 
 of nma ( figure 6b ) 
 with the charmm force field and formamide with opls - aa had previously 
 been found to underestimate neat liquid  values by 70% and 
 50% , respectively , indicating that the mean - field 
 approximation inherent in additive models is an inadequate representation 
 of these species . the drude model for nma reproduced a wide range 
 of gas - phase qm and condensed - phase experimental data and was the 
 first force field to yield dielectric values in near - quantitative 
 agreement with experimental values , including temperature dependence . cooperativity among dipole moments ( head - to - tail alignment of dipoles 
 such that their magnitudes are enhanced along the dipole vector ) between 
 hydrogen - bonded nma molecules explains the physical behavior , a property 
 that required accurate treatment of the molecular polarizability tensor 
 by using atom - based thole factors , as described above . in the drude 
 model , the kirkwood gk factor was considerably 
 larger ( 4.6  0.6 ) than the charmm additive model ( 3.0  
 0.2 ) , indicating a greater degree of cooperative dipole alignment . the decrease in dielectric as a function of temperature is rationalized 
 in the same way . individual nma molecules experience a minimal ( 0.1 
 d ) decrease in dipole moment , but the value of  drops considerably . 
 while the additive model experiences a small decrease in dielectric , 
 its values are well below those observed experimentally , but the drude 
 model reproduces experimental behavior nearly quantitatively . the elevated temperature 
 causes frequent breaking of hydrogen bonds , thus disrupting cooperativity 
 among dipole moments and ultimately a decreased value of dielectric . as a final note on dielectric values , the drude polarizable model 
 correctly reproduced the paradoxical trend among methylated amides . 
 rather than increasing as a function of hydrogen - bonding capacity , 
 the  values at 373 k for acetamide , nma , and n , n - dimethylacetamide are 66.3 , 100 , and 26 , respectively . the  values of nma and n , n - dimethylacetamide can be rationalized intuitively by differences 
 in hydrogen - bonding capacity , but the value of acetamide is somewhat 
 surprising , as it can participate in three hydrogen bonds per molecule , 
 while nma can participate in two . however , simulations with the drude 
 model revealed that hydrogen bonding to the cis amide 
 proton results in anticorrelated dipole moments , opposing the dipole 
 enhancement via the trans proton hydrogen bonding , 
 ultimately leading to an intermediate dielectric value . this work 
 on polarizable amides represented an 
 important advancement toward the construction of a polarizable protein 
 force field , as nma is a model of the polypeptide backbone . the ability 
 of the drude nma parameters to reproduce qm water interactions , ghydr , and dipole cooperativity has important 
 implications for simulations of proteins , as will be discussed later . building on this work , which involved primarily simulations of 
 neat acetamide liquid simulations , the properties of the acetamide 
 series in aqueous solution using the drude polarizable force field 
 have recently been assessed in greater detail . this work resulted in updated nma parameters based on an 
 initial reparametrization of internal terms to reflect polypeptide 
 geometric data taken from crystal surveys , with subsequent reoptimization 
 of nonbonded ( electrostatic and lj ) terms . this newer model of nma 
 serves as the basis for the drude polarizable protein force field 
 that will be described below . unlike the previous model of nma that 
 reproduced neat liquid properties well the resulting model reproduced many 
 experimental quantities over a range of nma and acetamide concentrations 
 in water , including kirkwood  buff integrals , water coordination numbers , solution density , partial molar 
 volumes , and molar enthalpy of mixing . these outcomes indicate that 
 the newer drude polarizable model for amides accurately balances solute  solvent , 
 solute  solute , and solvent  solvent interactions . such 
 balance is a crucial component of an accurate biomolecular force field , 
 since protein conformational equilibria depend , in part , on the balance 
 of intramolecular and protein  solvent interactions . several biologically important compounds 
 are comprised of nitrogen - containing 
 heterocyclic aromatic moieties , such as amino acids histidine and 
 tryptophan , all of the bases in dna and rna , and several enzyme cofactors , 
 and they occur in a wide range of medicinal compounds . accordingly , 
 a series of heterocyclic aromatic compounds ( pyrrole , imidazole , pyridine , 
 pyrimidine , indole , and purine ) were parametrized to create the building 
 blocks for these biologically relevant molecules . the heteroaromatic series was developed much in the same 
 manner as benzene and toluene described above and targeted many of 
 the same experimental data , including neat liquid and crystal simulations . 
 for heteroaromatics , the gas - phase polarizabilities were scaled by 
 0.85 to achieve good agreement for pyridine and pyrrole dielectric 
 constants . whereas both scaled and unscaled polarizabilities yielded 
 good agreement for pyridine , the unscaled parameters for pyrrole overestimated 
  , necessitating the scaling . achieving agreement with the liquid 
 dielectric constant is the guiding principle for applying polarizability 
 scaling , and the scaling factor can be an additional variable in parameter 
 fitting . prior to the development of the nitrogen - containing 
 heterocycles , lone pairs were only applied to carbonyl oxygen atoms 
 that were treated anisotropically . in 
 the case of the heterocycles , virtual sites were added to the nitrogen 
 atoms to represent the in - plane lone pairs of these atoms . the addition 
 of these virtual sites improves agreement with the qm esp maps and 
 is chemically intuitive , given the nature of these compounds . additional 
 details of the nitrogen - containing heterocycles will be discussed 
 below in the context of nucleic acid parametrization . thiols 
 and sulfides feature prominently in many chemical and biological systems , 
 particularly in proteins ( methionine and cysteine ) as well as in medicinal 
 chemistry applications . the sulfur atom is polar and strongly polarizable ; 
 thus , an accurate representation of these chemical groups containing 
 sulfurs is an important component of any force field . parameters for 
 sulfur - containing compounds using the drude oscillator model were 
 derived by zhu and mackerell in 2010 . similar to the polarizable force field for alcohols , the polarization of sulfur atoms was treated 
 anisotropically , yielding good agreement with orientation - dependent 
 qm water and ion interaction energies and a considerable improvement 
 over the additive charmm22 parameters for these same compounds . given 
 the sensitivity of sulfur to local electric fields , polarizability 
 scaling is a particularly important consideration . it was found that 
 a scaling factor of 0.7 for thiols led to good agreement with experimental 
 dielectric constants , while 0.85 was applied to dimethyl disulfide . 
 upon further investigation into ethyl methyl sulfide , it was determined 
 that a scaling factor of 0.6 was more suitable for reproducing its 
 condensed - phase properties and gas - phase dipole moment , and these 
 are the parameters that are used in biological moieties such as those 
 found in methionine . as in the case of alcohols , off - diagonal lj terms 
 ( nbfix ) between sulfur and water oxygen atoms were necessary to simultaneously 
 achieve agreement between the properties of the neat thiol and disulfide 
 liquids and ghydr . the assembled body of model compounds that has 
 been derived serves as the basis of the biomolecular drude polarizable 
 force field that will be described in the following sections . the 
 judicious choice of target data , from both qm calculations and condensed - phase 
 experimental data , and iterative refinement protocols that have been 
 used for many years in the development of the additive charmm force 
 field have shown that parameters are transferable across molecules . 
 freedom to tune thole screening factors on a per - atom basis , the use 
 of atom - pair - specific lj parameters , and scaling of gas - phase polarizabilities 
 based on the nature of the model compounds allow for robust models 
 of a variety of molecules . the inclusion of electronic polarizability 
 provides unique insight into microscopic properties at the small molecule 
 level , a feature that extends to biological macromolecules , as will 
 be shown below . monatomic ions 
 are prevalent throughout the aqueous cellular milieu and play important 
 roles in neuronal excitation , folding and structural stability of 
 proteins and nucleic acids , and catalysis . ions can be strongly polarizing , a property that can dramatically 
 influence their solvation and interactions with biomolecules but one that is impossible to model in additive 
 force fields . the first set of ion parameters developed as part of 
 the drude-2013 force field was a series of monovalent alkali and halide 
 ions , parametrized using the swm4-dp water model . this ion parameter 
 set became known as ah / swm4-dp (  ah  
 for  alkali and halides  ) . as monovalent ions have very 
 few tunable parameters in the polarizable model ( only atomic polarizability 
 and lj terms ) , the parametrization strategy employed in that work 
 was to take gas - phase atomic polarizabilities ( no scaling for cations , 
 but scaled by a factor of 0.7 for anions as described in section 3.2 ) and optimize 
 lj parameters to obtain accurate ghydr . given that the biological function of ions also relies on correct 
 coordination geometry in solution and when interacting with biomolecules , 
 the optimized parameters were validated against structural properties 
 of ion hydrates and bulk solution using available ab initio and experimental 
 data . the ion parameters produce nearly perfect results for neutral 
 salt ghydr ( which are more reliable 
 than the ghydr values for individual 
 ions ) , and the structures of the hydrated alkali ions ( which have 
 highly ordered water configurations ) agreed well with ab initio geometries 
 and energies . the structural properties of halide hydrates , which 
 have perturbed structures , showed some disagreement with ab initio 
 calculations , but much of the disparity could be attributed to limitations 
 in the water model , swm4-dp . further , both structural and thermodynamic 
 halide properties in bulk solvent agreed well with experimental data , 
 indicating that the ion parameters themselves were a reasonable representation 
 of physical reality . the same general parametrization strategy 
 was employed in the newest version of drude polarizable ions , which 
 include monovalent and divalent ions . the updated parameters were generated using the preferred swm4-ndp 
 water model . while the target data were the same as in the previous 
 parametrization , and comparable accuracy notably , divalent 
 ions ( ca , mg , zn , and sr ) were parametrized for the first time . divalent ions are 
 prone to  polarization catastrophe  due to a coulombic 
 singularity that can arise when the drude particle on the ion ( which 
 carries no lj terms ) interacts with the core charge on the water oxygen . this problem was circumvented by using a through - space thole screening 
 of the same form as that used for neighboring ( bonded ) dipoles described 
 above for selected pairs ( eq 7 ) . this algorithm , called nbthole ( for  nonbonded thole  ) 
 in the charmm parameter file , was an optimal solution to the overpolarization 
 problem , as additional restoring forces and/or assignment of lj terms 
 to drude oscillators of the ions led to stable simulations , but it 
 was not possible to achieve sufficient agreement with monohydrate 
 energies . the use of the nbthole formalism allowed for both stable 
 bulk simulations and good agreement with ab initio target data . it 
 should be noted that the agreement of monovalent ions in this parameter 
 set is somewhat better than that of the divalent ions , due to charge - transfer 
 effects that can not be modeled by the drude polarizable functional 
 form . despite this limitation , the divalent ion series still reproduces 
 experimental and ab initio target data well , including the hydration 
 structure in bulk solution . ion 
 interactions , especially in mg binding to rna , which 
 is known to facilitate folding and stabilize complex tertiary structures . finally , the diffusion constants of the monovalent and divalent 
 ion series were assessed and found to be in good agreement with experimental 
 values . the ion parameters thus represent a good physical model for 
 these species in terms of thermodynamic , transport , and gas - phase 
 hydrate properties . the existing drude polarizable ion series 
 has been parametrized 
 at infinite dilution ; thus , the behavior of the parameters in concentrated 
 ion solutions requires additional exploration . at higher concentration , 
 significant ion  ion interactions can affect dynamics and bulk 
 solution properties , though it is clear that there are interesting 
 physical insights to be obtained , as concentrated salt solutions lead 
 to augmented water dipole moments ( figure 5a ) . one approach for assessing the quality 
 of ion parameters in solution is the calculation of osmotic pressure , 
 for which a straightforward computational approach has been developed . this method has subsequently been applied in 
 an initial examination of concentrated solutions of drude polarizable 
 ions and also in the context of balancing 
 interactions among ions , water , and nucleic acid functional groups . these osmotic pressure calculations represent 
 an important advancement in force field parametrization and validation , 
 as the target data are well - defined and are good indicators of the 
 balance between ion  ion , ion  neutral solute , and ion  water 
 interactions , critical components of an accurate force field . building 
 upon the small molecule parametrization efforts described above , the 
 first protein force field based on the drude oscillator was completed 
 in 2013 . emphasis was placed on transferability 
 between compounds during the development of the drude-2013 force field 
 for small molecules , with the end goal being the ability to combine 
 these parameter sets into a full model for biomolecules . the initial 
 protein model ( called  drude-1  ) was produced using 
 the simplest means , by combining nma and ethane 
 into a model of the alanine dipeptide , which is a central model in 
 the force field parametrization of proteins . additional parameter 
 refinement is necessary when creating new connectivity between discrete 
 molecules ; to this end , internal terms were optimized , and 2-dimensional 
 cross - term energy correction maps ( cmap ) were generated on the basis of high - level gas - phase qm data . the 
 simple drude-1 model did not produce good agreement with j - coupling data for the model ( ala)5 peptide , as extended 
 ( c5 , figure 7 ) conformations were overly favored . this behavior arose because 
 the two neighboring peptide bonds are oriented in opposite directions , 
 such that their local dipole moments cooperatively enhanced each other , 
 thereby maximizing hydrogen bonding with water . these cooperative 
 dipole  dipole interactions also led to peptide conformations 
 that under - stabilized helical conformations , in favor of elongated 
 conformations . analysis of dipole moments of ( ala)5 indicated 
 that this behavior was due to improper treatment of the polarization 
 response associated with direct application of electrostatic parameters 
 based on nma to the polypeptide . as a result , further refinement of 
 this simple model was necessary . to address the proper polarization 
 response associated with peptide bond peptide bond cooperative 
 dipole interactions as a function of  and  values , a 
 model (  drude-2  ) was developed by individually targeting 
 esp maps of the alanine dipeptide in five different conformations 
 ( r , l , c5 , ppii , and 
 c7eq , figure 7 ) and taking the average of these parameter sets . the rationale for 
 this approach was to capture the changes in esp as a result of dipole 
 reorientation concomitant with changes in the relative orientation 
 of the two peptide bonds . while the drude-2 model yielded better conformational 
 properties relative to drude-1 , it suffered from defects in representing 
 gas - phase qm dipole moments , and the conformational properties of 
 a model ( ala)5 polypeptide in aqueous solution were inadequate ; 
 extended conformations were still overpopulated despite a reduction 
 in dipole enhancement relative to drude-1 . the r , l , c7eq , and ppii conformations are aligned along the cc 
 bond vector to illustrate the relative positioning of the backbone 
 as a function of  and . the final model , called  drude-3  by lopes 
 et al . 
 and ultimately  drude-2013 ( thus the basis of the name 
 of the force field described in this review ) , was produced by refining 
 the drude-2 model with a monte carlo / simulated annealing ( mc / sa ) protocol 
 that targeted a wide variety of target data , including the polarizability 
 and dipole moments of the alanine dipeptide , relative energetics of 
 the ( ala)5 peptide , and qm interaction energies between 
 water and the alanine dipeptide with different backbone conformations . this strategy was motivated to obtain an improved balance among intramolecular 
 electrostatic interactions ( from gas - phase target data ) and condensed 
 phase properties ( water interactions ) . the mc / sa protocol generates 
 random changes in an n - dimensional vector ( where n is the number of parameters being fit ) to move through 
 parameter space by minimizing an error function calculated as the 
 root - mean - square difference between all qm target data and mm values . 
 during mc / sa , the temperature is an important component of the 
 algorithm , as increases in the error function are accepted according 
 to the metropolis criterion to promote exhaustive sampling of parameter 
 space . targeting this wide range of data simultaneously , while also 
 including the longer ( ala)5 , an extension of the typical 
 additive charmm parametrization protocol using only the alanine dipeptide , 
 was a critical component in creating a viable polarizable model for 
 the protein backbone . the longer ( ala)5 peptide was utilized 
 to account for the conformational properties that result from communication 
 between neighboring peptide - bond dipoles , including i to i + 4 interactions occurring in -helices , 
 a property that was poorly represented by directly using nma . the drude-2013 protein force field was able to successfully model 
 relative energies between different peptide backbone conformations , 
 while also suitably reproducing water interaction energies . the balance 
 between intrapeptide and solvent interactions plays an important role 
 in determining the conformational ensemble of a protein in folded 
 and unfolded states , and the ability of the drude-2013 force field 
 to capture these effects is an important aspect of its parametrization . 
 to achieve this agreement , gas - phase polarizabilities were scaled 
 by a factor of 0.85 , similar to the scaling applied to small molecules 
 described above . a few additional details of the drude-2013 
 protein force field 
 merit discussion here . by convention , the peptide bond group ( including 
 the c and h atoms ) has a net charge of zero . side - chain 
 moieties ( from c onward ) also have a unit charge ( 1 , 
 0 , or 1 depending on the amino acid ) . this convention facilitated 
 combining the small molecule electrostatic parameters into the polypeptide 
 and , notably , allows for straightforward assessment of dipole moments 
 for these groups . these dipole moments have important implications 
 for the properties of the drude-2013 force field and the physical 
 insights it provides . compared to the most recent additive charmm 
 force field , charmm36 , the drude-2013 
 force field has larger backbone peptide - bond dipole moments . despite 
 the inherent overestimation of dipole moments in the additive force 
 field , the magnitude of the dipoles in the drude-2013 force field 
 is in better agreement with qm data . these dipoles also respond to peptide conformation , including longer - range 
 intramolecular interactions , such as the i to i + 4 hydrogen bonds in -helices ( see sections 4.1.1 and 4.1.2 ) . additionally , it was shown that the side - chain 
 dipole moments of tryptophan residues respond very strongly to conformation 
 ( by rotation about 1 torsions and changes in the 
 local environment ) , varying over a significantly wider range than 
 charmm36 . this observation indicates that the polarizable model exhibits 
 sensitivity to variations in the local electric field in response 
 to small conformational changes , a feature that is essentially absent 
 in additive force fields . similarly , the dipole moments of water molecules 
 around charged amino acids vary as a function of distance from these 
 groups , as shown in figure 5b , c . while the nearest water molecules have slightly decreased 
 dipole moments , at the maximum of the radial distribution function 
 around the charged groups , the water dipole moments increase , before 
 converging to the bulk value ( 2.46 d for swm4-ndp ) at longer 
 distance . the mutual polarization response between the charged residues 
 of the protein and solvating waters shows a small but nontrivial response 
 to the properties of microenvironments on the protein surface . the drude-2013 protein force field exhibits greater flexibility 
 than charmm36 , such that the rmsd values of a wide range of proteins 
 are typically larger than in the additive simulations . this outcome 
 suggests that the polarizable model allows for sampling of a broader 
 conformational ensemble , while still remaining in reasonable agreement 
 with available nmr data . the drude-2013 model represents the first 
 generation in protein force fields with the drude oscillator model , 
 and future refinements of the model will target better agreement with 
 solution nmr and other data to even more accurately represent conformational 
 dynamics of proteins . recently , the interactions between the 
 drude-2013 protein force 
 field and several biologically important ions ( na , k , ca , and cl ) have been refined . off - diagonal atom - pair - specific lj terms ( nbfix ) 
 and through - space thole screening ( nbthole ) parameters were derived 
 by targeting qm geometries and energies of ion  small molecule 
 complexes . these small molecules included nma ( as a representation 
 of the protein backbone , and side - chain amides of glutamine and asparagine ) , 
 ethanol ( a model of serine and threonine ) , and propionic acid ( to 
 model aspartic and glutamic acids ) . the newly parametrized interactions 
 were evaluated in terms of solvation free energies ( gsolv ) in water , nma , and ethanol ; both single 
 ion and salt gsolv values were 
 in near - quantitative agreement with experimental data . even more importantly , 
 the pair - specific interactions were able to reproduce transfer free 
 energies between solvents , indicating that the polarizable ion models 
 can respond to changes in dielectric of the surrounding solvent , a 
 key property of the polarizable model . the drude polarizable model 
 also manifested cooperativity in geometry changes as a function of 
 ligands . for instance , the oxygen  ion distances in nma complexes 
 increased as a function of the number of nma molecules in the complex 
 ( from one to four ) . this behavior reflects the outcome of the qm calculations , 
 suggesting that the empirical model can recover the cooperative dipole 
 behavior on the qm level . the final , and most significant , outcome 
 of this parameter refinement was shown in the behavior of protein  ion 
 complexes . by assessing 30 enzyme structures that have defined ion 
 binding sites , it was found that the polarizable model correctly reproduced 
 structural and thermodynamic data . subsets of each ion binding site 
 were extracted from structures following short md simulations for 
 qm and mm evaluation . whereas the additive charmm36 force field systematically 
 overestimated binding affinity , from 37 kcal mol in the case of k to as much as 80 kcal mol in the case of ca , the drude polarizable model was 
 in near - quantitative agreement for monovalent ions and within 10% 
 of qm values even for ca , which strongly polarizes its 
 coordinating ligands . this investigation also allowed for the first 
 quantitative assessment of many - body effects in protein the drude-2013 force field was compared with charmm36 and 
 the amoeba polarizable force field and was found 
 to best reproduce the energy associated with nonadditivity . amoeba 
 overestimated this phenomenon , and by definition the nonadditive contribution 
 in charmm36 was zero . these effects of polarization may account for 
 as much as 30% of the binding energy ( in the case of ca ) in protein ion complexes , highlighting how traditional fixed - charge 
 force fields can not account for ion binding in strongly polarizable 
 microenvironments such as those found in proteins ; even the most carefully 
 parametrized models are unlikely to ever correctly account for such 
 phenomena . in this regard , the development of parameters that accurately 
 describe such interactions in the drude polarizable model represents 
 a breakthrough in the pursuit of more realistic simulations of protein  ion 
 complexes . a rudimentary 
 nucleic acid 
 force field was proposed in the early work by anisimov et al . when 
 describing the method for deriving electrostatic parameters for the 
 drude polarizable force field . a short 
 ( 5 ns ) simulation of a small dna duplex was carried out using a topology 
 produced by combining parameters of nucleic acid bases derived in 
 that work with sugar and dimethyl phosphate ( dmp ) parameters . the 
 dna duplex remained intact over this short time period , suggesting 
 that the proposed parametrization approach would be suitable for complex 
 ( and highly charged ) systems such as dna , but the progressive increase 
 in rmsd of the structure indicated that further refinement would be 
 needed . building upon the work of lopes et al . in deriving parameters 
 for nitrogen - containing heterocycles , parameters for the five common nucleic acid bases ( adenine , guanine , 
 cytosine , thymine , and uracil ) were produced in 2011 . these nucleic acid base parameters became the 
 basis of the full dna force field , which was published in 2014 . development of the drude-2013 dna force field simultaneously 
 targeted the intrinsic conformational energetics of several model 
 compounds , as well as the condensed - phase behavior of dna . the performance 
 of the force field was assessed in terms of the energetics of these 
 model compounds , crystal survey data on a- and b - form dna , and nmr 
 order parameters , which provide crucial insight into the structure 
 and dynamics of dna . the most important of these local conformational 
 changes include sugar repuckering between north ( n ) and south ( s ) 
 conformations and the bi / bii equilibrium in b - dna that is a function 
 of backbone  and  torsions ( bi , the bi / bii 
 equilibrium modulates the major and minor groove widths in dna , so 
 proper balance of these states is important in establishing a physically 
 realistic model of dna , especially with the goal of a comprehensive 
 force field capable of simulating heterogeneous systems , such as protein  dna 
 complexes . as with the drude-2013 protein force field , combining 
 small - molecule 
 building blocks required additional optimization of bonded terms to 
 yield a stable force field , and suitable model compounds had to be 
 identified that could account , at least in part , for the correlation 
 that is inherent in dna backbone torsions , sugar repuckering , and 
 base orientation . several noncanonical geometries of sugar puckering 
 and base orientation relative to the sugar were assessed to act as 
 target data in the development of the drude model ; such species are 
 important in capturing the energetics during conformational transitions 
 between a- and b - form dna helices . in addition to model compounds 
 shown in figure 8 , 
 the correlation in dna dynamics had to be assessed by conducting simulations 
 of dna in solution . model compounds used in the parametrization of the drude-2013 
 dna 
 force field . for 2-deoxyribose ,  r indicates 
 one of the four standard nucleic acid bases ( adenine , guanine , cytosine , 
 and thymine ) . backbone dihedral angles (  ,  ,  , 
  ,  ) are shown for t3ps and t3pm , and the torsions defining 
 sugar puckering ( 0 , 1 , 2 , 3 , 4 ) and base orientation 
 relative to the sugar (  ) are labeled on 2-deoxyribose . the drude-2013 dna force field 
 better reproduces the two - dimensional 
 qm potential energy surface for the rotations of the  torsion 
 as a function of  than the additive charmm36 model . this agreement 
 was achieved through specific tuning of torsional terms on the small 
 molecule level using model t3pm in figure 8(182 ) and also empirical 
 optimization of electrostatic and lj terms in the phosphodiester backbone . while charmm36 has a pronounced minimum in 
 the bii state and a low barrier to the bi state , the drude model was 
 better able to capture the subtle features of the qm energy surface , 
 including better relative energetics between the bi and bii states . similarly , the drude model reproduced the conformational energies 
 of noncanonical nucleoside geometries , while the additive model generally 
 overestimated the energies . these outcomes suggest that the inclusion 
 of explicit polarization allows for a better representation of subtle 
 conformational change that is critical to the proper representation 
 of dna dynamics . as a result , the drude-2013 force field was able 
 to reasonably reproduce experimentally observed , sequence - specific 
 bii content and s order parameters from nmr for the ecori dodecamer . sugar repuckering was shown to occur on 
 a time scale of 200 ps , in agreement with experimental data . taken together , these data indicate that the 
 polarizable model more accurately describes the kinetics of local 
 conformational change , an attribute that derives from the greater 
 ability of the drude model to reproduce the qm conformational energies , 
 though further testing is required . another important test of 
 any dna force field is the response of 
 the dna to its environment . in solutions with low water activity nucleotide content can modulate 
 this effect , with gc - rich regions stabilizing the a - form to a greater 
 extent than at - rich sequences . the drude-2013 dna force field was 
 assessed by simulating a gc - rich a - form dna sequence in water and 
 a solution of 75% ethanol . the a - form was stable in ethanol and quickly 
 converted to the b - form in water . the results indicate that the drude 
 model is capable of correctly responding to changes in solvent polarity . however , we note that the drude-2013 dna force field does not yield 
 stable z - dna conformations ; ongoing efforts are being made to address 
 this and other limitations of the model as well as to extend the model 
 to rna . given that dna is a polyanion , with charge density much 
 larger 
 than that of most proteins , the electronic response of nucleic acid 
 bases to hydrating water and ions is of particular importance . over 
 the course of the full - length dna simulations , the dipole moments 
 of the nucleic acid bases were systematically larger and more variable 
 than those calculated from the additive charmm36 simulations , despite 
 the inherent overestimation of gas - phase dipole moments in the parametrization 
 of the additive force field . this outcome suggests that the bases 
 in the drude model are more sensitive to the electronic environment 
 within the dna helix and that the additive assumption of overestimating 
 dipole moments may be an inadequate representation of the condensed - phase 
 environment . given the fact that gas - phase polarizabilities were scaled 
 by 0.85 in the parametrization of the bases , the drude model is unlikely to be overpolarized . the functional 
 implications of nucleic acid base dipoles will be discussed later 
 in this review ( see section 4.2 ) . as was the case for water molecules solvating 
 charged residues 
 in the drude-2013 protein force field , water molecules in the major and minor grooves of dna also exhibited 
 distance - dependent variations in their dipole moments . differences 
 in dipole moments can also be seen as a function of geometry ; that 
 is , waters respond slightly differently in the major and minor grooves . 
 at larger distances ( beyond the first hydration shell ) , the water 
 dipole moments relax toward the bulk value ( 2.46 d ) . in contrast , 
 the dipole moments of the rigid tip3p water molecules in the additive 
 charmm36 simulations were invariant , by definition . the inclusion 
 of explicit polarization suggests that there is interesting behavior 
 in hydrating waters that is unaccounted for in simulations with additive 
 force fields . water 
 interactions is an important factor in determining structure and functional 
 dynamics . as a consequence , the polarizable model represents a significant 
 advancement in our understanding of these fundamental forces . in addition to interactions with water , the interactions between 
 polyanionic dna and mobile ions in aqueous solution are also very 
 important to the structure and dynamics of dna , and nucleic acids 
 in general . to tune the parameters for ion  dna interactions , 
 osmotic pressure calculations were conducted 
 on solutions of na  dmp and na qm interaction energies and geometries 
 of dmp  water complexes and the gsolv of dmp were also targeted in the refinement of the lj 
 parameters of dmp and , thus , the phosphodiester backbone of dna . these 
 attributes are important in characterizing the hydrated structure 
 of the dna backbone and the overall balance of interactions within 
 dna and between dna and the surrounding solvent . another critical 
 refinement in the drude-2013 dna force field was balancing interactions 
 between na ions and nucleic acid bases . although the parameters 
 for the ions and bases were developed using similar strategies , their interplay 
 had not been previously validated . it was observed that the na ions were too strongly attracted to nitrogen atoms in the 
 purine and pyrimidine ring systems . to better match qm interaction 
 energies and distances , off - diagonal nbfix terms had to be applied 
 to preserve the quality of the base parameters and simultaneously 
 produce reasonable interactions with ions . the completed drude-2013 
 dna force field was shown to quantitatively reproduce the extent of 
 charge neutralization by ions around dna , in accordance with counterion 
 condensation ( cc ) theory . while both 
 the charmm36 and drude-2013 force fields yielded counterion shells 
 around the dna , the level of condensation in the drude model was improved , 
 indicating that the nonbonded interactions in the polarizable model 
 provide a better physical description of the balance of interactions 
 among ions , water , and dna . carbohydrates are 
 ubiquitous in biological systems , playing a central role in metabolism 
 as well as serving as covalent adducts for proteins and lipids . the force field for linear polyols was introduced in 2013 , and parameter sets for hexopyranose and aldopentofuranose monosaccharides have recently been completed . the linear 
 polyol series was parametrized by targeting the conformational energetics 
 and molecular dipole moments of nearly 2000 unique conformers of compounds 
 with two to six alcohol groups ( ethylene glycol , glycerol , and a series 
 of tetritols , pentitols , and hexitols ) . the additive charmm parameters 
 for these molecules overestimated energies and dipoles in many of 
 these conformers due to the inherent overpolarization of the alcohol 
 groups . this convention produces considerable repulsion between neighboring 
 oxygen atoms due to the large partial charges . in contrast , the drude 
 model better reproduced the target qm data , demonstrating the role 
 of explicit polarization in recovering gas - phase conformational energetics 
 and dipole response . this behavior is a critical component of the 
 carbohydrate force field , as the balance of intra- and intermolecular 
 hydrogen bonding is crucial to the physicochemical properties of this 
 class of molecules . correctly accounting for the polarization effects 
 allowed the drude-2013 force field to improve upon the additive model 
 in all quantities considered in the validation , including dielectric 
 constants , vm , hvap , ghydr , crystal 
 volumes and associated lattice parameters , and aqueous solutions of 
 varying polyol concentration . as a notable example , the hvap for glycerol with the additive charmm force 
 field was underestimated by 14.2% because the intrinsic overpolarization 
 of the alcohol groups that led to reasonable interactions with water 
 caused very strong intramolecular hydrogen bonding between vicinal 
 hydroxyls in the gas phase , depressing the energy barrier relative 
 to the neat liquid . in contrast , hvap with the drude model of glycerol agreed with the experimental value 
 within 1.0% . the cooperative dipole enhancement of the three vicinal 
 hydroxyls is an important factor in achieving proper balance between 
 intra- and intermolecular nonbonded interactions and is unique to 
 the drude polarizable model . given the similarities in the parametrization 
 protocol and validation 
 of the hexopyranose and aldopentofuranose monosaccharide force fields , the remainder 
 of this section will be devoted to a unified discussion of both . electrostatic 
 parameters were derived using the mc / sa fitting algorithm described 
 above , starting with parameters transferred from the linear polyol and cyclic ether force fields . the target data for the fitting consisted of qm dipole 
 moments to tune partial charges , atomic polarizabilities , and atom - based 
 thole screening factors . for the hexopyranoses , 237 unique conformers 
 ( encompassing 16 diastereomers ) were used as target data , while 60 
 furanose conformers were used for the aldopentofuranoses . fitting 
 was carried out using three strategies , the first based on per - isomer 
 fitting ( termed  isomerfit  ) , yielding 16 parameter 
 sets for the hexopyranoses and four sets for the furanoses . an alternative , 
 per - anomer fitting protocol ( called  anomerfit  ) employed 
 a constraint such that parameters were equivalent between the two 
 anomers of each diastereomer , thus halving the total number of possible 
 parameter combinations . the final strategy , targeting all conformers 
 in one fitting run ( called  globalfit  ) and enforcing 
 equivalency among alcohol groups in the ring to avoid overfitting , 
 yielded a single parameter set that was chosen as the final parameter 
 set for both of the carbohydrate series . while fitting on a per - isomer 
 or per - anomer basis typically yielded better overall rmsd between 
 qm and mm dipoles , these parameters were not transferable between 
 monosaccharides . the global fitting protocol produced a unified , transferable 
 parameter set the quality of which was comparable to that of the other 
 methods and represented a significant improvement over the additive 
 charmm force field , which utilizes overestimated gas - phase dipole 
 moments . torsional fitting was carried out using the same global fitting 
 strategy as applied to the additive model , targeting the gas - phase 
 energetics of an extensive conformer series ( 1846 hexopyranose conformers 
 and 1096 aldopentofuranose conformers ) . again , significant improvements 
 in the reproduction of the qm data over the additive charmm force 
 field were obtained . achieving correct conformational energetics is 
 a critical component of the carbohydrate force field , as only low 
 energy barriers may separate different ring and hydroxyl conformers ; 
 relative energetics dictate the quality of the conformational sampling 
 in the force field . the polarizable carbohydrate parameter sets 
 were validated by examining 
 a number of condensed - phase properties , including crystal volumes , 
 nmr j - coupling constants , solution densities of varying 
 monosaccharide concentration , and sugar puckering during md simulations . crystal simulations of hexopyranoses and o - methylated furanosides 
 showed an improvement over the additive charmm force field ; despite 
 the fact that both force fields slightly overestimate crystal volumes , 
 the drude-2013 force field was in decidedly better agreement with 
 experimental values . aqueous simulations of varying monosaccharide 
 concentration were in better agreement with available experimental 
 data using the drude model , differing generally by < 1% and manifesting 
 correct temperature - dependent behavior . however , underestimation of 
 the diffusion constants of monosaccharides at higher concentrations 
 was observed , a problem also present in the additive model . cooperativity 
 in neighboring dipoles of the monosaccharide alcohol groups was observed 
 in aqueous solution , a unique feature of the drude polarizable model 
 that can not be obtained with an additive model . such behavior is responsible 
 for the environment- and conformation - dependent behavior of the force 
 field . the additive and polarizable models performed similarly 
 in terms 
 of the ability to reproduce sugar - puckering behavior , but the drude 
 model is notable in its ability to qualitatively reproduce the conformational 
 tendencies of exocyclic groups on the monosaccharide rings . in the 
 development of the hexopyranose force field , hamiltonian replica - exchange ( hrex ) simulations were used to comprehensively 
 sample conformational space of the  torsion angle ( o5c5c6o6 ) , 
 while long , standard md simulations were used in the case of the aldopentofuranoses 
 to sample the o4c4c5o5 exocyclic torsion angle . the high quality of the monosaccharide force 
 field derives from strong transferability between the parameters , 
 indicating that development of higher - order polysaccharide models 
 will be possible in the context of a comprehensive biomolecular force 
 field . as of the writing of this review , work on parameters for a 
 range of disaccharides was ongoing , along with further refinement 
 of the hexapyranose monosaccharide parameters . lipid membranes 
 serve as barriers between intra- and extracellular environments , around 
 organelles , and to regulate the transport of materials throughout 
 the cell . as such , the proper thermodynamic and transport properties 
 of lipids represent an important component of a force field . the dielectric 
 environment of a lipid membrane is complex , with strongly polarizable 
 head groups with large dipole moments at the membrane  water 
 interface , glycerol linkages , and a low - polarity ( hydrophobic ) alkane 
 interior . such properties represent a considerable challenge for additive 
 force fields , and though progress has been made toward accurate models 
 of membranes , multibody effects at interfaces 
 are significant in determining the membrane dipole potential , which 
 additive force fields overestimate . by 
 conducting simulations of lipid monolayers with a first - generation 
 dipalmitoylphosphatidylcholine ( dppc , figure 9 ) using the drude polarizable model at air / water 
 interfaces , harder et al . quantified the effects of induced polarization 
 of lipid headgroups , alkyl chains , and water molecules . the polarization 
 response was found to be dependent upon both distance and orientation , 
 thus illustrating the synergy between membrane structure and electronic 
 response . interfacial water molecules were strongly polarized by the 
 lipid headgroups , which themselves gave a large contribution to the 
 membrane potential . water dipole moments relaxed back toward a gas - phase 
 value at the air / water interface , demonstrating the sensitivity of 
 the polarizable model . additive simulations with the tip3p water model 
 did not , by definition , show such a response , and the contribution 
 to the membrane potential from water was significantly lower . as a 
 result , the net difference in membrane potential between the two interfaces 
 ( v ) was approximately 2 times larger than 
 experimental values using the additive model but in quantitative agreement 
 with the polarizable model . the polarization response in the alkyl 
 chains of the lipids was a significant factor in this outcome , as 
 it buffered the dipole potential . model compounds 
 used for constituent functional groups ( trimethylammonium , tma ; dimethyl 
 phosphate , dmp ; methyl acetate ; and butane ) and for torsional fitting 
 ( phosphatidylcholine , pc ; esterified glycerol , glyc ; phosphorylated 
 glycerol , glyp ) and the full structure of the dppc phospholipid . on the basis of this foundational 
 study , and the fact that the membrane 
 potential modulates the permeability 
 of ions and the association of proteins 
 with membranes and their resulting structures and dynamics , a polarizable description of lipids is desirable 
 for simulating complex , heterogeneous biological systems . to this 
 end , a refined drude polarizable model for dppc was developed by chowdhary 
 et al . in 2013 . the dppc model was initially 
 constructed from several model compounds that represent the various 
 chemical moieties in the full lipid ( figure 9 ) , including tetramethylammonium ( tma ) and 
 dmp for the phosphatidylcholine ( pc ) headgroup , esterified glycerol 
 and its phosphorylated form to represent the linkage between the headgroup 
 and glycerol backbone , and n - alkanes for the long 
 alkyl chains . the parameters for n - heptane were transferred 
 directly from the drude alkane force field , with dihedral parameters refined to accurately account for the delicate 
 balance between gauche and trans forms of the torsions , which are critical to a proper representation 
 of the structural features in the hydrophobic interior of the membrane . 
 the linkages between these alkyl chains and the glycerol backbone 
 required the parametrization of the ester linkage , which was not previously 
 available in the force field . methyl acetate was chosen as a model 
 compound for this functional group , and the accuracy of its parameters 
 was crucial to the success of the final dppc model , as its hydration 
 properties and conformational energetics give rise to important features 
 of the membrane / water interface . good agreement was achieved over 
 a wide range of qm and thermodynamic quantities of methyl acetate , 
 including molecular polarizability , dipole moment , water interaction 
 energies , hvap , neat liquid vm , ghydr , 
 and dielectric permittivity . transfer of the methyl acetate parameters 
 to larger esters yielded good agreement with available hvap and vm data , 
 indicating that the derived parameters were suitable for use with 
 larger compounds such as lipids . the fully constructed drude-2013 
 dppc model produced good agreement 
 with experimental data and the charmm36 force field , which also reproduces many structural properties of dppc 
 membranes . among the most extensively studied properties of dppc are 
 x - ray scattering form factors , deuterium order parameters , and area 
 and volume per lipid . the drude-2013 dppc model produced form factors 
 in slightly better agreement than the charmm36 simulations , though 
 both force fields perform well in this regard . the x - ray scattering 
 form factors are produced from an analytic fit of heavy - atom electron 
 densities , which are better reproduced by charmm36 , though the differences 
 with the drude model are small . the discrepancy between these two 
 observables may be related to the quality of the analytic fit and 
 the considerable scatter in experimental data for long wavelengths . the balance of forces between the hydrophobic effect in the membrane 
 core and the hydration of the lipid headgroups and glycerol region 
 gives rise to surface tension , and thus , this quantity is a measure 
 of the quality of the balance between nonbonded terms in the force 
 field . the balance of these interactions also manifests itself in 
 the area and volume per lipid in the membrane , as well as the membrane 
 thickness . together , these quantities are indicators of the overall 
 membrane structure and the cohesive forces within it . the drude model 
 slightly underestimated area per lipid and slightly overestimated 
 the volume per lipid and membrane thickness , but all results were 
 in reasonable agreement with the range of experimental data and were 
 of comparable quality to the charmm36 force field . refinement would 
 likely improve the overall agreement , but the present model reflects 
 a viable first - generation lipid force field for dppc . deuterium 
 order parameters along the dppc structure were used as 
 specific target data in the development of the charmm36 force field , and for this reason , good agreement is obtained 
 in simulations with the additive model . despite not being directly 
 targeted , comparable agreement was observed with the drude-2013 force 
 field , but several notable deviations in the headgroup carbon atoms 
 were observed , indicating that improvements are necessary . torsional 
 scans of several rotatable bonds revealed that the ability of the 
 drude model to reproduce gas - phase qm conformational energetics in 
 model compounds was lacking , though the large number of conformers 
 that needs to be assessed remains a considerable challenge . the locations 
 of energetic minima were in agreement , indicating a reasonable model 
 for conformational sampling , but future improvements will be needed 
 to yield better agreement with deuterium order parameters as well 
 as chemical shift anisotropy . the interactions of the drude-2013 
 dppc force field with water 
 also reveal important insights into the viability and application 
 of the polarizable membrane model . the drude model produced better 
 surface tension values than charmm36 as a function of area per lipid , 
 though this improvement is due in large part to the better surface 
 tension in swm4-ndp relative to tip3p . overall , the consistency achieved 
 indicates reasonable balance between the lipid and water parameter 
 sets , indicating that future refinements targeting other experimental 
 data should be straightforward . the interaction of dppc with swm4-ndp 
 also resulted in slower lateral diffusion of the lipids , again as 
 a consequence of the slower diffusion of swm4-ndp relative to tip3p ; 
 the former water model is in better quantitative agreement with experimental 
 data , while the latter diffuses nearly 3 times too fast . the final 
 assessment of the drude-2013 dppc model was to simulate dppc monolayers 
 at air / water interfaces , as was done in work by harder et al . as in that study , the membrane electrostatic 
 potential of the drude model was in decidedly better agreement than 
 the charmm36 force field , which overestimates it by a factor of 2 . the polarizability of the alkane moieties in the hydrophobic core 
 is crucial to the agreement with experimental data . as a final 
 note , it is important to highlight the fact that the 
 drude-2013 dppc force field ( and the entire drude-2013 force field 
 described in the present review ) was developed for use with a potential - switching 
 function for lj interactions , unlike the charmm36 lipid force field , which was developed using force - switching . proper treatment of lj forces is critical to 
 obtain reproducible values for membrane properties , and all investigators should take great care to correctly 
 employ these methods . to date , the drude-2013 force field 
 has been applied to a number 
 of systems , with results from our own laboratory presented below . inherent to the utility of a force field for condensed - phase simulations 
 of large molecular systems is its computational efficiency . being 
 able to access a time scale that is biologically relevant and comparable 
 to additive simulations allows for the acquisition of novel insights 
 concerning protein dynamics and conformational ensemble properties . notably , it also allows for extensive testing of the force field to 
 identify limitations in the model , allowing for improvements going 
 forward . the drude-2013 protein force 
 field has been applied to a variety of peptides and proteins , with 
 typical simulation times ranging from hundreds of nanoseconds to the 
 microsecond scale . while force field - based simulations have led to 
 numerous insights into the nature of protein folding , they have failed 
 to capture the level of cooperativity observed in the experimental 
 folding of polypeptides and proteins . one specific example of this 
 failure is the temperature dependence of secondary structure formation . recently , we showed that the drude-2013 protein force field can reproduce 
 the experimental temperature dependence of helix formation in the 
 acetyl-(aaqaa)3-nh2 peptide , and the electronic 
 induction of the dipole moments of peptide bonds drives such a cooperative 
 folding process . the ( aaqaa)3 peptide serves as a model system in the validation and development 
 of protein force fields . the average helical content at 300 k , which 
 was determined experimentally to be around 20% , was used as the target data in developing several of the most recent 
 additive force fields , including amber03 * , amber99sb * , and charmm36 . hrex simulations using the drude-2013 force field 
 led to a helical content of 25% at 300 k , indicating good balance 
 between the helix and coil conformations that arises from the quality 
 of the underlying parametrization and electrostatic model , rather 
 than as a function of having been directly targeted during force field 
 development . most importantly , good agreement with experimental helical 
 content was achieved with the drude-2013 force field at three other 
 temperatures considered ( 280 , 320 , and 340 k , with the 320 k data 
 not reported in the original study ) , reproducing the strong dependence 
 of the helical content on temperature observed in experimental measurements 
 ( figure 10a ) . folding cooperativity 
 of ( aaqaa)3 obtained with the 
 drude-2013 protein force field . ( a ) fraction helix as a function of 
 temperature determined from nmr experiments and replica exchange simulations 
 with the drude polarizable ( blue ) and the charmm36 additive ( purple ) 
 force fields . ( b ) probabilities of the peptide being in the unfolded 
 coil , partially folded , or folded states from 300 k simulations with 
 drude ( filled blue squares ) and charmm36 ( open red circles ) force 
 fields . long helix  is defined as at least eight consecutive 
 residues being in the helical segment , forming two or more helical 
 turns . having shown that the thermodynamic 
 aspect of cooperativity was 
 captured by the drude model , as reflected in the small temperature 
 interval over the transition from the unfolded to folded state , the 
 microscopic picture of folding cooperativity was investigated , which 
 involves the presence of partially folded states at very low populations 
 along with fully , or nearly fully , folded states . as illustrated in figure 10b , the drude-2013 
 polarizable force field produced a bimodal distribution in which both 
 unfolded coil and folded ,  long helix in contrast , the additive charmm36 protein force 
 field generated a more compact , unimodal conformational ensemble with 
 the more disordered or partially folded helices . roig 
 parameters also showed the agreement with the polarizable model to 
 be improved relative to the charmm36 force field , consistent with 
 the improved treatment of cooperative helix formation . the enhanced 
 dipole cooperativity associated with the drude-2013 
 force field was shown to be due to the explicit inclusion of electronic 
 polarizability instead of any particular aspect of the force field 
 parametrization . analysis of peptide bond dipole moments , with emphasis 
 on their variation as a function of peptide backbone conformation , 
 showed the peptide bond dipole moments to be enhanced by 0.9 
 d upon formation of a helical conformation . md trajectories of drude 
 simulations can be postprocessed by reminimizing the drude particles 
 frame - by - frame after any system modifications , such as removing all 
 the water molecules , with the difference in dipole moment reflecting 
 the impact of changing the molecular environment . such analysis showed 
 that roughly one - third of this enhancement comes from the hydration 
 environment , one - third from the directly adjacent residues , and the 
 remaining one - third from the i to i + 4 hydrogen bonding , which is the unique feature of helical structure . the 
 average peptide bond dipole moment of helical residues in the ( aaqaa)3 peptide was 4.91 d in simulations with the drude-2013 force 
 field , much larger than the value of 3.87 d obtained in simulations 
 with charmm36 . the charmm36 results also show an insensitivity to 
 peptide conformation , with dipole moment distributions for -helix 
 and random coil that are overlapping ; in contrast , the drude-2013 
 force field showed conformation - dependent variations in dipole moments 
 ( figure 11 ) . similar 
 analysis was carried out for the microsecond md trajectories of several 
 globular proteins , and the results are 
 summarized in table 2 . with the drude model , residues in different secondary structures 
 carry significantly different dipole moments , while with the additive 
 charmm36 there is little differentiation . the charmm36 simulations 
 also yield much smaller fluctuations . for a fixed - charge force field 
 such as charmm36 , the variation of dipole moment is only associated 
 with the change in geometry ( for example , in the case of backbone 
 dipole , only the elongation of the c = o bond ) . the nearly constant 
 values of peptide backbone dipole moment in additive force fields 
 might limit their ability to balance secondary structure elements 
 without bias toward -helix or -hairpin . instead , the 
 explicit treatment of induced polarization in the drude-2013 force 
 field can quickly respond to local variations in the environment . dipole 
 moment distributions for different conformations of the 
 model ( aaqaa)3 peptide conducted at 300 k using the charmm36 
 and drude-2013 force fields . -helical residues were defined 
 as those found in a sequence of at least three consecutive residues 
 with ( , ) values characteristic of an -helix ( 100 
 <  < 30 , 67 <  
 < 7 ) . the peak centered at  = 4.6 d for the 
 charmm36 nonhelix represents sampling of the l conformation 
 ( which are oversampled by the force field ) and are associated with 
 a change in geometry rather than any change in dipole moment associated 
 with electronic structure . qm calculations have indicated that the peptide bond 
 dipole moment 
 is approximately 5 d , consistent with 
 the values obtained with the drude-2013 force field . it is worth noting 
 that most theoretical arguments regarding the helix macrodipole ( arising 
 from the alignment of individual peptide bonds along the -helix 
 axis ) assume a peptide bond dipole of 3.5 d. however , taking into account a mutual induction effect , it should 
 be 4.7 d , which means the macrodipole should be 35% stronger and modeled 
 with a charge separation of 0.7 e , instead 
 of 0.5 e , at the two ends of the helix that 
 is typically assumed . the larger magnitude of dipoles observed 
 with the drude model is 
 notable in itself . the parametrization of additive force fields typically 
 involves systematic overestimation of dipole moments relative to the 
 gas phase , which in the case of the charmm force fields is about 20% , 
 to capture the induction effect of the environment according to a 
 mean - field approximation . in a polarizable 
 model , the alignment of peptide dipole moments during secondary structure 
 formation leads to favorable interactions enhanced with the larger 
 dipoles , which is counterbalanced by the positive self - energy for 
 polarizing the atoms ( eq 5 ) . the drude results indicate that amino acids carry larger dipole 
 moments that are also more responsive to the complex protein environment . using the drude-2013 force field , the role of peptide - bond dipole 
 moments in the stability of -helices has recently been explored 
 in the context of the amyloid -peptide ( a ) , the central pathological protein in alzheimer s 
 disease . by simulating a fragment of a that harbors several 
 familial mutations that alter the electrostatic nature of positions 
 22 and 23 in the amino acid sequence , polarizable simulations revealed 
 that side - chain dynamics have important implications for the stability 
 of the initial -helical structure due to alterations of the 
 dipole moments of peptide bonds . for example , in the wild - type peptide , 
 electrostatic repulsion between the glu22 and asp23 side chains led 
 to a perturbation of peptide - bond dipole moments up to four residues 
 away , through the hydrogen - bonding network of the -helix . backbone hydrogen 
 bonds in the d23n and e22q mutants led to destabilization , suggesting 
 that charge charge interactions are not unique in their ability 
 to destabilize the a helical structure ; different behaviors 
 exhibited by the various mutants could perturb the peptide - bond dipoles 
 and lead to unfolding . taken together , these results suggest subtle 
 but important details that are revealed by the use of a polarizable 
 force field , as these same behaviors were not observed when using 
 the additive charmm36 force field , which suggested only nonspecific 
 unfolding phenomena . the 
 fluctuation of dipole moments in biomolecules is related to their 
 dielectric properties . the specific 
 dielectric constants computed from kirkwood g - factors were recently compared between the drude polarizable 
 model and the additive charmm36 force field . the static dielectric constants are similar for whole proteins , 
 but for their hydrophobic cores , the drude-2013 dielectrics are systematically 
 larger than those from charmm36 . this observation indicates that charged 
 residues on protein surfaces dominate the dielectric properties of 
 entire proteins . for the buried residues , additional electronic degrees 
 of freedom of the charge distribution are important , an outcome related 
 to the better description of alkane dielectric properties described 
 above . the drude model allows 
 the calculation of ensemble - averaged optical dielectric constants 
 ( inf ) , based on molecular polarizability along simulation 
 trajectories . computed inf for six proteins range 
 from 1.76 to 2.10 with an average value of 2.0 , which equals the commonly assumed value . the variation across proteins suggests that a polarizable model 
 is necessary for studying properties dominated by the protein dielectric 
 relaxation , such as reorganization energy . due to the polyanionic nature 
 of nucleic acids , modeling their structure and dynamics using empirical 
 force fields considerable effort has gone into achieving such 
 a delicate balance , and recent studies 
 using the drude-2013 dna force field have focused on analyzing dna ion 
 interactions , including the impact of ions on dna conformation , and base flipping in dna . the 
 microscopic details of the interactions between dna and monovalent 
 ions have been examined in recent studies , including a description 
 of the ionic structure in solution and 
 ion competition . structural studies 
 often suffer from an inability to distinguish between biologically 
 relevant ions ; thus , simulations provide a useful tool in rationalizing 
 ion - specific behavior , even providing evidence that can resolve conflicting 
 experimental outcomes , such as those related to groove - specific ion 
 binding . experimental techniques such 
 as buffer equilibration  atomic emission spectroscopy ( be  aes ) 
 can be used to study cation competition in the ionic atmosphere around 
 dna using mixtures of salts with a common anion . in be  aes experiments , the excess of a given cation 
 and the depletion of the anion relative to their bulk concentrations 
 quantitatively describe the ionic atmosphere around dna . these data 
 provide a very powerful metric for evaluating the performance of empirical 
 force fields in describing the electrostatic environment of nucleic 
 acids in solution . ion competition simulations were carried 
 out using the drude-2013 dna force field for the monovalent cations li , na , k , and rb , following refinement of nonbonded parameters 
 against target gas - phase qm energetic data and osmotic pressure calculations , 
 as described previously . the first assessment 
 of the ionic environment around dna was a series of volume - jacobian - normalized 
 radial distribution functions ( rdf , figure 12 ) for each of the salts , establishing a 
 baseline for the ionic structure around dna . in this regard , both 
 the drude polarizable model and the additive charmm36 force field 
 produced qualitatively similar results , in that smaller ions ( li and na ) retained their hydration shells and participated 
 in water - mediated interactions with dna , whereas the larger ions ( k and rb ) could participate in short - range ( direct ) 
 interactions with dna . these results reflect the dehydration penalties 
 for this series of ions . in comparison with cc theory , however , the 
 drude-2013 force field produced near quantitative agreement with the 
 predicted 76% charge neutralization by the monovalent salts ; the charmm36 
 force field generally underestimated the amount of charge neutralization 
 by 510% . charge neutralization is largely a nonspecific and 
 macroscopic phenomenon ; that is , dna of any sequence accumulates positively 
 charged ions in the nearby ionic atmosphere . a representation of 
 the volume jacobian defined by an equidistant 
 shell around the central base pairs of a canonical b - form dna . competition between cations was 
 determined relative to na as a background cation ; the 
 ion / ion ratio within the conventional 
 9  manning radius reflects the relative propensity of a given 
 ion to associate with dna compared to na . in agreement 
 with both be  aes and dna mobility 
 experiments , the propensity of ions 
 to occupy the ionic atmosphere around dna follows the trend of li > na > k > rb in the 
 simulations , though the use of simulation systems of traditional sizes 
 ( e.g. , solvated 10 to 12  beyond the dna ) allows the results 
 to only be interpreted qualitatively , since the anion rdf does not 
 converge to unity in small simulation systems . to expand upon these 
 findings and obtain results that are quantitatively comparable to 
 be  aes experiments , large simulation systems are required to 
 avoid ambiguities arising from the distinction between neutralizing 
 cations and those constituting additional salt . by increasing the 
 size of the simulation cell ( e.g. , box lengths of 120  ) , 
 this ambiguity is mitigated , in that the fraction of specific neutralizing 
 counterions is considerably reduced at the expense of increased computational 
 cost . such larger simulation systems reveal quantitative improvement 
 of the drude-2013 force field over charmm36 . whereas the smaller systems 
 indicated that ion / ion ratios were comparable between the drude and 
 additive models , the larger systems revealed that , in fact , the drude 
 model was better able to capture ion competition effects that are 
 in near - quantitative agreement with be this was demonstrated 
 by performing simulations at the appropriate competing ion ( li , k , and rb ) concentrations that be aes 
 experiments show yield a ratio of 2.0 for the excess background cations 
 ( na ) in the absence and presence of competing cations . the drude-2013 dna force field reproduced this value almost exactly 
 across all of these ions , whereas the charmm36 force field produced 
 considerable variation as a function of competing cation identity . these findings reinforce the physical relevance of the specific ion 
 binding described above ; that is , the 
 improvement in the microscopic details of dna  ion interactions 
 manifests into a better macroscopic description of the ionic atmosphere 
 around dna . simulations are an essential technique 
 for studying the dynamics 
 of dna in aqueous solution , due to crystal packing effects that may 
 adversely affect structures determined using x - ray diffraction and potentially conflicting outcomes from the mathematical framework 
 provided by nmr that limit is applicability in describing the global 
 conformational properties of dna . small - angle 
 x - ray scattering ( saxs ) provides a powerful means to understand the 
 conformational properties of dna in solution , as the dominant properties 
 are highlighted in a one - dimensional trace that accounts for the ensemble 
 of conformations . saxs fingerprints can be further decomposed into 
 contributions from all the moieties in a biomolecule , for instance , 
 separating the contributions from the nucleic acid bases and the sugar phosphate 
 backbone in dna . to assess the ability of the drude-2013 dna 
 force field to describe the solution structure of dna , the saxs profiles 
 of two b - form dna sequences , d(cgcgaattcgcg)2 ( ecori ) and d(ccgctagcgg)2 ( pdb i d 1dcv ) , were computed from md simulations 
 in aqueous solution containing 100 mm nacl . these calculated profiles 
 were compared to two additive force fields that are widely used in 
 the simulation of nucleic acids , charmm36 and amber parmbsc0 . the drude-2013 
 dna force field best captured the features of the saxs profile , including 
 the positions of spectral peaks that correspond to structural features 
 on different length scales . both of the additive force fields were 
 in poorer agreement , indicating deviation from the experimental conformations 
 even for global features of both dna sequences . building upon 
 these simulations , the cation identity was varied 
 to compare the effects of li , na , k , and rb on the structural properties of both dna sequences , 
 using saxs profiles to quantify these effects . the different ions 
 caused distinct differences in these saxs profiles under the drude 
 polarizable model , producing shifts in saxs peaks and small variations 
 in peak amplitude . no such response was observed in the additive charmm36 
 force field , showing that the polarizable model is more sensitive 
 to the effects of different ions . specifically , while the different 
 ions did not produce differences in the geometric properties of the 
 dna major groove , the width of the minor groove was modulated as a 
 function of cation size , with spectral shifts related to the minor 
 groove following the order of cation size : li > na > k > rb . this effect was further 
 found to be sequence - specific , with the ecori sequence 
 ( which contains an aatt tract ) being more sensitive . as no experimental 
 data yet exist to corroborate these findings , they represent a true 
 prediction of the drude polarizable force field . thus , the use of 
 the drude-2013 dna force field is clearly capable of producing testable 
 hypotheses that arise due to the many - body effects of polarization 
 between the polyanionic dna and biologically relevant cations . a mechanism for this ion - mediated perturbation of dna minor groove 
 widths was recently proposed in a study of the microscopic effects 
 of ion - mediated hydrogen bonding . specifically , 
 hydrated ions in the minor groove of dna were able to form what were 
 called  strand(1)ion strand(2 ) hydrogen - bond 
 bridges  ( sis - hbb ) , which connected electronegative atoms in 
 the paired dna strands , leading to ion - dependent reduction in minor 
 groove width . that is , smaller ions like li and na had a more profound effect than larger ions like k and rb , which are also more likely to have hydrating 
 waters dissociate due to the lower affinity of these ions for water . interestingly , this behavior was only observed with the drude-2013 
 dna force field ; ions in the additive charmm36 force field were more 
 likely to partially dehydrate , leading to direct interactions of the 
 monovalent ions instead of water - mediated interactions . the drude-2013 dna force field 
 showed that the effects of sis - hbb on minor groove width were also 
 sequence - dependent , with the shorter 1dcv sequence being more susceptible to sis - hbb 
 modulation than the 2l8q and ecori sequences . in addition to the dynamics 
 of partial dehydration , the drude-2013 
 force field produced interesting insights into the variability of 
 water dipole moments . dipole moments of waters bound to li and na , which have a large charge density , increased 
 from their bulk value ( 2.45 d ) to 2.8 and 2.57 
 d , respectively . water molecules bound to these ions that also participated 
 in sis - hbb had dipole moments that were further enhanced by 0.1 
 d. this observation shows the interplay between changes in the electrostatic 
 environment , from diffusion through bulk medium to binding with an 
 ion , and ultimately interacting with dna . finally , by subdividing 
 the md trajectories between frames in which 
 one or more sis - hbb were formed and those in which none were present , 
 it was possible to analyze dna structural properties as a function 
 of sis - hbb formation . as described above ( see section 3.3.6 ) , b - dna populates both 
 bi and bii substates , dictated by the values of backbone  and 
  dihedrals . the formation of sis - hbb in the minor groove of 
 the dna sequences examined produced an increase in the bii substate , 
 concomitant with an increase in sampling of the  north  
 sugar pucker . taken together , these findings indicate that the inclusion 
 of explicit polarization is of particular importance in characterizing 
 the structure and dynamics of dna in solution , especially in the presence 
 of biologically relevant ions . in the initial 
 validation simulation of the drude-2013 dna force field , considerable variation in nucleic acid base dipole moments was observed , 
 reflecting sensitivity to the surrounding environment . to assess whether 
 these effects had implications on high - energy , extrahelical states 
 of the bases , simulations were carried out to determine the free energy 
 change associated with base flipping under the additive charmm36 and 
 drude-2013 dna force fields . simulations 
 were carried out using the umbrella sampling technique , in which a 
 harmonic biasing potential is applied along a reaction coordinate 
 to sample high - energy states along a given pathway . in the case of 
 base flipping , the reaction coordinate was a pseudodihedral defined 
 by the center of mass of groups of atoms in consecutive nucleotides 
 in the dna structure . the charmm36 force field yielded 
 large energy barriers and thus very low equilibrium constants for 
 base - open versus base - closed states . the drude-2013 dna force field 
 produced larger equilibrium constants that were in near - quantitative 
 agreement with nmr data , an improvement of 12 orders of magnitude 
 over the additive force field . indeed 
 the local electric field that a base experiences along the reaction 
 coordinate is highly variable , and the stabilization of the open state 
 leading to the larger equilibrium constants was shown to a be a result 
 of mutual polarization between water molecules in the first hydration 
 shell around the flipped base and the base itself . the dipole moments 
 of the flipping base increased by up to 12 d upon transfer 
 of a base from its watson  crick base - paired conformation into 
 the aqueous environment , with hydrating waters subsequently responding 
 on the order of 0.1 d. while this increase is small on a per - molecule 
 basis , the cumulative effect of dipole enhancement over several water 
 molecules , especially those participating in hydrogen bonds to the 
 base , led to more favorable solvation of the base in the open states 
 in the polarizable force field versus the additive model . these observations 
 further support the notion that induced polarization is a critically 
 important component of a proper physical model and is likely relevant 
 in many biological processes . in the present review , we have summarized 
 the theory , history , 
 parametrization strategy , and applications of an empirical force field 
 based on the classical drude oscillator model , called  drude-2013  . the force field has its origins in the additive charmm force field 
 and follows much of the same parametrization approach , an iterative 
 refinement of internal and nonbonded terms to match experimental and 
 theoretical ( qm ) target data . across a broad class of chemically diverse 
 small molecules , the inclusion of explicit polarization is clearly 
 an improved representation of electrostatic properties relative to 
 traditional fixed - charge ( additive ) models . electronic response as 
 a function of environment allows for broad applicability to varying 
 chemical environments , including gas and condensed phases , as well 
 as differential behavior between solvents of different polarity . currently , as of 2016 , work is ongoing to further refine parameters 
 for proteins , carbohydrates , lipids , and dna to achieve agreement 
 with a broader array of experimental data , and a complete nucleic 
 acid force field is nearing completion that will include rna ( lemkul 
 and mackerell , work in progress ) . parametrization of linkages between 
 monosaccharides will give rise to a fully functional polysaccharide 
 force field , which can then be combined with protein and lipid parameters 
 for representation of glycoproteins and glycolipids . all of these 
 efforts are simultaneously working toward a generalized drude force 
 field that will encompass druglike molecules , analogous to the cgenff 
 force field that complements the additive 
 charmm force field . the implementation 
 of algorithms for carrying out simulations with the drude-2013 force 
 field in charmm , namd , chemshell qm / mm , openmm , and gromacs , combined 
 with available input generation servers such as the  drude 
 prepper  in the charmm - gui and 
 automated parametrization in gaamp facilitates 
 broad use throughout the theoretical chemistry community . finally , 
 we note that period of time required for the development 
 of the model to its current , 2015 , state . work was initiated in march 
 2000 , with the initial swm4-dp water model published in 2003 . work on the small model compounds and atomic 
 ions occurred over the next 8 years up to 2011 . it took an additional 23 years to produce the first - generation 
 drude-2013 polarizable force field for proteins , dna , monosaccharides , and the membrane phospholipid dppc . this lag was associated with unexpected challenges inherently associated 
 with nonadditive effects encountered upon going from small molecules 
 to the larger , polymeric macromolecules . for example , small imbalances 
 in the interactions between positively charged basic residues and 
 negatively charged acidic residues that were not treated at the small 
 molecule level can lead to unphysical interactions , including polarization 
 catastrophe . dealing with these issues , which were much less problematic 
 with the charmm additive force field in our own experience , required 
 additional time , delaying the release of the drude protein force field 
 until 2013 . however , as presented above , the availability of a computationally 
 accessible , fully polarizable force field for a range of systems including 
 biomolecules has already yielded quantitative improvements over the 
 additive force fields that are associated with the presence of polarizability , 
 allowing for novel insights on the contribution of polarizability 
 to the structure and dynamics of macromolecules to be obtained ; we 
 anticipate a large number of further interesting findings as the drude-2013 
 force field is applied more widely in the simulation community .
molecular 
 mechanics force fields that explicitly account for induced 
 polarization represent the next generation of physical models for 
 molecular dynamics simulations . 
 several methods exist for modeling 
 induced polarization , and here we review the classical drude oscillator 
 model , in which electronic degrees of freedom are modeled by charged 
 particles attached to the nuclei of their core atoms by harmonic springs . 
 we describe the latest developments in drude force field parametrization 
 and application , primarily in the last 15 years . 
 emphasis is placed 
 on the drude-2013 polarizable force field for proteins , dna , lipids , 
 and carbohydrates . 
 we discuss its parametrization protocol , development 
 history , and recent simulations of biologically interesting systems , 
 highlighting specific studies in which induced polarization plays 
 a critical role in reproducing experimental observables and understanding 
 physical behavior . as 
 the drude oscillator model is computationally 
 tractable and available in a wide range of simulation packages , it 
 is anticipated that use of these more complex physical models will 
 lead to new and important discoveries of the physical forces driving 
 a range of chemical and biological phenomena .