FULL PAPER

WWW.C-CHEM.ORG

Polarizable Coarse-Grained Models for Molecular Dynamics Simulation of Liquid Cyclohexane Oliwia M. Szklarczyk,† Eirini Arvaniti,‡ and Wilfred F. van Gunsteren* Force field parameters for polarizable coarse-grained (CG) supra-atomic models of liquid cyclohexane are proposed. Two different bead sizes were investigated, one representing two fine-grained (FG) CH2r united atoms of the cyclohexane ring, and one representing three FG CH2r united atoms. Electronic polarizability is represented by a massless charge-on-spring particle connected to each CG bead. The model parameters were calibrated against the experimental density and heat of

Introduction Investigating biological systems with the aid of biophysical methods on the nanoscale level of spatial resolution presents a great challenge to the scientific community. The development of computational techniques to simulate such systems has been an ongoing process for many decades,[1] since the properties of Lennard–Jones liquids were investigated with the aid of molecular dynamics (MD) simulation.[2] Very accurate models are often costly in terms of computational power. For example, quantum-chemical calculations in which molecular systems are treated at the level of nuclear and electronic degrees of freedom are prohibitively costly when applied to a protein or chain of DNA, and only very short time scales and small systems can be considered using this technique. Approaches to develop approximate methodologies that would allow for faster simulation of large biomolecular systems have been the focus of many research groups in the field. Many studies aim at replacing the available atomic-level biomolecular force fields by supra-atomic ones based on simpler and yet possibly physically correct representations of biomolecules, such as proteins,[3–11] DNA,[12–15] or lipids.[16,17] These approaches include both multilevel resolution representations of systems, that is, treating part of the system in a more detailed manner, usually the part which is crucial to some biomolecular process, for example, the quantummechanical/molecular-mechanical (QM/MM) approach[18–20] which can be used to treat active sites of enzymes on the quantum-chemical level and the remainder of the protein and the solvent on a classical–dynamical level, as well as coarsegrained (CG) supra-atomic or supramolecular models with beads representing sets of atoms or even sets of molecules.[21] Coarse-graining is an efficient approach because of the possibility of reducing the number of degrees of freedom in the system.[22] This may lead to a loss of essential interactions, for example, hydrogen bonding between water and proteins,[23] which results in an overstabilization of the protein conforma-

vaporization of liquid cyclohexane, and the free energy of cyclohexane hydration. Both models show good agreement with thermodynamic properties of cyclohexane, yet overestimate the self-diffusion. The dielectric properties of the polarizC 2015 Wiley able models agree very well with experiment. V Periodicals, Inc. DOI: 10.1002/jcc.23929

tion. When coarse-graining a system, it is thus crucial to maintain a realistic representation of the interactions while simplifying it as much as possible to maximize the computational efficiency. To be able to simulate large systems and still keep the essential degrees of freedom, one may develop CG models for the solvent that reproduce structural properties, such as density and thermodynamic properties, for example, heat of vaporization and excess free energy, of the pure liquid, as well as dielectric properties, for example, static dielectric permittivity, which governs the screening of long-range Coulomb interactions which depend on and thus can influence the conformations adopted by the solute molecule. Using CG solvent models, one can simulate larger systems with an abundance of solvent molecules, which allows for a reduction of finite-size effects. The importance of explicitly including polarization effects has broadly been discussed in the literature, with examples in several review articles.[24–28] Protein stability and function is influenced by the polarizable environment.[29–35] Electronic polarizability affects the conformational properties of peptides near membranes[36] and plays a role in ion transport through channels. Yet, the majority of available molecular force fields [a] O. M. Szklarczyk, Eirini Arvaniti, Wilfred F. van Gunsteren Laboratory of Physical Chemistry, Department of Chemistry and Applied Biosciences, Swiss Federal Institute of Technology ETH, 8093, Z€ urich, Switzerland E-mail: [email protected] Contract grant sponsor: The National Center of Competence in Research (NCCR) in Structural Biology of the Swiss National Science Foundation; Contract grant number: 200020-137827; Contract grant sponsor: The European Research Council; Contract grant number: 228076 † Present address: Department of Electrical Engineering and Information Technology Swiss Federal Institute of Technology ETH, 8006, Z€ urich, Switzerland ‡ Present address: Department of Biology Swiss Federal Institute of Technology ETH, 8093, Z€ urich, Switzerland C 2015 Wiley Periodicals, Inc. V

Journal of Computational Chemistry 2015, 36, 1311–1321

1311

FULL PAPER

WWW.C-CHEM.ORG

and MD software packages accounts for electronic many-body effects in a mean-field manner by enhanced charges and a pairwise additive Coulomb potential energy function. Such an implicit treatment of electronic polarization may well be justified for spatially homogeneous systems such as neat liquids, but may be a poor approximation when considering spatially electronically inhomogeneous systems containing proteins, membranes, interfaces between different compounds, or mixtures of solvents. To include electronic polarizability explicitly in force fields, mainly three approaches are used: the fluctuating charge (FQ) model,[37–39] the inducible point-polarizable dipole (PPD) approach,[18,40] and the charge-on-spring (COS) model.[41,42] In the FQ model, the polarizability is introduced by changing the size of the atomic partial charges in response to changes in the local electric field. The fluctuating charges can be assigned fictitious masses and treated as additional dynamic degrees of freedom, or their values may be obtained at every time step of the simulation by minimization of the electrostatic energy of the system for the given configuration. The FQ approach was applied to obtain polarizable finegrained (FG) models for various compounds,[43] for example, aliphatic hydrocarbons[44–47] or lipids and membrane proteins.[48,49] In the PPD approach, each polarizable atom or site is assigned an induced ideal dipole moment proportional to the electric field at the site and to the atomic or site polarizability. The electric field can be determined through a self-consistent field calculation in which the field equations are solved iteratively until a certain level of convergence is achieved. An alternative method to compute the electric field is the extended Lagrangian method,[50] in which the polarizable degrees of freedom are included as dynamic variables. An example of a polarizable force field for FG alkanes based on the PPD model can be found in Ref. [51]. The COS model has been applied to develop several force fields for small compounds such as water,[52,53] ions,[54] ethylene glycol,[55] methanol,[56] alcohols,[57] carbon tetrachloride,[58] chloroform,[59] ethers,[60,61] mercaptans,[62] or hydrocarbons,[63,64] such as cyclohexane,[65] as well as lipid membranes[66,67] or nucleic acids.[68,69] The COS model allows to mimic a redistribution of the electron density around a given nuclear configuration due to the local electric field, and generally assumes a linear response of the induced dipole moment lvi of the COS site vi to the electric field Ei at this site, lvi 5ai Ei ;

(1)

where ai is the polarizability of site i. The induced dipole moment is represented by two point charges of equal size and opposite sign connected via a spring, one located at an atom or a fixed site in a molecule, and the other a massless virtual particle v. The force constant of the spring is determined by the given atom or site polarizability, and charge and virtual charge do not interact with each other through a Coulomb force, only through the spring force. The COS model has the very useful feature of avoiding the complex evaluation of 1312

Journal of Computational Chemistry 2015, 36, 1311–1321

