Implicit Solvation Using the Superposition Approximation (IS-SPA): Extension to Polar Solutes in Chloroform
IImplicit Solvation Using the Superposition Approximation (IS-SPA):Extension to Polar Solutes in Chloroform.
Peter T. Lake, a) Max A. Mattson, a) and Martin McCullagh b) Department of Chemistry, Oklahoma State University Department of Chemistry, Colorado State University (Dated: 6 October 2020)
Efficient, accurate, and adaptable implicit solvent models remain a significant challenge in the field of molecu-lar simulation. A recent implicit solvent model, IS-SPA, based on approximating the mean solvent force usingthe superposition approximation, provides a platform to achieve these goals. IS-SPA was originally developedto handle non-polar solutes in the TIP3P water model but can be extended to accurately treat polar solutesin other polar solvents. In this manuscript, we demonstrate how to adapt IS-SPA to include the treatmentof solvent orientation and long ranged electrostatics in a solvent of chloroform. The orientation of chloroformis approximated as that of an ideal dipole aligned in a mean electrostatic field. The solvent–solute force isthen considered as an averaged radially symmetric Lennard-Jones component and a multipole expansion ofthe electrostatic component through the octupole term. Parameters for the model include atom-based solventdensity and mean electric field functions that are fit from explicit solvent simulations of independent atomsor molecules. Using these parameters, IS-SPA accounts for asymmetry of charge solvation and reproducesthe explicit solvent potential of mean force of dimerization of two oppositely charged Lennard-Jones sphereswith high fidelity. Additionally, the model more accurately captures the effect of explicit solvent on themonomer and dimer configurations of alanine dipeptide in chloroform than a generalized Born or constantdensity dielectric model. The current version of the algorithm is expected to outperform explicit solvent simu-lations for aggregation of small peptides at concentrations below 150 mM, well above the typical experimentalconcentrations for these materials.
I. INTRODUCTION
The solvation of polar molecules changes the con-formations they sample and the higher order struc-tures they form relative to that found in vacuum.Self-assembling peptides, for example, achieve differentmacroscopic structure and properties under different sol-vent environments. Computational tools used to predictthis behavior are limited due to (1) the length and timescales achievable using all-atom simulations and (2) thethe lack of solvent transferability in coarse-grained mod-els. A computationally cheap, thermodynamically accu-rate, and solvent adaptable implicit solvent model willhave an immediate impact on studying these processes.Molecular simulations play an important role at dis-cerning the mechanism of aggregation and self-assemblyof peptide-based materials. Experimental approaches of-ten lack the multiscale resolution necessary to determinemacroscopic structure as well as the underlying molecu-lar driving forces. All-atom molecular dynamics (aaMD)simulations have been employed to help elucidate the lat-ter component. aaMD, however, cannot readily sam-ple the large length- and time-scales necessary to capturemacroscopic self-assembly. Additionally, these methodsstruggle with the low concentration of solutes typicallyused in self-assembly experiments. Top-down coarse-grained (CG) simulations and mesoscale models have a) These authors contributed equally to this work. b) [email protected] been used to investigate peptide aggregation but theirtie to atomistic detail is often obscured. Bottom-upCG models provide a rigorous tie between scales but of-ten lack transferability. Recent efforts have attemptedto alleviated the transferability issue but have not beenwidely adopted in molecular simulations.
Recently,a multiscale model for peptide assembly has been pro-posed that uses a combination of top-down CG modelsand aaMD simulations but a rigorous connections be-tween scales is not guaranteed. Thus, it remains anopen challenge in the field of molecular simulation toachieve computationally efficient and thermodynamicallyconsistent bridge between all-atom and CG models. Implicit solvent models provide a necessary bridge be-tween aaMD simulations and CG models and an appeal-ing platform for self-assembly prediction. Previous sim-ulation studies of peptide aggregation using implicit sol-vation are limited.
The limited application of thesemodels is due to the inaccuracies of the computationallyfeasible models and the significant computational ex-pense of more accurate models. Recent developmentshave attempted to bridge the gap between the compu-tationally expensive and cheap models but, as of yet,the Goldilocks implicit solvent model does not exist.
Additionally, much of the development of implicit solventmodels focus on water as the solvent and the portabilityto other solvents is not considered.Our recently developed implicit solvent model, im-plicit solvation using the superposition approximation(IS-SPA), provides a platform to achieve an accurateand computationally feasible implicit solvent model. This model builds on previous approaches that uti- a r X i v : . [ phy s i c s . c h e m - ph ] O c t lize the Kirkwood superposition approximation (SPA) to estimate the mean solvent force for a given soluteconfiguration. Initially developed for non-polar so-lutes in water, nothing about the underlying theory isspecific to water as a solvent. Indeed, the SPA is pre-dicted be more accurate for a less polar solvent, such aschloroform.In this paper, we extend the development of IS-SPA toaccurately and efficiently capture peptide monomer con-figurations and dimerization behavior in chloroform. Inthe subsequent theory section, we set forth the underly-ing theory of IS-SPA and how it is adapted for (1) non-spherically symmetric solvent forces and (2) polar solutesin a polar solvent. The accuracy and limitations of theapproach are demonstrated on the solvation and dimer-ization of spherical charges in chloroform followed by themodeling of the monomeric conformations and dimeriza-tion behavior of alanine dipeptide (ADP) in chloroform.
II. THEORY
IS-SPA relies on estimating the average solvent forceon atom i given the position R N of the N solute atoms,which is given exactly when the solvent is sphericallysymmetric as h f i i solv = ρ Z f i, solv ( r − R i ) g ( r ; R N ) d r , (1)where ρ is the solvent density, g is the solvent distributionfunction, and f i, solv is the solvent force. The first gen-eration of IS-SPA introduced two main approximationsto estimate Equation 1 on the fly: (1) the many-bodysolvent distribution function, g ( r ; R N ), is approximatedusing SPA with atomic radial distribution functions fitfrom the results of a single molecular configuration and(2) Monte Carlo (MC) integration is used to approxi-mate the integral. This second point yields a simula-tion with the correct force on average, with the ther-mostat being responsible for removing any excess heatproduced by the inexact calculation. This approach toan implicit solvent model was shown to capture physicsof solvent interactions for non-polar association in waterbetter than other implicit solvent models such as SASAand RISM. Subsequent sections detail the extensions ofIS-SPA to account for the non-spherically symmetric sol-vent potentials, the long range electrostatic forces, andhow the associated parameters are calculated.
A. Inclusion of Non-spherically Symmetric SolventPotentials
The previous version of IS-SPA only considered thespherically symmetric Lennard-Jones solvent forces fromthe TIP3P model of water. In the case of polar solutes ina solvent of chloroform, both the Coulomb and Lennard-Jones solvent–solute forces also depend on the orientation
FIG. 1. Schematic of a chloroform molecule with the dipolemoment labeled as p and the displacement of the center offorce of the molecule labeled as d . of the solvent relative to the solute atoms. The meansolvent force found in Equation 1 is changed to be h f i i solv = ρ Z Z f i, solv ( r − R i , Ω ) g ( r , Ω ; R N ) d Ω d r , (2)where Ω represents the internal coordinates of the solventmolecule. We use a series of approximations to analyt-ically integrate over Ω such that only the position r ofthe solvent needs to be sampled.The first approximation we introduce is to presumethat the only important internal and orientational de-gree of freedom of a solvent molecule is the alignment ofits dipole moment and that the molecule is axially sym-metric. This reduces the twelve degrees of freedom in Ω for a chloroform molecule to two, signified by ˆ p . Next,the SPA is modified to account for the orientation of thedipole moment. Namely, g ( r , ˆ p ; R N ) ’ Q Ni =1 g i ( | r − R i | ) P i ( ˆ p ; | r − R i | ) R Q Ni =1 P i ( ˆ p ; | r − R i | ) d ˆ p , (3)where P i ( ˆ p ; | r − R i | ) is the probability distribution ofthe dipole moment at a given distance from atom i .Next, we approximate the distribution of the dipolemoment to be that of a thermally ideal dipole in a con-stant electric field. Each atom is presumed to produce aradially symmetric mean field along its separation vector, E i ( r ) ˆ r . Thus, P i ( ˆ p ; | r − R i | ) ’ exp h − pT ˆ p · ˆ r E i ( | r − R i | ) i , (4)where p is the static dipole moment of a chloroformmolecule and T is the temperature in units of energy. Inthis sense, the orientation is dependent on the superpo-sition of electric fields generated by each atom. The onlyparameters that are needed for the model to predict thesolvent distribution function are the atomic radial distri-bution functions g i ( r ) and the magnitude of the effectiveelectric field from each atom E i ( r ). The latter is mea-sured in simulation through the average polarization andrelated to the effective electric field using the Langevinfunction.The solute–solvent force functions depend on all of thesolvent degrees of freedom so must be simplified in a man-ner analogous to the distribution function. We use theexpansion of the axial multipole moments for the elec-trostatic force. A description of the equations associatedwith the solvent force are found in the SupplementaryInformation. For this, the multipole moments are cal-culated for the minimum energy structure of the solventmolecule in vacuum, assuming that the distribution ofthe rotation around the dipole moment is isotropic. Afree parameter, d , is introduced relating to the positionof the molecular center relative to the carbon atom alongthe dipole of the molecule, with positive value being to-wards the chlorine atoms (Figure 1). We choose a valueof -0.21 ˚A from the carbon atom in order to minimizethe magnitude of the moments and to have the hexade-capole moment be zero, as shown in Figure SI2. We thencalculate the solvent electrostatic force to fourth order.The simplification of the Lennard-Jones potential isnot as immediately obvious. It is possible to expandthe potential analogous to the multipole expansion. Theproblem with this approach is that the functions do notconverge as quickly due to the rapidly varying potentialenergy as a function of orientation. Instead, we considerthe force to be radially symmetric, effectively truncatingthe expansion at zeroeth order. We use the same defini-tion of the solvent molecular center as with electrostaticforce. Instead of developing a function to describe theLennard-Jones force, the average force that the solventputs on each atom at a given distance along the separa-tion vector is measured and histogrammed to be used inthe simulations. B. Treatment of Long Ranged Electrostatics
The inclusion of electrostatic forces requires handlingtheir long range nature. We calculate an analytic resultfor the long ranged interactions instead of proposing anyapproximations to simplify the calculation of the elec-trostatic solvent force. We limit the discussion to non-periodic systems to avoid needing to use methods suchas particle mesh Ewald. This limitation still allows forcalculating the monomer configurations and dimer PMFsdiscussed herein.The approach to calculating the integral of the solventforce over all space is to divide the space into two parts.One part is the union of spherical volumes within somecut off distance from the solute atoms, named the interac-tion volume. The MC sampling is only performed withinthe interaction volume and the force from any MC pointis calculated for all solute atoms. Outside of the interac-tion volume, the Lennard-Jones force is set to zero, thesolvent density is presumed to be that of the bulk fluid,and the polarization density is taken to be that found ina constant density dielectric (CDD). A cut off distanceof 12 ˚A is chosen to define the interaction volume basedon when the variations in density and polarization tendto these bulk values.An analytic result for the force due to a CDD contin-uum is not immediately accessible for an arbitrary shapefound for the interaction volume created by a molecule.What is calculated readily is the result of the solventforce on atom i from the solvent polarization producedby the second j outside the interaction volume of the twoatoms. Namely, f CDD ij ( r ij ) = i = j − q q π(cid:15) (cid:0) − (cid:15) (cid:1) r ij (8 r cut − r ij )24 r ˆ r ij i = j, r ij < r cut − q q π(cid:15) (cid:0) − (cid:15) (cid:1) r i j ˆ r ij i = j, r ij > r cut . (5)This pairwise force is added to the calculation of the di-rect solute–solute forces. Since the law of superpositionis an exact relation in the case of the CDD, adding allthese pairwise interactions correctly calculates the forceoutside of the interaction volume, at the expense of alsocalculating the forces inside the interaction volume of therest of the molecule but outside that of the given ij pair.This overcounted force is sampled along with the MCsampling within the molecular interaction volume andsubtracted from the final forces. Note that this force isnon-zero even within the excluded volume of the atomssuch that the MC integration needs to sample volumesthat have zero density.This exact approach to calculating the long range forceintroduces many inefficiencies to the model, including theneed to sample within the excluded volume of the solutewhere the solvent forces are identically zero and that the force of each MC sampled point needs to be calculated foreach solute atom regardless of separation distance. Theseinefficiencies are a target for future approximations tosimplify the calculations while introducing minimal error. C. Fit Parameters for Molecular Systems
There are two parameter sets discussed in the abovetheory that need to be measured for the system: theatomic radial distribution functions and effective elec-tric fields. In the previous work on IS-SPA, it is shownthat the SPA is not valid to describe the solvent distri-bution function of a molecule built upon the radial dis-tribution functions of each separate atom. Instead, it isbetter to fit the parameters to the solvent distributionfunction observed for the molecule. The input for thismethod is achieved by running a simulation with the so-lute molecule pinned in space and sampling the solventdistribution. The mean density and polarization of sol-vent is measured in a cubic grid with cells much smallerthan the size of a solvent molecule, cubes of linear length0.25 ˚A in the case of chloroform. Instead of restrictingthe fit to some functional form with a minimal number offree parameters, the functional form of the fit is a binnedfunction of distance with bins of 0.1 ˚A.Since the underlying statistics for measuring the dis-tribution function from a molecular dynamics simulationis counting statistics, the goodness of fit function thatis minimized in the case of finding parameters for theatomic radial distribution function is related to the Pois-son distribution. Finding the most likely set of param-eters to describe the explicit solvent simulation resultsleads to minimizing the function Υ g with respect to eachfit parameter,Υ g = X n ∈ cells " N Y i =1 g i ( | r n − R i | ) − g obs n N X i =1 ln g i ( | r n − R i | ) , (6) where the outer sum is over all cells in the measureddistribution function, the inner sum and product are overall solute atoms, and g obs n is the observed value of thedistribution function of the cell. The Newton-Raphsonmethod is applied to iteratively converge to the minimumin Υ g .A different goodness of fit is needed for the case of find-ing the effective electric field generated by each atom inthe solute. The nature of the statistics associated withthe distribution of polarizations of a solvent molecule isnot as simple, but it can be approximated as being thatof a dipole in a constant external field. Finding the mostlikely set of parameters to describe the aaMD resultsleads to minimizing the function Υ EF with respect toeach fit parameter,Υ EF = X n ∈ cells N n h ln (cid:16) pT | E mod n | (cid:17) − ln sinh (cid:16) pT | E mod n | (cid:17) + pT h ˆ p i n · E mod n i , (7)where N n is the number of solvent molecules observedin the cell, h ˆ p i n is the average unit vector of the dipolemoment in the cell, and the model effective electric field E mod n is equal to E mod n = N X i =1 E i ( | r n − R i | ) | r n − R i | ( r n − R i ) . (8)In this sense, the polarization is the observable in thesimulation and the effective electric field per atom is thefit parameter. The total effective electric field in a givenpoint of space is the sum over all atomic effective fieldswhich is converted to a polarization using the Langevinfunction. Again, Newton’s method is used to iterativelyconverge to the minimum of Υ EF .Using a goodness of fit function related to the underly-ing statistics of the measurement gives the most reliableset of parameters. In fact, artifacts in the fit parameterscan be found if a typical Gaussian χ method is used,especially for small distances in the radial distributionfunctions. It is also be noted that the parameters as-sociated with atoms with identical underlying force fieldparameters are fit to a single function. The results ofthese fitting procedures are found in Figure SI3 for theADP system. III. METHODS
Explicit solvent molecular dynamics simulations alongwith the Generalized Born (GB) and constant CDD sim-ulations are performed with the AMBER 18 softwarepackage. The GB model used is the GBneck model us-ing the modified Bondi radii for the atoms. The im-plicit solvent models use an external dielectric constantof 2.3473, as measured in the from the bulk chloroformmodel.
The custom-made single ion solutes are givenLennard-Jones parameters of (cid:15) = 0 .
152 kcal/mol and r min = 7 ˚A. The parameters for alanine dipeptide arefrom the ff14SB force field. Simulations were run in theNPT ensemble at 298 K maintained using a Langevinthermostat with a collision frequency of 2 ps − and apressure of 1 bar maintained using a Berendsen barostat.A direct interaction cut off of 12 ˚A is used with particlemesh Ewald being used to account for long range elec-trostatic forces. SHAKE is employed to allow the use ofa 2 fs integration time step.The IS-SPA simulations are performed using our ownFortran code. These simulations maintain the NVT en-semble with an Andersen thermostat with a collision fre-quency of 16.67 ps − . The mass of the hydrogen atomsis set to 12 u to reduce the frequency of those bonds andallow the use of a 2 fs integration time step. No cutoff isused in the IS-SPA simulations and 100 MC points per . . . . . r (Å) g ( r ) . 00 . . . . . g ( r ) − − . .
51 0 2 4 6 8 10 12 14 16 r (Å) h ˆ r · ˆ p i − − . . h ˆ r · ˆ p i − r (Å) f L J ( k ca l / m o l / Å ) − f L J ( k ca l / m o l / Å ) q = . q = − . q = − . (d) q = . q = + . q = + . (a) q = . q = − . q = − . (e) q = . q = + . q = + . (b) q = . q = − . q = − . (f) q = . q = + . q = + . (c) FIG. 2. IS-SPA parameters for solvation from explicit solvent simulations of neutral and charged Lennard-Jones spheres inchloroform. The top row shows behavior around the neutral and positive solute and the bottom row depicts behavior around theneutral and negative solute. (a,d) The solute–solvent radial distribution function, (b,e) The solvent polarization as a functionof separation and (c,f) the mean solute–solvent Lennard-Jones force. Dashed lines in (b) and (e) are the constant densitydielectric results. atom are used for estimating the integral in Equation 1within a distance of 12 ˚A of all solute atoms.Additional simulation details are provided in the Sim-ulation Methods section of the Supplementary Material.These details include the umbrella sampling protocolsand the magnitudes of the multipole moments of chlo-roform used in the IS-SPA simulations.
IV. RESULTS AND DISCUSSIONA. Charged Lennard-Jones Solutes in Chloroform
The solvation and dimerization of neutral and chargedLennard-Jones spheres are used as test cases for the IS-SPA procedure. For the ions, point charges are placed atthe center of the Lennard-Jones spheres with | q | = 0,0.5, and 1.0 e. The solute–solvent radial distributionfunction, polarization of the solvent around the solute,and mean Lennard-Jones solvent–solute force are mea-sured from explicit solvent simulation of each of the fiveLennard-Jones spheres independently. These parametersare then used to numerically integrate the IS-SPA equa-tions for the dimerization of oppositely charged Lennard-Jones species.The solute–chloroform radial distribution functions forthe five different Lennard-Jones spheres as measuredfrom explicit solvent simulation are shown in Figure 2(a) and (d). The q = 0 e, q = 0 . q = 1 . r = 6 ˚A with g <
2, followed by an inter-stitial region, and a slight second solvation shell beforeachieving bulk density. The only major discrepancy be-tween these three species is the shoulder in the q = 0e curve (purple curve) at r = 4 . q = 0 e, q = − . q = − . q = − . r = 5 ˚A followed bya shoulder at r = 6 ˚A before the chloroform starts tobehave similarly to that found in around the q = 0 especies. The q = − . g ( r ) demon-strates even more dramatic discrepancy to the neutralwith g > This asymmetry is captured in IS-SPA sincethe g ( r )’s plotted in Figure 2 (a) and (d) are used di-rectly as parameters in the model.The polarization of chloroform around the solute isan additional parameter that is measured from explicitsolvent simulation and utilized in the IS-SPA equations.From simulation of each solute, the average orientationof the dipole moment of chloroform, ˆ p , relative to thesolute–solvent separation vector, ˆ r , is measured as afunction of solute–solvent separation distance. These R (Å) u P M F ( k ca l / m o l ) . R (Å) δ u L J P M F ( k ca l / m o l ) . R (Å) δ u C P M F ( k ca l / m o l ) .CDD + LJ | q | = q = q = + q = − q = − − − − − − −
505 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 (a) LJ q = + q = − q = q = + q = − q = − − (b) CDD q = + q = − q = + q = − − − − − −
100 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 (c)
FIG. 3. The potential of mean force (PMF) of dimerization of two Lennard-Jones spheres in chloroform. The red and blueIS-SPA curves represent the integration of the positive and negative ion forces respectively. (a) The total PMF of dimerizationfor neutral and oppositely charged ions in various models of chloroform. (b) The Lennard-Jones component of the PMF forthese same models and (c) the Coulombic component of the PMF for these same models. functions are plotted for the five charged Lennard-Jonesspheres in Figure 2 (b) and (e). Given the orientation ofthe dipole moment in chloroform depicted in Figure 1, avalue of h ˆ r · ˆ p i = 1 indicates that the hydrogen of thechloroform is pointing towards the solute and a value of h ˆ r · ˆ p i = − < r < r < r > r > r ≈ | q | = 1 e pair, plotted in Figure 3(a). The green curvesare for the neutral pair with the explicit curve in dashedlines and the IS-SPA curve in solid. The neutral specieshas a minimum in the free energy of − .
63 kcal/mol and − .
79 kcal/mol from explicit and IS-SPA respectivelyat a separation distance of 6 ˚A. Both IS-SPA and ex-plicit solvent demonstrate a slight desolvation barrier at9.5 ˚A with minimal correlations beyond this distance.The charged Lennard-Jones spheres demonstrate a largepropensity to dimerize with the explicit (dashed purple)and IS-SPA (solid blue and red) binding free energies be-ing close to 30 kcal/mol. The CDD result under predictsthe stabilization of the contact pair by approximately 5kcal/mol. The IS-SPA free energy curve for the oppo-sitely charged solutes shows high fidelity with explicitover the entire domain plotted, lending strong supportthat IS-SPA is capturing the correct physics of ion solva-tion in chloroform.The Lennard-Jones and Coulombic components of thefree energy are computed separately in the IS-SPA frame-work and can be assessed from explicit solvent using afree energy decomposition scheme. Figure 3(b) showsthe Lennard-Jones component of the PMF and Fig-ure 3(c) shows the Coulombic component of the PMF.The oscillatory nature of the Lennard-Jones compo-nent is qualitatively captured by IS-SPA; both explicit(dashed lines) and IS-SPA (solid lines) demonstrate anattractive component for the Lennard-Jones componentto the PMF of the negative ion (red) and a repul-sive component to the PMF of the positive ion (blue).
FIG. 4. Chloroform densities and polarization around a single configuration of alanine dipeptide. (a) Explicit solvent density,(b) IS-SPA density fit to explicit density, and (c) free energy difference between explicit and IS-SPA chloroform densities. (d)Explicit chloroform polarization, (e) IS-SPA polarization fit to explicit, and (f) difference in explicit and IS-SPA polarizations.
The discrepancies between IS-SPA and explicit solventLennard-Jones components of the PMF are outweighedby the dominant Coulombic contribution depicted in Fig-ure 3(c). Of particular importance is how the positiveand negative ions have slightly different Coulombic con-tributions, which is captured by IS-SPA. Additionally,the CDD result agrees at large distance, as expected,but deviates from the explicit and IS-SPA results of bothions at distances shorter than 5 ˚A.The Lennard-Jones and Coulombic components of thePMF need not be identical for the positive and negativeion, as they are not state functions, but the sum of thetwo must be equivalent. To demonstrate this, we inte-grate the mean forces from IS-SPA on the positive andnegative ions separately and plotted them in solid red andblue curves in Figure 3(a). Despite the different param-eters for the cation and anion, IS-SPA produces nearlyidentical PMFs for both species during their dimeriza-tion. That IS-SPA accurately reproduces the explicitPMF for both the positive and negative ions despite hav-ing no closure relation to restrict such a equivalence is animportant test of the method.
B. Alanine Dipeptide in Chloroform
Alanine dipeptide (ADP) is chosen as a model molecu-lar system since it is a well studied model system for sol-vent models as well as enhanced sampling methods. ADP has been mostly studied in aqueous environmentbut there have been studies of how solvent, includingchloroform, affects the observed configurations.
Asin these previous studies, we are concerned with themonomer configurations observed in a solvent of chlo-roform as quantified by the φ and ψ backbone dihedrals.The solvent density around ADP is used to fit atom-based radial density functions. For molecular systems, aparticular solute configuration is chosen and explicit sol-vent simulation is run to measure the solvent distributionaround the solute. The chosen solute configuration anda 2D slice of the 3D measured explicit solvent density areshown in Figure 4(a). In this plot, white pixels indicatebulk density, red indicates below bulk density, and blueindicates above bulk density. The excluded volume of thesolute molecule is evident as the contiguous red area inthe center of the figure. The blue to black regions justoutside the excluded volume are the first solvation shellof the molecule. Two subsequent rings of low and highdensity are evident as the solvent gets radially fartheraway from the solute followed by the noisy region of bulkdensity. The complete 3D data set is used to fit radialatomic densities (results provided in Figure SI1) usingthe Poisson regression in Equation 6. The resulting SPAfit densities are shown in Figure 4(b) and a quantitativecomparison of the relative free energies is shown in Fig-ure 4(c). Overall, agreement is observed in comparisonwith the explicit density. While there are discrepancies φ (radians) ψ (r a d i a n s ) φ (radians) φ (radians) ∆ A ( k ca l / m o l ) − π π − π π (a) − π π (b) − π π (c) FIG. 5. Ramachandran plots of alanine dipeptide in (a) vacuum, (b) explicit chloroform, and (c) IS-SPA chloroform. between the SPA and explicit, in particular in the im-mediate vicinity of the solute molecule, the free energydifferences are within the thermal energy of the system.Thus, we conclude that SPA and the fitting procedureperformed here are sufficient at reproducing the solventdensity around a given configuration of the solute.Similarly, the solvent polarization around ADP is usedto fit atom-based radial electric mean fields. A 2D sliceof the 3D chloroform polarization around ADP measuredfrom explicit solvent simulation is presented as a vectorfield in Figure 4(d). These are fit to atomic radial meanfields using the regression in Equation 7, resulting in thepolarization depicted in Figure 4(e). A difference map ispresent in Figure 4(f) with little quantitative differenceobserved. Thus, we conclude that the SPA and fittingprocedure performed here are sufficient to capture thepolarization of chloroform around this single orientationof ADP.The atom-based parameters for density and polariza-tion used to generate Figure 4(b) and (d), in combina-tion with the force functions, can be used to computethe mean solvent force as a function of ADP internalcoordinates. Unlike the results for the Lennard-Jonesspheres, the molecular systems necessitate a samplingprotocol for which we developed our own IS-SPA sim-ulation code. The results from the IS-SPA simulationare compared to vacuum and explicit chloroform simu-lations by looking at the relative free energy from eachsimulation as a function of two backbone dihedral an-gles, φ and ψ . Also known as a Ramachandran plot, thefree energy plots are depicted in Figure 5 for (a) vac-uum, (b) explicit solvent, and (c) IS-SPA. A low relativefree energy is depicted in black to purple and indicatesa high propensity for the simulation to populate thatstate; a low propensity is indicated in yellow with re-gions in white never being sampled. All three systems(vacuum, explicit, and IS-SPA) have free energy wells infour regions: C ( π < φ < − π , π < ψ < − π ), P II ( − π < φ < π < ψ < − π ), C eq7 ( − π < φ < < ψ < π ), and C axial7 (0 < φ < π , − π < ψ < eq7 witha probability of 64.4% which is depleted in both explicit(47.9%) and IS-SPA (36.7%) models of chloroform. Thedominant increase due to solvation is seen in the C pop-ulations of both explicit (32.20%) and IS-SPA (41.18%)as compared to vacuum (22.38%). Based on the proba-bility of these four states, IS-SPA captures the impact ofchloroform solvation on the ADP monomer.In addition to vacuum, explicit and IS-SPA simula-tions of ADP monomer, we also ran simulations of ADPwith a CDD and GB solvent with the dielectric constantof the solvent set to 2.3473 to match the explicit chlo-roform. The percent populations of the four states dis-cussed above are provided in Table I for CDD and GBmodels in addition to vacuum, explicit and IS-SPA. TheCDD model dramatically overstabilizes the C configura-tion (48.5%) and understabilizes the C eq7 configuration ascompared to explicit solvent. The GB model also under-stabilizes the C configuration but overtabilizes the P II configuration as compared to explicit. To more conciselycompare these distributions, we consider the relative en-tropy of the Ramachandran distribution of each model(IS-SPA, CDD and GB) as compared to explicit solvent.We include vacuum compared to explicit as a control inthe values tabulated in Table II. Of all the models, IS-SPA has the smallest relative entropy value of 0.072(5),indicating that it has the most similar Ramachandrandistribution to explicit solvent. Vacuum has a relativeentropy to explicit solvent of 0.114(5), demonstrating theIS-SPA is more similar to explicit than vacuum is to ex-plicit. Surprisingly, the CDD and GB results are actuallyworse than vacuum. C. Dimerization of Alanine Dipeptide in Chloroform
The IS-SPA parameters developed for the monomer ofADP can also be used to simulate a dimer of ADP. Toquantitatively compare with explicit solvent, we investi-
TABLE I. Percent population of monomer dihedral states in alanine dipeptide in different solvent models. The number inparentheses is the error in the last digit(s).Vacuum Explicit IS-SPA CDD GBC . . . . . II . . . . . eq7 . . . . . axial7 . . . . . S rel Vacuum 0 . . . . − − − − − R (Å) u P M F ( k ca l / m o l ) ExplicitIS-SPAIS-SPA / ExplicitvacuumGBCDD
FIG. 6. The potential of mean force (PMF) as a function ofcenter-of-mass separation. R , for the dimerization of alaninedipeptide (ADP). Five different solvation models are consid-ered: Explicit, IS-SPA, vacuum, GB, CDD. Additionally, thePMF is computed by integrating the mean force from IS-SPAalong the Explicit solvent configurations and this is denotedIS-SPA/Explicit. gate the dimerization behavior along the center-of-massseparation distance. We utilize umbrella sampling simu-lations of five different models: vacuum, explicit, IS-SPA,GB and CDD. Additionally, we compute the PMF by in-tegrating the mean force of IS-SPA using the configura-tions sampled in the explicit solvent simulations and referto this data as ‘IS-SPA/Explicit’. The resulting PMFsas a function of center-of-mass separation between themonomers are shown in Figure 6.IS-SPA correctly captures the effect of chloroform sol-vation on the dimerization of ADP. This is evident in twoaspects of the PMFs shown in Figure 6. The first aspectis the dimerization free energy of explicit and IS-SPA as compared to vacuum. ADP dimerization has a min-imum in the PMF of 5.0 kcal/mol in explicit solvent ascompared to 9.5 kcal/mol in vacuum. This demonstratesthat solvation destabilizes dimerization by 4.5 kcal/mol.IS-SPA (green curve) destabilizes dimerization of ADPrelative to gas phase by 2.7 kcal/mol. The second aspectof the PMFs that indicate IS-SPA is correctly capturingthis effect of chloroform is the position at which the PMFgoes to zero. Chloroform screens the interaction of oneADP with the other such that the attractive componentof the PMF is not present until R <
R <
10 ˚A but that thesevalues are discrepant for 5 < R <
R > < R < − . . . . .
62 0 2 4 6 8 10 12 14 16 R (Å) f ( k ca l / m o l / Å ) . − − − − − − f ( k ca l / m o l / Å ) − − − −
101 0 2 4 6 8 10 12 14 16 R (Å) f ( k ca l / m o l / Å ) − f ( k ca l / m o l / Å ) ExplicitIS-SPAIS-SPA / ExplicitGBCDD(c) Solvent Coulomb ExplicitIS-SPAVacuumGBCDD(a) Solute Coulomb ExplicitIS-SPAIS-SPA / Explicit(d) Solvent LJ ExplicitIS-SPAVacuumGBCDD(b) Solute LJ
FIG. 7. Mean forces as a function of center-of-mass separation distance, R for alanine dipeptide dimerization. The totalmean force is separated into: (a) solute–solute Coulomb, (b) solute–solute Lennard-Jones, (c) solute–solvent Coulomb, and(d) solute–solvent Lennard-Jones. Five different solvation models are considered: Explicit, IS-SPA, vacuum, GB, and CDD.Additionally, the mean solute-solvent forces from the IS-SPA model computed along the explicit solvent configurations areshown under label IS-SPA/Explicit separated into Coulomb, (c), and Lennard-Jones, (d), components. IS-SPA better reproduces the explicit dimerization ofthe PMF than the other solvation models tested hereas seen in Figure 6. Unlike Explicit and IS-SPA, GBand CDD models demonstrate finite attraction out to
R >
10 ˚A. Additionally, CDD and GB understabilizethe ADP contact dimer and oversimplify the curvaturenear the minimum. These discrepancies can be quantifiedin a single parameter, χ , defined as the integral of thesquared difference between the PMFs, χ = 1 T Z ∆∆ A ( R ) d R, (9)where ∆∆ A ( R ) = ∆ A ( R ) model − ∆ A ( R ) explicit . Thesevalues are computed for each model depicted in Figure 6as compared to explicit solvent and are tabulated in Ta-ble III. The least discrepant model is IS-SPA/Explicitwith a value of χ = 6 . < R < χ = 14 . TABLE III. Differences between PMFs shown in Figure 6 asquantified by the χ = β R ∆∆ A ( R ) d R .Model χ (˚A)IS-SPA 14 . . . (5)GB 19 . . than the CDD ( χ = 16 . χ = 19 . D. IS-SPA Algorithm Scaling
An IS-SPA simulation follows typical all-atom molec-ular dynamics simulation protocols expect for the inclu-sion of an additional IS-SPA routine to estimate the mean1solvent force on each atom. Currently, this algorithm isperformed at every step of the simulation. For a systemof N solute particles and N MC MC points per particle,the algorithm is as follows:1. Generate N · N MC MC points (sampled from pre-determined distribution).2. For each MC point, loop over all atoms and useSPA to determine density and mean field at MCpoint.3. For each atom, loop over all MC points and com-pute Lennard-Jones and Coulomb force from MCpoint on atom.The algorithm involves two loops of size N · N MC . Thefirst to compute the density and mean field at each MCpoint and the second to push the force from that MCpoint onto each atom. This is similar to what is done inGB except that we have the additional N MC points persolute atom in IS-SPA.Considering a system of N solute particles solvatedwith either M solvent atoms per solute atom in explicitsolvent or N MC MC points per solute atom in IS-SPA, weachieve the following naive (no neighbor list, no PME)performance scaling relationships for non-bonded andsolvent calculations P IS-SPA ≈ a · N MC · N + b · N (10) P Explicit ≈ b · ( M + 2) M · N + b · N , (11)where a and b are performance coefficients in front of theIS-SPA and non-bonding loops, respectively. We expect a ≥ b since IS-SPA involves two loops over all-pairs. Inpractice, we get b ≈ . a in our IS-SPA code due to theadditional algebraic steps to compute the solvent forces.Notice that IS-SPA does not have an N term in theexpected scaling relationship since each MC point is in-dependent of the other. From these scaling relationships,it is apparent that IS-SPA will perform better than ex-plicit when N MC < ba M ( M + 2). If we set N MC = 100as used in the current work and ba = 0 .
27, we find thatIS-SPA is expected to perform better than explicit forADP concentrations of lower than 150 mM in chloroform.This is significantly higher than the 54 mM high experi-mental concentrations of diphenylalanine in self-assemblyexperiments. V. CONCLUSIONS
In the current work, we adapt the previous IS-SPAmodel to accurately account for polar solutes in a po-lar, non-spherical solvent. We consider the solvent to bea dipole orienting in a superposition of mean fields em-anating from each solute atom. These mean fields andradial solute–solvent densities are determined from a sim-ulation of the monomer in explicit solvent. Combinedwith a long-ranged electrostatic term and solvent–solute Coulombic forces through the octupole term, IS-SPA sim-ulations can be performed in a procedure analogous toour previous non-polar version.Using this model, it is demonstrated that polar so-lute solvation and association is accurately captured. Asa test case, the dimerization of charged Lennard-Jonesspheres in chloroform using IS-SPA is found to be in highfidelity with explicit solvent simulations. It should alsobe noted that the asymmetry of charge solvation is builtinto the parameters fed into IS-SPA. Thus, IS-SPA cap-tures the asymmetry of charge solvation in a polar solutesuch as chloroform. More importantly, the asymmetry ofthe charge solvation still amounts to the same PMF onthe anion and cation for the dimerization of the oppositecharges. This behavior is not guaranteed due to the lackof closure of the SPA.These additions to IS-SPA also accurately capture thesolvation behavior of chloroform around alanine dipep-tide. The Ramachandran plot of the monomer is wellreplicated by the model, with IS-SPA netting the lowestrelative entropy to explicit solvent out of the three sol-vation models tested. The dimerization of ADP is alsowell captured by IS-SPA with the lowest integrated freeenergy difference relative to explicit solvent. Here thereis still room for improvement. The quantitative discrep-ancy in the dimer PMFs between 5 and 8 ˚A stem fromsubtle inaccuracies in the solvent force. We hypothesizethis is mainly due to the Lennard-Jones force since we seethat this force is less quantitative in the Lennard-Jonesspheres. We will pursue a variety of ways to improve thisincluding accounting for the non-radial component of theLennard-Jones force from chloroform.Finally, using the na¨ıve performance scaling behaviorwe determined that the current IS-SPA model should out-perform explicit solvent simulations at 150 mM peptidesolute concentration. This predicted behavior does notaccount for neighbor lists or long-ranged electrostatic cor-rections that will impact performance. Next steps in thedevelopment of the method for broad use will be to de-termine how to implement cut offs, and thus the abilityto use neighbor lists, while still accurately accounting forthe solvent forces.
SUPPLEMENTARY MATERIAL
Supporting information is provided including: simula-tion methods, multipole moments of chloroform, equa-tions for the electrostatic solvent forces in IS-SPA, andexample atomic parameters for alanine dipeptide.
AUTHORS’ CONTRIBUTIONS
P.T. Lake and M.A. Mattson contributed equally tothis work.2
ACKNOWLEDGMENTS
We gratefully acknowledge the Army Research Officefor financial support (ARO Fund number W911NF-17-1-0383). This work used the Extreme Science and Engi-neering Discovery Environment (XSEDE), which is sup-ported by National Science Foundation grant numberACI-1548562. We would specifically like to acknowledgethe San Diego Supercomputer Center Comet used underXSEDE allocation CHE160008 awarded to MM.
DATA AVAILABILITY
The data that support the findings of this study areavailable from the corresponding author upon request.All code used to parameterize and run IS-SPA for thispaper are provided in the GitHub repository: https://github.com/mccullaghlab/IS-SPAQ.git . C. L. Shen and R. M. Murphy, Biophys. J. , 640 (1995). M. McCullagh, T. Prytkova, S. Tonzani, N. D. Winter, and G. C.Schatz, J. Phys. Chem. B , 10388 (2008). O.-S. Lee, S. I. Stupp, and G. C. Schatz, J. Am. Chem. Soc. , 3677 (2011). J. Anderson, P. T. Lake, and M. McCullagh, J. Phys. Chem. B , 12331 (2018). R. Weber and M. McCullagh, J. Phys. Chem. B , 1723 (2020). X. Yan, P. Zhu, and J. Li, Chem. Soc. Rev. , 1877 (2010). R. A. Mansbach and A. L. Ferguson, J. Phys. Chem. B , 1684(2017). P. W. J. M. Frederix, R. V. Ulijn, N. T. Hunt, and T. Tuttle, J.Phys. Chem. Lett. , 2380 (2011). P. W. J. M. Frederix, G. G. Scott, Y. M. Abul-Haija, D. Kalafa-tovic, C. G. Pappas, N. Javid, N. T. Hunt, R. V. Ulijn, andT. Tuttle, Nature Chem , 1 (2015). T. Tuttle, Can. J. Chem. Eng. , 724 (2015). T. Sanyal and M. S. Shell, J. Phys. Chem. B , 5678 (2018). J. Jin, A. J. Pak, and G. A. Voth, J. Phys. Chem. Lett. , 4549(2019). X. Zhao, C. Liao, Y.-T. Ma, J. B. Ferrell, S. T. Schneebeli, andJ. Li, J. Chem. Theory Comput. , 1514 (2019). G. Colombo, P. Soto, and E. Gazit, Trends Biotechnol. , 211(2007). M. S. Shell, J. Chem. Phys. , 144108 (2008). R. Ishizuka, G. A. Huber, and J. A. McCammon, J. Phys. Chem.Lett. , 2279 (2010). R. C. Remsing, J. M. Rodgers, and J. D. Weeks, J. Stat. Phys. , 313 (2011). R. C. Remsing, S. Liu, and J. D. Weeks, Proc. Natl. Acad. Sci.U.S.A. , 2819 (2016). A. Gao, R. C. Remsing, and J. D. Weeks, Proc. Natl. Acad. Sci.U.S.A. , 201918981 (2020). P. T. Lake and M. McCullagh, J. Chem. Theory Comput. ,5911 (2017). J. G. Kirkwood, J. Chem. Phys. , 300 (1935). M. Pellegrini and S. Doniach, J. Chem. Phys. , 2696 (1995). M. Pellegrini, N. Gronbech-Jensen, and S. Doniach, J. Chem.Phys. , 8639 (1996). D. Case, R. Betz, D. Cerutti, T. Cheatham III, T. Darden,R. Duke, T. Giese, H. Gohlke, A. Goetz, N. Homeyer, S. Izadi,P. Janowski, J. Kaus, A. Kovalenko, T. Lee, S. LeGrand, P. Li,C. Lin, T. Luchko, R. Luo, B. Madej, D. Mermelstein, K. Merz,G. Monard, H. Nguyen, H. Nguyen, I. Omelyan, A. Onufriev,D. Roe, A. Roitberg, C. Sagui, C. Simmerling, W. Botello-Smith,J. Swails, R. Walker, J. Wang, R. Wolf, X. Wu, L. Xiao, andP. Kollman,
Amber 2018 (University of California, San Fran-cisco., 2018). J. Mongan, C. Simmerling, J. A. McCammon, D. A. Case, andA. Onufriev, . Chem. Theory Comput. , 156 (2007). T. Fox and P. A. Kollman, Langmuir , 8070 (1998). P. Cieplak, J. Caldwell, and P. Kollman, J. Comput. Chem. ,1048 (2001). J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E.Hauser, and C. Simmerling, J. Chem. Theory Comput. , 3696(2015). P. T. Lake and M. McCullagh, in
Intra- and IntermolecularInteractions Between Non-covalently Bonded Species , edited byE. R. Bernstein (Elsevier, Oxford, 2020) Chap. 3, pp. 71–89. W. Wu, L. Xing, B. Zhou, and Z. Lin, Microbial Cell Factories , 9 (2011). A. N. Drozdov, A. Grossfield, and R. V. Pappu, J. Am. Chem.Soc. , 2574 (2004). K. Cai, F. Du, J. Liu, and T. Su, Spectrochimica Acta Part A:Molecular and Biomolecular Spectroscopy , 701 (2015). J. Rubio-Martinez, M. S. Tomas, and J. J. Perez, J.Mol. Graph.Model. , 118 (2017). upplementary Material for: Implicit Solvation Using the SuperpositionApproximation (IS-SPA): Extension to Polar Solutes in Chloroform. Peter T. Lake, Max A. Mattson, and Martin McCullagh a)1) Department of Chemistry, Oklahoma State University Department of Chemistry, Colorado State University (Dated: 6 October 2020) a) [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] O c t . SIMULATION METHODS The model parameters used for chloroform in aaMD are for a completely flexible model.
The custom-made single ion solutes were given Lennard-Jones parameters of (cid:15) = 0 . r min = 7 ˚A and charges spanning 0.0 e to 1.0 e of equal and opposite sign.The parameters for alanine dipeptide are from the ff14SB force field. ? Explicit solventmolecular dynamics simulations along with the Generalized Born (GB) and constant CDDsimulations are performed with the AMBER 18 software package. The GB model used is theGBneck model using the modified Bondi radii for the atoms. The implicit solvent modelsuse an external dielectric constant of 2.3473, as measured in the model solvent. The IS-SPAsimulations are performed using our own Fortran code.All simulations are performed at a temperature of 298 K using a 2 fs time step in thesimulation. A Langevin thermostat with a collision frequency of 2 ps − is used in thesimulations using AMBER and an Andersen thermostat with a collision frequency of 16.67ps − is used for IS-SPA simulations. The explicit solvent simulations are run at a pressure of1 bar using a Berendsen barostat. For the explicit solvent simulations, a direct interactioncut off of 12 ˚A is used with particle mesh Ewald being used to account for long rangeelectrostatic forces. No cut offs are used in the implicit solvent simulations. Simulationsrun in AMBER use SHAKE to restrain the bonds containing hydrogen atoms in the solutemolecule while the IS-SPA code allows for fluctuations in these bonds but changes the massof the hydrogen atoms to 12 u to reduce the frequency of those bonds.Initial frames for all the implicit solvent simulations are collected from the explicit solventsimulations. The input for the explicit solvent simulations are produced using tleap. Theenergy of the initial configuration is first minimized for 20000 steps. The system is thenheated to 298 K over 50 ps and then the volume is allowed to change for another 50 ps. Thesystem is evolved for 3 ns before collecting data in order to remove any artifacts from theinitial configuration.The solvent forces in the IS-SPA simulations are calculated by MC integration of theintegral in Equation (1) from the Theory section. The focus of the current study is not onthe efficiency of calculating this integral. Instead, the sampling is performed in each stepwith 100 MC sampled points per solute atom uniformly chosen in a spherical volume arounda solute atom with a 12 ˚A radius. The long range electrostatic force calculation discussed2n the Theory section requires sampling even within the excluded volume of the solute. Forthe electrostatic force, the axial multipole moments of chloroform are calculated using theminimum energy structure of the model system at a position of d = − .
21 ˚A as described inthe Theory section. These values are 0.2291 e˚A for the dipole moment, 0.09275 e˚A for thequadrupole moment, 0.5140 e˚A for the octupole moment, and 0 e˚A for the hexadecapolemoment. The calculation of the electrostatic forces is truncated at the octupole term.The parameters for IS-SPA of the single Lennard-Jones ion are calculated from a simula-tion of a single ion solvated with 33,000 chloroform molecules resulting in a cubic box withan approximate linear length of 165 ˚A. A total of 10 ns of simulation is run for each systemand the average solvent density, polarization, and Lennard-Jones force is measured in 0.1 ˚Abins.The mean forces and PMFs of the Lennard-Jones ions are calculated by performingumbrella sampling (US) simulations as a function of the separation distance between thetwo solutes of equal and opposite charges. The two ions are solvated with 1,800 chloroformmolecules resulting in a cubic box with an approximate linear length of 62 ˚A. A harmonicrestraint with a spring constant of 20 kcal/mol/˚A is applied to the two ions. The directinteraction between the two ions is not calculated allowing for calculation of the PMFbetween 0 ˚A and 25 ˚A, using windows separated by 0.5 ˚A. Each window is simulated for 100ns. Only the explicit solvent model is simulated since the IS-SPA results can be obtainedby direct integration of the mean force for a given ion separation.The distribution of internal orientations of alanine dipeptide is calculated by simulatinga single molecule. For the explicit solvent simulation, the molecule is solvated with 5,000chloroform molecules resulting in a cubic box with an approximate linear length of 65 ˚A. Abias potential is applied to the φ -dihedral angle in order to fully sample the dihedral statesof the molecule. The bias potential is found by calculating the minimum energy states ofthe system in vacuum in all φ - and ψ -dihedral space and finding the Fourier decompositionto third order. The resulting potentials have spring constants of k = 3 . k = 3 . k = 3 . ψ = 4 . ψ = 1 . ψ = 3 . is applied to the two solute molecules. The PMF is calculated via windowsbetween 3.5 ˚A and 16.5 ˚A separated by 0.5 ˚A. Each window is simulated for 100 ns. II. MULTIPOLE MOMENTS OF CHLOROFORM − − . − − . . . − − . . d (Å) m o m e n t ( e Å ℓ ) ℓ = ℓ = ℓ = ℓ = ℓ = FIG. 1. Magnitude of the multipole moments for chloroform as a function of distance from thecarbon atom. A positive value is towards the chlorine atoms. II. ELECTROSTATIC SOLVENT FORCES IN IS-SPA
The force of the chloroform solvent electrostatic force is expanded in a multipole expan-sion. The assumption is that the distribution of the dipole moment of a chloroform moleculeat a given point in space is that of an ideal dipole. Given that distribution, the averageforce is then expanded. The force is calculated as f Csolvent = ∞ X ‘ =1 q i M ‘ h P ‘ ( ˆ p · ˆ r ) i π(cid:15) r ‘ +2 (cid:16) P ‘ +1 ( ˆ E · ˆ r ) ˆ r − P ‘ ( ˆ E · ˆ r ) ˆ E (cid:17) , (1)where q i is the charge on the solute atom, M ‘ are the multipole moments, ˆ r and r are theunit vector and distance between the solute and solvent, ˆ E is the unit vector of the electricfield, and P ‘ are the Legendre polynomials. The angle brackets is the ensemble average overthe distribution of dipole moments. This requires calculating h ( ˆ r · ˆ p ) ‘ i for a dipole in auniform electric field, which is calculated recursively as h ( ˆ r · ˆ p ) ‘ i = − ‘TpE h ( ˆ r · ˆ p ) ‘ − i + ‘ even pET ‘ odd , (2)with h ( ˆ r · ˆ p ) i = 1. IV. ATOMIC PARAMETERS FOR ALANINE DIPEPTIDE − − . . . .
53 0 2 4 6 8 10 12r (Å) A ( k ca l / m o l ) − − − − E M F ( k ca l / m o l / e / Å ) − f L J ( k ca l / m o l / Å ) (a) (b) (c) FIG. 2. Plots of the atomic parameters for each atom type in alanine dipeptide. a) The freeenergy of the solvent to a given atom. b) The effective electric field from a given atom. c) Theaverage Lennard-Jones force between a given atom and the solvent. EFERENCES T. Fox and P. A. Kollman, Langmuir , 8070 (1998). P. Cieplak, J. Caldwell, and P. Kollman, J. Comput. Chem. , 1048 (2001). D. Case, R. Betz, D. Cerutti, T. Cheatham III, T. Darden, R. Duke, T. Giese, H. Gohlke,A. Goetz, N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T. Lee, S. LeGrand,P. Li, C. Lin, T. Luchko, R. Luo, B. Madej, D. Mermelstein, K. Merz, G. Monard,H. Nguyen, H. Nguyen, I. Omelyan, A. Onufriev, D. Roe, A. Roitberg, C. Sagui, C. Sim-merling, W. Botello-Smith, J. Swails, R. Walker, J. Wang, R. Wolf, X. Wu, L. Xiao, andP. Kollman,
Amber 2018 (University of California, San Francisco., 2018). J. Mongan, C. Simmerling, J. A. McCammon, D. A. Case, and A. Onufriev, . Chem.Theory Comput.3