dipole–dipole interactions in the PPD approach—only pairwise Coulombic interactions are calculated. Thus, it is straightforwardly implemented in most (bio-)molecular force fields and software. The computational cost of the additional chargecharge interactions can be kept low by assigning COS sites to a subset of atoms, for example, only to atoms that are strongly polarizable, as is done in the work presented here. The induced COS dipoles are calculated from the total electric field at their charge sites which is determined self-consistently at each time step.[70] Each atom i with a permanent charge qi that serves as a COS site has in the COS model a charge equal to qi 2qvi , where qvi is the charge of the massless COS particle (28e). The displacement Drvi of the massless virtual particle from the COS site determining the induced dipole moment is calculated from the electric field Ei , Drvi 5rvi 2ri 5

ai Ei : qvi

(2)

The electric field can be calculated from the derivative of the electrostatic potential with respect to atomic positions ri , Nat X

1 ECi ðrN Þ5 4pE0 Ecs

j 6¼ i

  3 2  qj 2qvj ri 2rj qvj ri 2rvj 4 5; 1 jri 2rvj j3 jri 2rj j3

j inside cut-off i j not excluded

(3) where 0 is the dielectric permittivity of vacuum, cs is the dielectric permittivity inside the cut-off sphere, Nat denotes the number of atoms or COS sites in the system, and i and j refer to atom i and atom j, respectively. To obtain the electric field at positions ri of the atomic centers an iterative procedure for solving eqs. (2) and (3) is usually used, in which the COS massless particle positions rvj are optimized until the electric field fulfills a convergence criterion. The Coulombic part of the electrostatic term in the potential energy function becomes V C ðrN ; sM Þ5

at 21 1 NX 4pE0 Ecs i51

Nat X

j>i j inside cut-off i j not excluded

2 4

  qi 2qvi qj 2qvj jri 2rj j

 1

 qi 2qvi qvj 1 jri 2rvj j

  qj 2qvj qvi jrvi 2rj j

(4) 3

1

qvi qvj 5: jrvi 2rvj j

The three described methods for explicit inclusion of electronic polarizability are mostly used in all-atom and unitedatom force fields. Polarizable CG force fields are relatively recent and include mostly models for water.[71–74] In the work presented here, we develop two polarizable CG force field parameter sets for cyclohexane, a molecule used as solvent in industry,[75] in studies involving proteins,[76] or in experiments involving living organisms.[77] A nonpolarizable CG model of cyclohexane can be found in Ref. [78]. It has three so-called WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

dielectric permittivity of the COS FG model agrees very closely with the experimental value.

Parameterization strategy of the CG models

Figure 1. CG models of cyclohexane: 3-bead model pCG3 (left, shown in green) with each bead representing two CH2r united atoms (shown in black) and 2-bead model pCG2 (right, shown in orange) with each bead representing three CH2r united atoms.

SC1 beads, each approximating two CH2 ring groups. The nonbonded interaction of the CG molecule is represented by a purely Lennard–Jones potential energy function. Because of its nonpolar nature dielectric properties are not represented, as in most FG models of this compound. The polarizable CG models of cyclohexane developed here are compatible with the GROMOS force field.[79–81] Two sizes of the CG beads are investigated: a smaller one used to model cyclohexane by three CG beads, as in Ref. [78], and a larger one used to model cyclohexane by two beads (Figure 1). The bonded and nonbonded van der Waals force field parameters of these models are calibrated against experimental values for the density and heat of vaporization of liquid cyclohexane, and its Gibbs free energy of hydration. Subsequently, each of the two models was extended by explicitly including electronic polarizability. Explicit electronic polarization of CG cyclohexane is modeled by a COS massless virtual charge added to each CG bead. Neither off-atom sites[70] nor damping of the electric field[82] was used.

Methods All simulations presented here were performed with the aid of the GROMOS11 software package for biomolecular simulation.[83,84] The analysis of the simulations was done using the GROMOS11 software package.[85,86] FG united (CH2) atom level models of cyclohexane Simulations at 298.15 K and 1 atm of previously developed FG (united-atom [CH2] level) nonpolarizable and polarizable models of cyclohexane[64] were used as a basis for the development of the CG models of cyclohexane. The FG cyclohexane consists of six GROMOS united atoms of type CH2r[79] connected by rigid bonds of length bFG 0 . The polarizable COS FG model of cyclohexane reported in Ref. [64], which will be further referred to as the pFG6 model, consists of united atoms of type pCH2r which have the same van der Waals C6 and C12 parameters as the nonpolarizable FG model,[79] and each of the united carbon atoms serves as COS site with polarizability a corresponding to the value for an aliphatic CH2 group in vacuum taken from Ref. [87]. Both the FG and COS FG models reproduce structural and thermodynamic experimental properties of liquid cyclohexane very well.[64,79] Additionally, the static

Two polarizable CG models were studied: a two-bead model of cyclohexane with each bead representing three CH2 groups, which will be further referred to as the pCG2 model, and a three-bead model with each bead representing two CH2 groups, which will be further referred to as the pCG3 model. The CG cyclohexane models were developed in a stepwise manner (Figure 2). The CG initial configurations were created based on single-molecule configurations of nonpolarizable FG cyclohexane. The position of a bead was taken as the center of geometry of the two or three CH2 groups that were to be CG into one CG bead. Bead masses (Table 1) were chosen equal to the sums of the atomic masses of the atoms represented by the beads. The CG bond lengths bCG 0 were taken from FG simulations by mapping the configurations onto CG cyclohexane molecules CG3 and CG2, and calculated as the average over the bCG values of the FG to CG mapped distribution. In the GROMOS force field, the van der Waals interactions between atoms are represented by a Lennard–Jones potential energy term V LJ , for an atom pair or bead pair (i, j), VijLJ ðrij Þ5

! "   6 # C12 ði; jÞ C6 ði; jÞ rij 12 rij 2 6 2 54‹ij ; rij rij rij12 rij

(5)

where ij describes the depth of the V LJ minimum, and rij corresponds to a zero LJ potential energy. Interactions between different types of atoms or beads i and j are defined by LJ repulsion and dispersion parameters obtained using the combination rules pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C12 ði; iÞC12 ðj; jÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C6 ði; jÞ5 C6 ði; iÞC6 ðj; jÞ:

C12 ði; jÞ5

(6)

In the rigid CG cyclohexane models, no intramolecular van der Waals interactions or Coulomb interactions are evaluated. The initial values of the LJ force field parameters for pairs (I, J) of CG beads were estimated from the FG-to-CG mapping as an ensemble average, VIJLJ ðrIJ Þ5

1 2

X NI X NJ

VijLJ ðrij Þ ;

(7)

i51 j51

where NI is the number of atoms i per bead I and VijLJ ðrij Þ is calculated from eq. (5). In the absence of an external electric field, a liquid of polarizable molecules without permanent partial charges has zero electrostatic interactions. However, in mixtures with polar molecules, a neutral, polarizable molecule with polarizability a will acquire an induced dipole moment [eq. (1)]. The polarizability a of each CG bead was calculated as the sum of the atomic polarizabilities of CH2 aliphatic groups in vacuum, that is, the Journal of Computational Chemistry 2015, 36, 1311–1321

1313

FULL PAPER

WWW.C-CHEM.ORG

Figure 2. Parameterization strategy for developing a polarizable CG model of a compound.

sum of polarizabilities of the FG united atoms in the COS FG model of cyclohexane (pFG6[64]), apCG2

5

2apCH2

apCG3

5

3apCH2

(8)

Different properties of liquid cyclohexane were used for calibration of the nonbonded interaction parameters of the CG models at steps 3–5 of the parameterization procedure (Figure 2). Simulation set-up Unless otherwise specified, the simulations were performed at constant temperature and pressure. Cubic boxes containing 1024 CG cyclohexanes were generated at the experimental density at 298.15 K and 1 atm and the motion was simulated under periodic boundary conditions. The temperature was kept constant at a reference value Tref 5 298.15 K using weak coupling to a temperature bath[88] with a coupling time of 0.1 ps, separately for the translational and rotational degrees of freedom. The pressure was kept at 1 atm by weak coupling to a pressure bath with a coupling time of 0.5 ps. MD simulations.

The isothermal compressibility was set to the experimental value for cyclohexane.[89,90] The SHAKE algorithm[91] was applied to maintain the geometry of cyclohexane with a relative geometric tolerance of 1024. The integration time-step was 2 fs. We used a molecular virial and a charge-group based nonbonded interaction pair list. Each CG bead was treated as a separate charge group. In pure FG simulations, triple range cut-off radii of 0.8/1.4 nm were used to treat van der Waals and electrostatic interactions. These values have been chosen[92] such that the cut-off error in the van der Waals energy is maximally 1%. A twin-range cut-off radius of 2.0 nm was used in the simulations of pure CG systems, and triple-range radii of 1.4 and 2.0 nm were used when simulating mixed FG/ CG systems. When using supramolecular CG beads, a longer cut-off radius is required to keep the cut-off error in the van der Waals energy low.[74] The interactions within the shortrange cut-off were calculated every time step, while the intermediate-range interactions were kept constant for five steps. The long-range interactions beyond 2.0 nm were approximated by a reaction field with dielectric permittivity 2.028, the dielectric permittivity of the FG polarizable cyclohexane model FG6.[64] The center-of-mass motion was removed

pffiffiffiffiffi pffiffiffiffiffiffiffi Table 1. Values of the force field parameters: ideal bond length b0, Lennard–Jones dispersion C 6 and repulsion C 12 parameters for interactions with pffiffiffiffiffiffiffiffiffiffiffiffiffi nonpolar atom types, Lennard–Jones repulsion parameter C 12 ð2Þ for interactions with CG polarizable water (CGW[74]) water, and particle or bead charge q and polarizability a.

System pFG6 pCG3 pCG2

Bead Composition

Bead mass (amu)

0

b (nm)

1 pCH2r 2 pCH2r 3 pCH2r

14.027 28.054 42.081

0.153 0.218 0.194

q (e)

a ((4p0) 1023 nm3)

pffiffiffiffiffi C6 ((kJ mol21 nm6)1/2)

pffiffiffiffiffiffiffi C12 (1023 (kJ mol21 nm12)1/2)

pffiffiffiffiffiffiffiffiffiffiffiffiffi C12 ð2Þ (1023 (kJ mol21 nm12)1/2)

0.000 0.000 0.000

1.835 3.670 5.505

0.085640 0.188194 0.312803

5.2970 16.1401 36.4096

5.56185 40.35025 91.02408

Shown are values for the polarizable FG model of cyclohexane (pFG6) consisting of 6 pCH2r united atoms, for the CG cyclohexane model of bead size 2 (pCG3) and for the CG cyclohexane model of bead size 3 (pCG2).

1314

Journal of Computational Chemistry 2015, 36, 1311–1321

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

every 20 ps to avoid generation of motion of the whole system due to numerical noise. For evaluating the LJ interaction parameter sets, we performed 0.5 ns equilibration followed by 10-ns MD production simulation. For the parameter set, best reproducing the experimental density and heat of vaporization, the simulation was prolonged for another 10 ns. The configurations and energies during all the CG simulations were saved every 5 ps for analysis. Thermodynamic integration calculations. The Helmholtz excess free energy of liquid cyclohexane was calculated using thermodynamic integration (TI) by switching off all nonbonded interactions.[64] The system was changed from a fully interacting pure liquid of Nmol 5 1024 cyclohexane molecules in a cubic box (k 5 0) to 1024 molecules without nonbonded interactions (k 5 1), at constant temperature and volume. The Gibbs free energy or free enthalpy of solvation DGsolv of a cyclohexane molecule in a box of water and of a water molecule or a methanol molecule in a box of cyclohexane were also calculated using TI.[64] The systems were changed from one solute molecule in a box of solvent (k 5 0) to one solute molecule without nonbonded interactions in a box of solvent molecules (k 5 1) under isothermal–isobaric conditions. The dependence of the Hamiltonian on the coupling parameter used to switch off or on particular interactions is given in Refs. [93,94]. The free enthalpy of solvation DGsolv was calculated for cyclohexane in water and for water in cyclohexane combining the polarizable FG and CG models for cyclohexane with the polarizable FG (COS/G2[53]) model and a polarizable CG (CGW[74]) model for liquid water. For the calculations of solvation of methanol in cyclohexane, we used either the FG nonpolarizable methanol model (“MeOH” from the 45A3 GROMOS force field[79]) or the polarizable FG methanol model (“COSM”[56]) in combination with the CG liquid cyclohexane models. Single solute molecules were solvated in a box of solvent. The pFG6 cyclohexane solute molecule was solvated in 2043 FG water molecules. The pCG2 cyclohexane molecule was solvated in 489 CGW water molecules, and the pCG3 cyclohexane was solvated in a box of 491 CGW water molecules. In case of solvation of water or methanol in cyclohexane, each system consisted of a single water or methanol molecule or water bead and 1023 cyclohexane molecules. The boxes were thermalized for 200 ps with a linear increase of temperature from 1 to 298.15 K, and subsequently equilibrated at constant temperature and volume for 200 ps. Finally, a 400 ps equilibration was performed at constant temperature and pressure for each system. The isothermal compressibility was set to the experimental value. Triple-range cut-off radii were used for van der Waals and electrostatic interactions for mixed FG/CG systems with the short-range cut-off of 1.4 nm and the longrange cut-off equal to 2.0 nm. The interaction was updated every step within the short-range cut-off radius and every five steps in the intermediate range. In the case of pure CG systems, a twin-range cut-off of 2.0 nm was used. The long-range interactions were approximated by a reaction field with the

dielectric permittivity of the medium outside the cut-off radius equal to the value of each model solvent, 2.028 for cyclohexane,[64] 87.8 for COS/G2 water,[53] and 78.5 for CGW water.[74] The dielectric permittivity of the cut-off sphere, CG , for calculating CGW–CGW interactions was set to 2.5, and the dielectric permittivity for calculating the interactions between pCGn and CGW water, mix , was set to 2.3.[95] In the TI simulations, 21 k-values were used, 0.05 apart in the range [0, 1]. At each k-value, 100 ps equilibration and 2 ns sampling were performed. The energies and configurations were saved every 2 ps. The systems were weakly coupled to a temperature bath[88] with a coupling time of 0.1 ps, separately for the solute and solvent. We used a two-step method for the TI calculations of the polarizable systems, in which the electrostatic interactions and the Lennard–Jones interactions were switched off separately. In a first step (first TI simulation from k 5 0 to 1), no soft-core interaction was used (aLJ 5 0, aC 5 0), and only the electrostatic interactions were switched off. In a second step (second TI simulation from k 5 0 to 1), the Lennard–Jones interactions were switched off with a softness parameter aLJ 5 0.5.[92,96] After integrating over h@H=@ki, the two integrals were added. Analysis The enthalpy of vaporization DHvap , the self-polarization energy Upol, the radial distribution function g(r), the averageinduced molecular dipole moment lind;mol , the static dielectric permittivity (0), the self-diffusion constant D, the surface tension c, the isothermal compressibility jT, the thermal expansion coefficient aT, and the shear viscosity g were calculated as described in Ref. [64]. Heat capacity Cp. The isobaric molar heat capacity was calculated using a finite-difference approximation,[97,98]

Cp 5Cpsim 1dCp ;

(9)

where Cpsim is calculated from the simulations as Cpsim 5

Utot;2 2Utot;1 : T2 2T1

(10)

Utot is the total energy per mole, T is the temperature, and dCp is an empirical correction to the heat capacity calculated per functional group. To evaluate the first term in eq. (9), for each system, 10 ns constant pressure simulations (preceded by a 1 ns equilibration) were performed at two additional temperatures: 278 and 318 K. The temperatures were chosen such that the considered systems would remain in the liquid phase. Thus, the systems were simulated at three different temperatures, from which two values of the heat capacity were obtained. The final value of the first term in eq. (9) was the average of the two intermediate values. The second term in eq. (9), dCp , is a first-order empirical correction for the neglect of hydrogen atoms in the united atom model, for the rigid treatment of bonds within a molecule, as well as for the neglect of quantum-mechanical features of the Journal of Computational Chemistry 2015, 36, 1311–1321

1315

FULL PAPER

1316

Shown are values for the polarizable FG model of cyclohexane (pFG6) consisting of 6 pCH2r united atoms, taken from Ref. [64], values for the polarizable CG cyclohexane model of bead size 2 (pCG3) and values for the polarizable CG cyclohexane model of bead size 3 (pCG2). Experimental values are shown in the first row. The heat capacity values Cpsim are those obtained directly from the simulations; no corrections are included.

0.88[89,90] 0.57 0.42 0.46 1.4[100] 2.2 2.7 2.4 33.0[99] 32.9 33.1 33.0 1 pCH2r 2 pCH2r 3 pCH2r Exp pFG6 pCG3 pCG2

773.9[89,90] 768.5 772.6 771.2

17.3[89] 17.2 18.0 18.6

24.6[89,90] 23.9 29.2 30.7

153.0[89,90] 120.6 65.4 57.0

11.5[89,90] 9.1 9.1 9.3

1.23[89,90] 1.05 0.95 0.88

2.016[99] 2.028 2.030 2.034

g (cP) D (1025 cm2 s21)  (0) aT (K21) jT (1025 atm21) Cp (J K21 mol21) c (mN m21) DFexc (kJ mol21) DHvap (kJ mol21) q (kg m23) Bead composition System

Table 2. Summary of all liquid properties of the cyclohexane models: density q, heat of vaporization DHvap, excess free energy DFexc, surface tension c, heat capacity C sim p , isothermal compressibility jT, thermal expansion coefficient aT, static dielectric permittivity (0), self-diffusion coefficient D, and shear viscosity g.

WWW.C-CHEM.ORG

Journal of Computational Chemistry 2015, 36, 1311–1321

liquid in the classical approximation. In case of the CG model, this correction term should also account for the neglect of the degrees of freedom due to approximating multiple CH2 groups by a single CG bead. For the pFG6 model, an empirical additive correction was calculated in Ref. [64]. For the pCG3 and pCG2 models, we do not calculate dCp . It is expected that the simulated heat capacity Cpsim will decrease with an increasing level of coarse-graining.

Results and Discussion Calibration of the force field parameters for CG cyclohexane The ideal bond lengths b0 between CG beads in Table 1 were obtained from mapping FG simulations to 2-bead and 3-bead CG cyclohexane. The C12 and C6 Lennard–Jones parameters (Table 1) were calibrated against the experimental density and heat of vaporization of the system (Table 2). Simulations with the initial values of these parameters obtained from FG to CG mapping did not show a good agreement with the experimental data, and a parameter search appeared necessary (Supporting Information Tables S1–S4). The C6 and C12 parameters of the 3-bead cyclohexane were obtained from scaling up r by 5% compared to the initial values from the mapped FG simulations. The value of  was not changed. The resulting repulpffiffiffiffiffiffiffi sive Lennard–Jones parameter C12 is more than three times pffiffiffiffiffi the FG equivalent, and the dispersion parameter C6 is twice as large. The CG model of cyclohexane consisting of only two beads was calibrated with smaller changes to the  and r parameters, to obtain agreement with experimental data (Supporting Information Tables S3 and S4). The values of the nonbonded Lennard–Jones parameters were obtained by scaling  down by 12 % and r up by 13 %. The resulting pffiffiffiffiffiffiffi C12 parameter is almost seven times the FG equivalent, and pffiffiffiffiffi the C6 value almost four times as large. The intermolecular Lennard–Jones potential energy functions of the calibrated model consisting of three CG beads did not change significantly in comparison to the VLJ (r) of the initial model (black dotted-dashed and solid curves in Figure 3). The minimum of the final potential energy function is only slightly deeper, and its zero value is shifted toward distances r about 0.02 nm larger than the initial value of r. For the model consisting of two beads, the change is more drastic (gray dotted-dashed and solid curves in Figure 3). The potential energy minimum shifts up by about 0.4 kJ mol21 (decreased ) and toward longer distances by nearly 0.1 nm, while the r parameter shifts toward longer distances by 0.07 nm. The structural properties of the CG models were analyzed and compared to those of the FG model in terms of the radial distribution function (Figure 4). The positions of the maxima and minima of g(r) for the CG models roughly correspond with the ones calculated from the mapped FG-to-CG distributions. The minima and maxima of the g(r) curves for both 2- and 3-bead CG models were observed at the same distances as for the FG system. With increasing bead size, the WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

FULL PAPER

Figure 3. Intermolecular Lennard–Jones potential energy VLJ (r) of liquid cyclohexane as a function of the interbead distance r for the CG cyclohexane models [eq. (5)]. Shown are the VLJ (r) functions calculated from the mapped FG simulations (dashed curves), from the CG simulations with initial Lennard–Jones parameters (dotted-dashed curves), and from the CG simulations with calibrated Lennard–Jones parameters (solid curves). Color coding: 3bead cyclohexane (pCG3)—black, 2-bead cyclohexane (pCG2)—gray.

extrema become more pronounced and agree less with the FG g(r) (black curve with filled squares in Figure 4). In all three models, FG and 2- or 3-bead CG, the first neighbor peak consists of at least two contributions. The 3-bead CG model shows even a splitting into two peaks. This may be due to two molecules of cyclohexane aligning either perpendicularly or in parallel with respect to each other. Performance of the CG cyclohexane models: liquid properties To investigate the most relevant thermodynamic, dielectric, and dynamic properties of the CG models for cyclohexane in the liquid phase, the excess free energy, the static dielectric permittivity, and the self-diffusion coefficient were calculated

Figure 4. Radial distribution functions g(r) of liquid CG cyclohexane calculated for all intermolecular distances from the mapped FG simulations (dashed curves) and the CG simulations (solid curves). g(r) of the polarizable FG (pFG6) cyclohexane is shown as black solid curve with filled squares. Color coding: 3-bead cyclohexane (pCG3)—black (see figure legend), 2-bead cyclohexane (pCG2)—gray.

Figure 5. Electronic polarization Pz (top panel) in the z-direction as a function of the applied electric field Ez ext in the simulations of liquid CG cyclohexane. The inset graph shows a close-up of the values for low field strengths, up to 0.5 e nm21, which were considered for the linear leastsquares fitting procedure to obtain the static dielectric permittivity. Shown are also the induced dipole moment lind (middle panel) and selfpolarization energy Upol (bottom panel) as a function of the applied electric field Ezext. Color coding: FG cyclohexane (pFG6)—black with filled squares, 3-bead cyclohexane (pCG3)—black with open circles, and 2-bead cyclohexane (pCG2)—gray with open circles (see figure legend).

(Table 2). The excess free energy DFexc deviates from experiment by 4 % for the 3-bead CG model and by twice that value for the case of the 2-bead model, while the FG model reproduces DFexc within 0.6% of the experimental value. This illustrates the price to be paid when coarse-graining to gain computational speed. The 3-bead model is computationally four times faster than the FG model, and the 2-bead CG model is even seven times faster. Still, all three models reproduce the experimental excess free energy within 1.3 kJ mol21. The surface tension, c, is correlated with the excess free energy and is overestimated for both CG models. The pCG3 model yields a slightly better value than the pCG2 model. The experimental value of the dielectric permittivity is reproduced well by the CG models, within 0.7% for the model consisting of three beads and within 0.9% for the model consisting of two beads. The values of the electrostatic polarization in the z-direction which were taken into account for the least-squares fit to obtain (0) are shown in the inset of the top panel of Figure 5. Values up to 0.5 e nm22 were considered for the fitting. The induced molecular dipole moment lind,mol and the self-polarization energy Upol show the expected dependence on the applied electric field (Figure 5, middle and bottom panels, respectively). Journal of Computational Chemistry 2015, 36, 1311–1321

1317

FULL PAPER

WWW.C-CHEM.ORG

solute=solvent

Table 3. Gibbs free energies of solvation DGsolv solute water or methanol in cyclohexane. Cyclohexane model Bead composition water=cyclohexane

DGsolv

pFG6 1 pCH2r

SPC=pXGn

pCG2 3 pCH2r

–

21.64 (0.30)

22.46 (0.29)

20.34 (0.32)

21.54 (0.30)

22.30 (0.30)

DGsolv

COSG2=pXGn

pCG3 2 pCH2r 21.63[101]

(exp)

DGsolv

(in kJ mol21) of a

CGW=pXGn

; mix 52.3

–

20.75 (0.11)

20.56 (0.10)

CGW=pXGn

; mix 52.0

–

20.94 (0.10)

20.74 (0.11)

DGsolv DGsolv

methanol=cyclohexane

DGsolv

26.2[102,103]

(exp)

MeOH=pXGn

–

28.01 (0.44)

29.06 (0.46)

COSM=pXGn

26.4 (0.4)[64]

26.65 (0.45)

28.13 (0.41)

DGsolv DGsolv

“pXGn” denotes a polarizable FG (X@F) or CG (X@C) cyclohexane model with n beads per molecule. “SPC” indicates that the atomistic nonpolarizable model of water[104] was used, “COSG2” indicates that the polarizable all-atom (FG) COS/G2 water model[53] was used, “CGW” indicates the polarizable CG CGW water model of Ref. [74]. Because one CGW CGW=pCGn water bead represents five water molecules, the DGsolv values obtained from the TI calculations were divided by 5. “MeOH” indicates the FG nonpolarizable methanol model from the 45A3 GROMOS force field,[79] and “COSM” indicates the FG polarizable methanol model.[56] mix indicates the values of the dielectric permittivity of the cavity within the cut-off radius for calculating the interactions between CGW and pCGn beads. Experimental values are shown in bold. Uncertainty estimates are shown in parentheses. The values for the CG models of cyclohexane were calculated with the C12 parameters of Table 1.

A dynamic property of CG cyclohexane, the self-diffusion coefficient D is overestimated by the CG models compared to the available experimental value, slightly more than by the FG model for cyclohexane. This may be due to the absence of hydrogen atoms and the small CAH dipoles in the FG and CG models based on CHn united atoms. Another dynamic property, the shear viscosity g, was also calculated and, as in the case of the FG model, the pCGn models underestimate it. Some secondary thermodynamic properties of the CG cyclohexane models were calculated (Table 2) to test the capability solute=solvent

Table 4. Gibbs free energies of hydration DGhydr solute cyclohexane in water. Cyclohexane model Bead composition cyclohexane=water

DGhydr

pFG6 1 pCH2r

pXG =SPC

pXG =COSG2 DGhydrn pXG =CGW ; mix 52.3 DGhydrn pXG =CGW ; mix 52.0 DGhydrn

pCG3 2 pCH2r

pCG2 3 pCH2r

5.15[105]

(exp)

DGhydrn

(in kJ mol21) of a

Performance of the CG cyclohexane models: solvation of polar FG molecules To investigate the performance of the CG cyclohexane models as a solvent for FG polar solutes, the Gibbs free energy of solvation in cyclohexane was calculated for a nonpolarizable and a polarizable water and methanol molecule, and for a polarizable CG water bead representing five water molecules (Table 3). An FG water molecule solvates slightly better in CG cyclohexane than a CG water bead, which can be understood from its smaller size. For the water solutes solvated in pCG3 and pCG2 cyclohexane, the experimental value is reproduced within 1 kJ mol21. We note, however, that for a CG water bead, representing five water molecules, no experimental value to compare with is available. The nonpolarizable methanol solute is slightly too well solvated, 1.8 kJ mol21 in pCG3 and 2.9 kJ mol21 in pCG2, while these values are reduced to 0.4 and 1.9 kJ mol21 for a polarizable methanol molecule. This suggest that, in particular, the pCG3 model may be used as a solvent for polar biomolecules represented at the FG united atom level of resolution. Use of the CG cyclohexane models as CG solutes

–

22.28 (0.68)

25.47 (0.67)

6.1 (2.1)[64]

21.19 (1.20)

25.98 (0.82)

–

216.17 (0.49)

216.24 (0.50)

–

216.68 (0.55)

216.42 (0.45)

The symbols are the same as in Table 3. Experimental values are shown in bold. Uncertainty estimates are shown in parentheses. The values for the CG models of cyclohexane were calculated with the C12 parameters of Table 1.

1318

of the models to mimic cyclohexane’s response to a change of thermodynamic state point. The heat capacity, Cp, is underestimated by more than 50%, which is an expected result, since Cp is generally proportional to the number of degrees of freedom in the system. The pCG3 model has 6 degrees of freedom and the pCG2 model has 5 degrees of freedom to be compared to 12 degrees of freedom of the pFG6 model. This implies that the CG model is applicable in a limited range of temperature and pressure values. We note, however, that at 298 K and 1 atm, the experimental heat of vaporization and excess free energy of cyclohexane are rather well reproduced by the CG models, which indicates a proper enthalpy/entropy ratio for the liquid. The isothermal compressibility, jT, is rather well reproduced compared to experiment for both models, with jT equal to the one of the pFG6 model for the 3-bead cyclohexane and slightly larger for the 2-bead cyclohexane. The thermal expansion coefficient, aT, is underestimated for all models, including FG cyclohexane. aT is lower but close to the pFG6 value in case of the pCG3 model, and the pCG2 model yields the lowest value. Overall, the 3-bead model of cyclohexane performs better compared to the 2-bead model for cyclohexane in the condensed phase.

Journal of Computational Chemistry 2015, 36, 1311–1321

Although the CG cyclohexane models were developed for use as solvent to solvate solutes modelled either at the FG or at the CG level of resolution, one could think of using it as a CG solute in CG water, for example, representing cyclic alkane moieties in solutes. To this end, we calculated the Gibbs free energy of solvation of CG cyclohexane in water, modelled at the FG and CG level of resolution (Table 4). We note that solvating a CG solute in an FG solvent makes computationally no sense, but we wanted to explore the limits of application of WWW.CHEMISTRYVIEWS.COM

FULL PAPER

216.68 (0.55) 215.07 (0.51) 213.49 (0.50) 4.77 (0.59) 5.10 (0.61) 5.29 (0.63) 6.63 (0.50)

5.15[105] 216.17 (0.49) 214.84 (0.38) 213.23 (0.42) 4.41 (0.52) 5.02 (0.49) 6.01 (0.49) 7.10 (0.37) 1.00 1.05 1.10 2.45 2.50 2.60 2.70

Exp

ffiffiffiffiffiffiffiffiffiffi is the scaling factor for the square root of the original C12 Lennard–Jones repulsion parameter to calculate interactions between cyclohexane and polar atom types. The C12 was scaled for the interacfpC12ð2Þ CGW=pCGn values obtained from the TI calculations were divided by 5. Valtion with both CW and DW particles of the CGW water model. Because one CGW water bead represents five water molecules, the DGsolv ues from simulation which agree with the experimental value within the statistical error are indicated with violet shading. The parameters indicated with green shading yield a best fit. The symbols are the same as in Tables 3 and 4. Experimental values are shown in bold in the first row. Values corresponding to the simulations with an unscaled C12 are shown in bold in the second row. Error estimates are shown in parentheses.

216.24 (0.50) 214.74 (0.48) 213.58 (0.42) 4.78 (0.42) 5.30 (0.54) 6.35 (0.45) 6.80 (0.54)

pCG2=CGW

DGsolv mix 52:3

pCG3=CGW

DGsolv mix 52:0

pCG3=CGW

DGsolv mix 52:3 ffiffiffiffiffiffiffiffiffiffi fpC12ð2Þ

Table 5. Gibbs free energies of solvation of CG cyclohexane in CG water DGsolv

pCGn=CGW

.

DGsolv mix 52:0

pCG2=CGW

216.42 (0.45) 214.85 (0.48) 213.70 (0.47) 4.76 (0.38) 5.73 (0.51) 5.97 (0.52) 6.89 (0.67)

DGsolv mix 52:3

CGW=pCG3

20.75 (0.11) 20.37 (0.10) 20.11 (0.12) – – – –

DGsolv mix 52:0

CGW=pCG3

20.94 (0.10) 20.62 (0.11) 20.26 (0.11) – 4.64 (0.12) – –

21.63[101]

DGsolv mix 52:3

CGW=pCG2

20.56 (0.10) 20.20 (0.11) 20.12 (0.10) – 5.12 (0.11) – –

DGsolv mix 52:0

CGW=pCG2

20.74 (0.11) 20.32 (0.13) 20.03 (0.12) – – – –

WWW.C-CHEM.ORG

the CG cyclohexane models. Tables 4 and 5 show that both types of CG cyclohexane are too well, by about 21 kJ mol2[1], solvated in CG water. Solvation of CG cyclohexane in nonpolarizable or polarizable water is also too favorable, by 7 and 6 kJ mol21 for pCG3 and by about 11 kJ mol21 for pCG2. This means that the C12 repulsive parameter (Table 1) calibrated using liquid cyclohexane properties (Table 2) is not suitable for use in combination with the (non)polarizable FG and CG water models for the solvent. pffiffiffiffiffiffiffi Variation of the C12 parameter value for the interaction between the CG cyclohexane beads and CG water beads to reproduce the experimental value of solvation of cyclohexane pffiffiffiffiffiffiffiffiffiffiffiffiffi cyclohexane=water in water (DGsolv ) resulted in a C12 ð2Þ Lennard– Jones parameter scaled up by a factor 2.5 (Table 5). Unfortunately, increasing the C12 parameter causes an increase in the CGW=pCGn Gibbs solvation free energy DGsolv of a CG water bead in CG cyclohexane, to about 4.6 kJ mol21 for pCG3 and 5.1 kJ mol21 for pCG2. So the CG cyclohexane models developed should only be used as solvent models.

Conclusion Two CG rigid models of cyclohexane were developed, a molecule consisting of three beads each representing two CH2 groups, and a molecule consisting of two beads each representing three CH2 groups. The CG bond lengths were obtained from average values of FG-to-CG mapped distributions from simulations of FG cyclohexane modelled using united (CH2) atoms. The Lennard–Jones nonbonded force field parameters were calibrated as to reproduce the experimental density and heat of vaporization of liquid cyclohexane at ambient temperature and pressure. Both models account for the electronic polarizability of cyclohexane by the addition of a COS site to each CG bead with an electronic polarizability a corresponding to the vacuum polarizability of the CH2 aliphatic groups represented by a CG bead. Structural, thermodynamic, and dielectric properties of liquid cyclohexane were reproduced by the CG models. The self-diffusion coefficient was overestimated and the shear viscosity was underestimated, as was the case for the FG model. The heat capacity was underestimated due to the inevitable loss of degrees of freedom when coarsegraining. The 3-bead model of cyclohexane performs better than the 2-bead model. The Gibbs free energy of solvation of a water molecule, nonpolarizable or polarizable at the FG level of resolution, or a polarizable CG water bead in CG cyclohexane deviates less than 1 kJ mol21 from the experimental value. For a nonpolarizable or polarizable methanol solute, the corresponding deviation is about 2 kJ mol21 for the pCG3 cyclohexane solvent and about 3 kJ mol21 for the pCG2 cyclohexane solvent. This suggests that the CG models of cyclohexane may be used as solvent models in combination with FG models, for example, of GROMOS type, for a solute.

Acknowledgments The authors would like to thank Andreas Eichenberger for his advice in regard to developing coarse-grained models for alkanes. This Journal of Computational Chemistry 2015, 36, 1311–1321

1319

FULL PAPER

WWW.C-CHEM.ORG

work was financially supported by The National Center of Competence in Research (NCCR) in Structural Biology of the Swiss National Science Foundation, and The European Research Council, which are gratefully acknowledged. Keywords: molecular dynamics  polarization  cyclohexane  coarse-graining  charge-on-spring model  GROMOS  dielectric permittivity

How to cite this article: O. M. Szklarczyk, E. Arvaniti, W. F. van Gunsteren J. Comput. Chem. 2015, 36, 1311–1321. DOI: 10.1002/jcc.23929

]

Additional Supporting Information may be found in the online version of this article.

[1] [2] [3] [4] [5]

W. F. van Gunsteren, J. Dolenc, Mol. Simul. 2012, 38, 1271. A. Rahman, Phys. Rev. 1964, 136, A405. T. X. Hoang, M. Cieplak, J. Chem. Phys. 2000, 112, 6851. V. Tozzini, Curr. Opin. Struct. Biol. 2005, 15, 144. M. Neri, C. Anselmi, M. Cascella, A. Maritan, P. Carloni, Phys. Rev. Lett. 2005, 95, 218102(1). A. Y. Shih, A. Arkhipov, P. L. Freddolino, K. Schulten, J. Phys. Chem. B 2006, 110, 3674. P. J. Bond, J. Holyoake, A. Ivetac, S. Khalid, M. S. P. Sansom, J. Struct. Biol. 2007, 157, 593. P. Derreumaux, N. Mousseau, J. Chem. Phys. 2007, 126, 025101(1). W. Han, C. K. Wan, Y. D. Wu, J. Chem. Theory Comput. 2008, 4, 1891. L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman, S.-J. Marrink, J. Chem. Theory Comput. 2008, 4, 819. R. D. Hills Jr., C. L. Brooks, III, Int. J. Mol. Sci. 2009, 10, 889. M. L. Huertas, S. Navarro, M. C. Lopez Martinez, J. Garcıa de la Torre, Biophys. J. 1997, 73, 3142. A. Morriss-Andrews, J. Rottler, S. S. Plotkin, J. Chem. Phys. 2010, 132, 035105. P. D. Dans, A. Zeida, M. R. Machado, S. Pantano, J. Chem. Theory Comput. 2010, 6, 1711. T. E. Ouldridge, A. A. Louis, J. P. K. Doye, J. Chem. Phys. 2011, 134, 085101. S. J. Marrink, A. H. de Vries, A. E. Mark, J. Phys. Chem. B 2004, 108, 750. M. Orsi, D. Y. Haubertin, W. E. Sanderson, J. W. Essex, J. Phys. Chem. B 2008, 112, 802. A. Warshel, M. Levitt, J. Mol. Biol. 1976, 103, 227. M. J. Field, P. A. Bash, M. J. Karplus, J. Comput. Chem. 1990, 11, 700. H. M. Senn, W. Thiel, Angew. Chem Int. Ed. 2009, 48, 1198. K. Meier, A. Choutko, J. Dolenc, A. P. Eichenberger, S. Riniker, W. F. van Gunsteren, Angew. Chem. Int. Ed. Engl. 2013, 52, 2820. H. J. C. Berendsen, Simulating the Physical World: Hierarchical Modeling from Quantum Mechanics to Fluid Dynamics, Cambridge University Press, Cambridge, New York, 2007. S. Riniker, A. P. Eichenberger, W. F. van Gunsteren, Eur. Biophys. J. 2012, 41, 647. T. A. Halgren, W. Damm, Curr. Opin. Struct. Biol. 2001, 11, 236. T. D. Rasmussen, P. Ren, J. W. Ponder, F. Jensen, Int. J. Quantum Chem. 2007, 107, 1390. P. J. Dyer, H. Docherty, P. T. Cummings, J. Chem. Phys. 2008, 129, 024508. P. Cieplak, F. Y. Dupradeau, Y. Duan, J. Wang, J. Phys.: Condens. Matter 2009, 21, 333102. P. E. M. Lopes, B. Roux, A. D. MacKerell, Theor. Chem. Acc. 2009, 124, 11. X. Song, J. Chem. Phys. 2002, 116, 9359. P. Soto, A. E. Mark, J. Phys. Chem. B 2002, 106, 12830. C. G. Ji, J. Z. H. Zhang, J. Phys. Chem. B 2009, 113, 16059. Y. Tong, Y. Mei, Y. L. Li, C. G. Ji, J. Z. H. Zhang, J. Am. Chem. Soc. 2010, 132, 5137.

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

[23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

1320

Journal of Computational Chemistry 2015, 36, 1311–1321

[33] Z. Lin, N. Schmid, W. F. van Gunsteren, Mol. Phys. 2011, 109, 493. [34] Y. Xiang, L. Duan, J. Z. H. Zhang, J. Chem. Phys. 2011, 134, 205101. [35] L. L. Duan, Y. Gao, Y. Mei, Q. G. Zhang, B. Tang, J. Z. H. Zhang, J. Phys. Chem. B 2012, 116, 3430. [36] P. G. Pascutti, K. C. Mundim, A. S. Ito, P. M. Bisch, J. Comput. Chem. 1999, 20, 971. [37] A. K. Rappe, W. A. Goddard, III, J. Phys. Chem. 1991, 95, 3358. [38] S. W. Rick, S. J. Stuart, B. J. Berne, J. Chem. Phys. 1994, 101, 6141. [39] S. W. Rick, S. J. Stuart, J. S. Bader, B. J. Berne, St. Phys. Theor. Chem. 1995, 83, 31. [40] D. Van Belle, I. Couplet, M. Prevost, S. J. Wodak, J. Mol. Biol. 1987, 198, 721. [41] P. Drude, The Theory of Optics, Dover Publications, New York, 1900. [42] T. P. Straatsma, J. A. McCammon, Mol. Simul. 1990, 5, 181. [43] S. A. Patel, C. L. Brooks, III, Mol. Simul. 2006, 32, 231. [44] S. W. Rick, B. J. Berne, J. Phys. Chem. B 1997, 101, 10488. [45] R. Chelli, P. Procacci, J. Chem. Phys. 2002, 117, 9175. [46] S. A. Patel, C. L. Brooks, III, J. Chem. Phys. 2006, 124, 204706. [47] J. E. Davis, G. L. Warren, S. Patel, J. Phys. Chem. B 2008, 112, 8298. [48] T. R. Lucas, B. A. Bauer, J. E. Davis, S. Patel, J. Comput. Chem. 2011, 141. [49] T. R. Lucas, B. A. Bauer, S. Patel, Biochim. Biophys. Acta 2012, 1818, 318. [50] D. van Belle, M. Froeyen, G. Lippens, S. J. Wodak, Mol. Phys. 1992, 77, 239. [51] W. Xie, J. Pu, A. D. MacKerell, Jr., J. Gao, J. Chem. Theory Comput. 2007, 3, 1878. [52] H. Yu, T. Hansson, W. F. van Gunsteren, J. Chem. Phys. 2003, 118, 221. [53] H. Yu, W. F. van Gunsteren, J. Chem. Phys. 2004, 121, 9549. [54] H. Yu, T. W. Whitfield, E. Harder, G. Lamoureux, I. Vorobyov, V. M. Anisimov, A. D. MacKerell, Jr., B. Roux, J. Chem. Theory Comput. 2010, 6, 774. [55] D. P. Geerke, W. F. van Gunsteren, Mol. Phys. 2007, 105, 1861. [56] H. Yu, D. P. Geerke, H. Liu, W. F. van Gunsteren, J. Comput. Chem. 2006, 27, 1494. [57] V. M. Anisimov, I. Vorobyov, B. Roux, A. D. MacKerell, Jr., J. Chem. Theory Comput. 2007, 3, 1927. [58] A. P. E. Kunz, A. P. Eichenberger, W. F. van Gunsteren, Mol. Phys. 2011, 109, 365. [59] Z. Lin, A. P. E. Kunz, W. F. van Gunsteren, Mol. Phys. 2010, 108, 1749. [60] I. Vorobyov, V. M. Anisimov, S. Greene, R. M. Venable, A. Moser, R. W. Pastor, A. D. MacKerell, Jr., J. Chem. Theory Comput. 2007, 3, 1120. [61] C. M. Baker, A. D. MacKerell, Jr., J. Mol. Model. 2010, 16, 567. [62] X. Zhu, A. D. MacKerell, Jr., J. Comput. Chem. 2010, 31, 2330. [63] I. Vorobyov, V. M. Anisimov, A. D. MacKerell, Jr., J. Phys. Chem. B 2005, 109, 18988. [64] O. M. Szklarczyk, S. J. Bachmann, W. F. van Gunsteren, J. Comput. Chem. 2014, 35, 789. [65] C. M. Baker, R. Best, J. Chem. Theory Comput. 2013, 9, 2826. [66] D. Robinson, J. Chem. Theory Comput. 2013, 9, 2498. [67] J. Chowdhary, E. Harder, P. E. M. Lopes, L. Huang, A. D. MacKerell, Jr., B. Roux, J. Phys. Chem. B 2013, 117, 9142. [68] V. M. Anisimov, G. Lamoureux, I. Vorobyov, N. Huang, B. Roux, A. D. MacKerell, Jr., J. Chem. Theory Comput. 2005, 1, 153. [69] C. M. Baker, V. M. Anisimov, A. D. MacKerell, Jr., J. Phys. Chem. B 2011, 115, 580. [70] H. Yu, W. F. van Gunsteren, Comput. Phys. Commun. 2005, 172, 69. [71] T. Ha-Duong, N. Basdevant, D. Borgis, Chem. Phys. Lett. 2009, 468, 79. [72] Z. Wu, Q. Cui, A. Yethiraj, J. Phys. Chem. B 2010, 114, 10524. [73] S. O. Yesylevskyy, L. V. Sch€afer, D. Sengupta, S.-J. Marrink, PLoS Comput. Biol. 2010, 6, 1. [74] S. Riniker, W. F. van Gunsteren, J. Chem. Phys. 2011, 134, 084110. [75] Mannsville, Chemical Products Synopsis, Cyclohexane, Mannsville Chemical Products Corporation, New York, 1993. [76] G. H. Peters, D. M. van Aalten, O. Edholm, S. Toxvaerd, R. Bywater, Biophys. J. 1996, 71, 2245. [77] S. Uribe, P. Rangel, G. Espınola, G. Aguirre, Appl. Environ. Microbiol. 1990, 56, 2114. [78] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, A. H. de Vries, J. Phys. Chem. B 2007, 111, 7812. [79] L. D. Schuler, X. Daura, W. F. van Gunsteren, J. Comput. Chem. 2001, 22, 1205.

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

[80] C. Oostenbrink, A. Villa, A. E. Mark, W. F. van Gunsteren, J. Comput. Chem. 2004, 25, 1656. [81] N. Schmid, A. P. Eichenberger, A. Choutko, S. Riniker, M. Winger, A. E. Mark, W. F. van Gunsteren, Eur. Biophys. J. 2011, 40, 843. [82] A. P. E. Kunz, W. F. van Gunsteren, J. Phys. Chem. A 2009, 113, 11570. [83] W. F. van Gunsteren. Available at: http://www.gromos.net. Last accessed 1 May, 2015. [84] N. Schmid, C. D. Christ, M. Christen, A. P. Eichenberger, W. F. van Gunsteren, Comput. Phys. Commun. 2011, 183, 890. [85] A. P. Eichenberger, J. R. Allison, J. Dolenc, D. P. Geerke, B. A. C. Horta, K. Meier, C. Oostenbrink, N. Schmid, D. Steiner, D. Wang, W. F. van Gunsteren, J. Chem. Theory Comput. 2011, 7, 3379. [86] A. P. Eichenberger, Molecular Dynamics Simulation of Alkanes and Proteins: Methodology, Prediction of Properties and Comparison to Experimental Data, Ph.D. thesis, ETH Z€ urich 2012. [87] A. I. Vogel, J. Chem. Soc. 1948, 1833. [88] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem. Phys. 1984, 81, 3684. [89] M. Frenkel, X. G. Hong, R. C. Wilhoit, TRC Thermodynamic Tables: Hydrocarbons, Thermodynamics Research Center, Texas Engineering Experiment Station, Texas A & M University System, Boulder, 2000. [90] M. Frenkel, Q. Dong, R. C. Wilhoit, K. R. Hall, Int. J. Thermophys. 2001, 22, 215. [91] J. P. Ryckaert, G. Ciccotti, H. J. C. Berendsen, J. Comput. Phys. 1977, 23, 327. [92] X. Daura, A. E. Mark, W. F. van Gunsteren, J. Comput. Chem. 1998, 19, 535. [93] A. P. E. Kunz, J. R. Allison, D. P. Geerke, B. A. C. Horta, P. H. H€ unenberger, S. Riniker, N. Schmid, W. F. van Gunsteren, J. Comput. Chem. 2012, 33, 340.

[94] S. Riniker, C. D. Christ, H. S. Hansen, P. H. H€ unenberger, C. Oostenbrink, D. Steiner, W. F. van Gunsteren, J. Phys. Chem. B 2011, 115, 13570. [95] S. Riniker, W. F. van Gunsteren, J. Chem. Phys. 2012, 137, 044120. [96] W. F. van Gunsteren, H. J. C. Berendsen, Groningen Molecular Simulation (GROMOS) Library Manual, Biomos, Groningen, The Netherlands, 1987. [97] J. P. M. Postma, A Molecular Dynamics Study of Water, Ph.D. thesis, Rijksuniversiteit Groningen 1985. [98] P. J. Gee, W. F. van Gunsteren, Mol. Phys. 2006, 104, 477. [99] D. R. Lide, CRC Handbook of Chemistry and Physics, CRC press, Boca Raton, Fla, 2012. [100] M. Holz, S. R. Heil, A. Sacco, Phys. Chem. Chem. Phys. 2000, 2, 4740. [101] J. Li, T. Zhu, G. D. Hawkins, P. Winget, D. A. Liotard, C. J. Cramer, D. G. Truhlar, Theor. Chem. Acc. 1999, 103, 9. [102] A. Radzicka, R. Wolfenden, Biochemistry 1988, 27, 1664. [103] J. Li, T. Zhu, G. D. Hawkins, P. Winget, D. A. Liotard, C. J. Cramer, D. G. Truhlar, Theor. Chem. Assoc. 1999, 103, 9. [104] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, J. Hermans, In Intermolecular Forces, vol. 11; Ed. Bernard Pullman; Reidel: Dordrecht, The Netherlands, 1981, pp. 331–342. [105] C. McAuliffe, J. Phys. Chem. 1966, 70, 1267.

Received: 17 December 2014 Revised: 9 March 2015 Accepted: 12 March 2015 Published online in Wiley Online Library

Journal of Computational Chemistry 2015, 36, 1311–1321

1321

Copyright of Journal of Computational Chemistry is the property of John Wiley & Sons, Inc. and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.