A GROMOS-Compatible Force Field for Small Organic Molecules in

Article pubs.acs.org/JCTC

A GROMOS-Compatible Force Field for Small Organic Molecules in the Condensed Phase: The 2016H66 Parameter Set Bruno A. C. Horta,*,†,§ Pascal T. Merz,† Patrick F. J. Fuchs,∥ Jozica Dolenc,†,‡ Sereina Riniker,† and Philippe H. Hünenberger† †

Laboratory of Physical Chemistry, ETH Zürich, CH-8093 Zürich, Switzerland Chemistry, Biology and Pharmacy Information Center, ETH Zürich, CH-8093 Zürich, Switzerland § Instituto de Química, Universidade Federal do Rio de Janeiro, Rio de Janeiro 21941-909, Brazil ∥ Institut Jacques Monod, UMR 7592 CNRS, Université Paris-Diderot, Sorbonne Paris Cité, F-75205 Paris, France ‡

S Supporting Information *

ABSTRACT: This article reports on the calibration and validation of a new GROMOScompatible parameter set 2016H66 for small organic molecules in the condensed phase. The calibration is based on 62 organic molecules spanning the chemical functions alcohol, ether, aldehyde, ketone, carboxylic acid, ester, amine, amide, thiol, sulfide, and disulfide, as well as aromatic compounds and nucleic-acid bases. For 57 organic compounds, the calibration targets are the experimental pure-liquid density ρliq and the vaporization enthalpy ΔHvap, as well as the hydration free energy ΔGwat and the solvation free energy ΔGche in cyclohexane, at atmospheric pressure and at (or close to) room temperature. The final root-mean-square deviations (RMSD) for these four quantities over the set of compounds are 32.4 kg m−3, 3.5 kJ mol−1, 4.1 kJ mol−1, and 2.1 kJ mol−1, respectively, and the corresponding average deviations (AVED) are 1.0 kg m−3, 0.2 kJ mol−1, 2.6 kJ mol−1, and 1.0 kJ mol−1, respectively. For the five nucleic-acid bases, the parametrization is performed by transferring the final 2016H66 parameters from analogous organic compounds followed by a slight readjustment of the charges to reproduce the experimental water-to-chloroform transfer free energies ΔGtrn. The final RMSD for this quantity over the five bases is 1.7 kJ mol−1, and the corresponding AVED is 0.8 kJ mol−1. As an initial validation of the 2016H66 set, seven additional thermodynamic, transport, and dielectric properties are calculated for the 57 organic compounds in the liquid phase. The agreement with experiment in terms of these additional properties is found to be reasonable, with significant deviations typically affecting either a specific chemical function or a specific molecule. This suggests that in most cases, a classical force-field description along with a careful parametrization against ρliq, ΔHvap, ΔGwat, and ΔGche results in a model that appropriately describes the liquid in terms of a wide spectrum of its physical properties. Condensed-phase force fields such as CHARMM,5−15 AMBER,16−31 OPLS,32−44 TraPPE,45−52 and GROMOS53−68 mainly aim at the description of solids, liquids, and solutions, including solvated biomolecules as a particularly relevant special case. They usually involve a large number of atom types accounting for atoms in different chemical environments, a relatively simple functional form, empirical parametrization procedures, the extensive use of parameter combination and transferability assumptions, and parameters calibrated mainly against experimental thermodynamic and spectroscopic data concerning liquids and molecules in solution. Because their principal aim is an accurate description of condensed-phase properties including solvation effects and conformational equilibria, the emphasis during parametrization is placed on the proper description of nonbonded and torsional-angle interactions. The GROMOS force field is one of the widely used force fields of this type. The main principles underlying its construction can be summarized as follows:55,67,69 (1) United-atom representation70 of aliphatic CH, CH2, and CH3 groups.56,59,60

1. INTRODUCTION Molecular dynamics (MD) simulation represents a powerful tool for investigating the properties of molecular systems relevant in physics, chemistry, and biology.1−4 The usefulness of this method in the context of condensed-phase systems results in particular from a favorable trade-off between model resolution and computational cost. Although classical atomistic models represent an approximation to quantum-mechanical (QM) ones, they can still provide a realistic description of molecular systems at spatial and temporal resolutions on the order of 0.1 nm and 1 fs, respectively, whereas their computational cost remains tractable at present for system sizes and timescales on the order of 10 nm and 1 μs, respectively. These scales are sufficient to enable in many cases: (1) a proper description of solvation by explicit treatment of the solvent molecules within a sufficiently large solvation volume; (2) the calculation of converged thermodynamic properties via statistical mechanics; (3) a direct comparison with experimental data, namely structural, thermodynamic, transport, and dynamic observables measured on similar spatial and temporal scales. However, the ability of MD simulations to represent accurately the properties of a given molecular system depends crucially on the quality of the underlying force field. © 2016 American Chemical Society

Received: February 20, 2016 Published: June 1, 2016 3825

DOI: 10.1021/acs.jctc.6b00187 J. Chem. Theory Comput. 2016, 12, 3825−3850

Article

Journal of Chemical Theory and Computation

The set labeled68 54A8 is based on 54A7 and differs in the atomic partial charges within charged amino acid side-chains. The set labeled91,92 53A6GLYC is based on 53A6 and includes modified torsional-angle parameters for hexopyranoses. The set labeled93 56A6CARBO represents a full reparametrization of hexopyranose-based carbohydrates. It is also based on 53A6 and includes three additional atom types for ring atoms. A revised version94 56A6CARBO_R was also proposed recently to improve the representation of ring conformational properties in carbohydrate chains. Although the five variants mentioned above represented important changes relative to 53A6 in the context of specific biomolecules, none of them significantly altered the basis of the force field as defined by the condensed-phase (pure-liquid and hydration) properties of small organic molecules. In 2011, however, we introduced a parameter set labeled77 53A6OXY with the aim of reproducing simultaneously the experimental data concerning the pure-liquid as well as the aqueous and nonpolar solvation properties of small organic molecules. The calibration of this set was based on: (1) a simultaneous adjustment of the Lennard-Jones interaction parameters and of the atomic partial charges (vs separate adjustments in the design63 of 53A6); (2) a slight adjustment of the dispersive C1/2 6 Lennard-Jones interaction coefficients along with that of the corresponding repulsive C1/2 12 coefficients (vs an adjustment of 63 of 53A6). This additional flexibility C1/2 12 only in the design permitted to resolve77 the apparent incompatibility between the reproduction of pure-liquid and hydration properties that had led to the splitting63 between 53A5 and 53A6. Only oxygen compounds were considered in this reoptimization, namely, alcohols, ethers, aldehydes, ketones, carboxylic acids, and esters. The refinement involved changes in the Lennard-Jones parameters for interactions involving oxygen atoms and in the atomic partial charges within the above chemical functions, along with marginal adjustments in the reference bond lengths and angles within these functions. The 53A6OXY parameter set reproduced the target experimental properties of 35 organic compounds at atmospheric pressure and at (or close to) room temperature, with root-mean-square deviations of 22.4 kg m−3 for the liquid densities ρliq, 3.1 kJ mol−1 for the vaporization enthalpies ΔHvap, 3.0 kJ mol−1 for the hydration free energies ΔGwat, and 1.7 kJ mol−1 for the solvation free energies ΔGche in cyclohexane. The 53A6OXY set was extended to a 53A6OXY+A set95 including the N-monosubstituted amide function and to a 53A6OXY+D set96 including the vicinal diether function, and validated in different contexts.97−114 The goal of the present work is to extend the parametrization strategy applied in the design of 53A6OXY, 53A6OXY+A, and 53A6OXY+D to other organic molecules (Table 1), namely, amines, amides (unsubstituted and N-disubstituted), thiols, sulfides, disulfides, aromatic compounds, and nucleic-acid bases. The resulting parameter set is labeled 2016H66, where 2016 stands for the year, H corresponds to the initials of the first and last authors, and 66 denotes the number of atom types in the set. It encompasses 53A6OXY as well as the additions of 53A6OXY+A and the modifications of 53A6OXY+D, i.e., it relies on the same force-field parameters for the compounds considered in the corresponding articles.77,95,96 This is true with a small 1/2 coefficient for the thirdproviso: the Lennard-Jones C12 neighbor interactions of an aliphatic CH2 group within a ring was reduced to the value adopted in93 56A6CARBO, which affects in particular cyclohexane and thus slightly alters the calculated ΔGche values. Taken together, the calibration set of 2016H66 includes 62 small organic molecules: 35 oxygen compounds from 53A6OXY

(2) Covalent force-field terms including quartic (or harmonic, or constrained) bond stretching, cosine-harmonic (or angle-harmonic) bond-angle bending, harmonic improper dihedral-angle distortion, and cosine-series torsional dihedral-angle rotation. (3) Lennard-Jones representation71 of the van der Waals interactions. (4) Nonbonded exclusion of first and second (in some cases, third) covalent neighbors. (5) Lennard-Jones interaction reduction for third covalent neighbors using a special set of parameters, along with unaltered electrostatic interactions. (6) Application of a geometric-mean combination rule72,73 for pairwise Lennard-Jones interaction parameters, distinguishing between non-hydrogen-bonding, uncharged hydrogen-bonding, and oppositely charged interactions. (7) Freely adjustable atomic partial charges, i.e., the charge is not determined by the Lennard-Jones atom type. (8) Compatibility with the simple point charge (SPC) water model.74 (9) Compatibility with reaction-field electrostatics75,76 and Lennard-Jones interactions excluding a long-range correction based on an effective cutoff distance of 1.4 nm. Note that (4)−(6) represent default rules, which can be overridden in exceptional cases. A complete description of the functional form of the GROMOS force field can be found in refs 55, 57, 58, 63 (see also the GROMOS Web site69). The pure-liquid and solvation properties of small organic molecules have played a central role in the parametrization of the recent versions of the GROMOS force field.56,60,63,77 Along the successive revisions, the calibration has been performed in the first place against primary experimental data (observed thermodynamic and spectroscopic properties) considering small molecules. An assumption of transferability was then invoked to justify the use of parameters calibrated in this context for chemically similar fragments within larger molecules (e.g., biomolecules). The validity of this assumption was tested explicitly in numerous simulations of complex molecular systems, including solvated biomolecules.67,78−83 The most recent major reparametrization of the GROMOS force field63 has led to two distinct parameter sets labeled 53A5 and 53A6, exclusively differing in the atomic partial charges within polar functional groups. Set 53A5, which was calibrated to reproduce the densities and vaporization enthalpies of pure liquids, was recommended for the simulation of organic liquids and liquid mixtures. Set 53A6, in which the atomic partial charges were readjusted so as to reproduce the hydration free energies of polar amino acid side-chain analogs, was recommended for the simulations of biomolecular systems. The decision of providing two distinct parameter sets resulted from the apparent impossibility of reproducing pure-liquid and hydration properties simultaneously with a sufficient accuracy. The parameter set included in the latest release of the GROMOS program69,84−88 in 2011 is labeled67 54A7. This set is based on 53A6 and encompasses limited changes affecting specific classes of biomolecules (modified Lennard-Jones selection rule and torsional-angle parameters involved in the backbone of peptides, new atom type in the choline moiety of phospholipids89) as well as new sodium and chloride ion parameters (taken from set L in ref 90). Recently, a number of other modifications or extensions of 53A6 have been proposed. 3826


Article


Table 1. Organic Compounds Considered in the Calibration and Validation of 2016H66 along with Experimental Values of Their Target Propertiesa chemical function alcohols

ethers

aldehydes

ketones

carboxylic acids

esters

amines

amides

thiols

sulfides/disulfides

aromatic compounds

code

name

MTL ETL 1PL BTL PTL HXL HPL OTL 2PL 2BTL 2PTL 3PTL CHXL 2M2P 2M2B DME DEE MPH DXE EAL PAL BAL PPN BTN 2PN 3PN 2HN 3HN ACA PPA BTA EAE MPE PAE BAE EAN PAN BAN DEAN TEAN EDAN AMD LNMA DAMD ETSH PRSH BTSH DMS DES EMS DMDS BZN TOL PHL PHT ANLN PYRI

methanol ethanol propanol butanol pentanol hexanol heptanol octanol propan-2-ol butan-2-ol pentan-2-ol pentan-3-ol cyclohexanol 2-methylpropan-2-ol 2-methylbutan-2-ol methoxymethane ethoxyethane 1-methoxypropane 1,2-dimethoxyethane acetaldehyde propionaldehyde butyraldehyde propanone butanone pentan-2-one pentan-3-one hexan-2-one hexan-3-one acetic acid propionic acid butyric acid ethylacetate methylpropionate propylacetate butylacetate ethylamine propylamine butylamine diethylamine triethylamine ethylenediamine acetamide N-methylacetamide dimethylacetamide ethanethiol propanethiol butanethiol dimethylsulfide diethylsulfide ethylmethylsulfide dimethyldisulfide benzene toluene phenol p-hydroxytoluene aniline pyridine

calibration × ×

× × × ×

× ×

×

× × ×

× × × × × ×

× × × × × × × × ×

ε

Tl+g (K)

33.50b 24.00b 20.00b 17.70b 15.10b 13.00b 11.50b 10.10b 19.10b 16.70b 13.80b 13.40b 16.40b 11.50b 5.70b 6.20b 4.20b 3.70j 7.30b 21.10b 18.40b 13.40b 20.80b 17.70b 15.40b 16.60b 14.50b 10.57j 6.20b 3.40b 2.90b 6.00b 6.00b 5.60b 5.10b 8.70b,n 5.10b,o 4.90b,k 3.90b,k 2.40b 13.80b,k 65.00b,q 139.50b 37.80b,o 6.70b 5.70b 5.10b 6.70b,w 5.70b 5.87j 9.60b 2.30b 2.40b 13.40b,k 13.50b 6.70b,o 12.90b

254

364 333

333 309

ρliq (kg m−3)

ΔHvap (kJ mol−1)

ΔGwat (kJ mol−1)

ΔGche (kJ mol−1)

784.00c 784.93c 799.60c 805.75c 810.80c 815.34c 820.00e 821.57c 780.00c 802.41c 805.40c 816.00c 968.40c 781.20c 805.00c 722.00h 707.82c 735.60f 863.70c 778.00c,k 791.20c 796.40c 784.40c 799.70c 801.50c 809.45c 806.70c 810.00e 1043.92c 988.08c 953.20c 894.55c 909.00e 883.03c 876.36c 677.00e 712.10c 736.83c 701.60c 723.05c 893.10c 989.20c 926.00t 936.34c 833.00e 841.10f 836.74c 842.28c 831.18c 832.00e 1062.50c,k 873.60c 862.19c 1054.46c,x 1018.50e,z 1017.50c 978.24c

37.43c 42.31c 47.49c 52.34c 56.94c 61.85c 66.81f 70.98c 45.52c 49.66c 53.10c 53.10c 62.01c 46.82c 50.20c 21.24f,i 27.10c 27.60f 36.39c 26.11c 29.63c 33.68c 31.30c 34.51c 38.40c 38.56c 42.90c 42.47f 51.60l 57.30l 58.00l 35.62c 35.85m 39.83c 43.64c 26.60l 31.26c 35.74c 31.32c 34.81c 46.00l 56.10c,r 54.39u,v 49.15c 27.30l 31.89f 36.53c 27.65c 35.77c 31.85f 38.33c 33.84c 37.99c 55.40y 64.96y 55.83c 40.41c

−21.37d −20.95d −20.58d −19.75d −18.70d −18.24d −17.70d −17.10g −19.33d −19.16g −18.36g −18.20g −22.90g −18.87d −18.53g −8.03d −7.36d −6.95d −20.25d −14.60g −14.40g −13.30g −16.10d −15.73d −14.77d −14.26d −13.76d

−5.40d −10.84d −11.42d −14.73d −15.10d −22.22d −25.18d

−28.30d −27.07d −26.60d −13.05d −12.27d −11.96d −10.66d −18.83d −18.37d −17.95d −17.03d −13.47p

−7.24d −15.80d

−40.63d −41.84d −35.56g −5.18d −4.39d −4.14p −6.44d −5.98d −6.28p −7.66d −3.64d −3.72d −27.70d −25.69d −22.97d −19.66d

−12.60s

−9.92d

−12.26d

−12.68d

−11.17d −14.56d −17.53d −18.00d −19.96d

−14.89d −15.53d −18.23d −20.66d −8.54d −15.56d −15.10d

−13.05d

−15.80s −17.53d −20.50d −23.30d −24.64d −23.10d −17.99d

a

The abbreviations (code) used in this article and the official IUPAC names115 (name) are provided along with the experimental values of the liquid static relative dielectric permittivity (ε), liquid density (ρliq), vaporization enthalpy (ΔHvap), hydration free energy (ΔGwat) and solvation free energy 3827


Article

Journal of Chemical Theory and Computation Table 1. continued

in cyclohexane (ΔGche). These quantities are standard, referring to 1 bar and, unless otherwise specified, 298.15 K. For the solvation free energies, the standard-state definition involves identical reference concentrations for the gas-phase and solute standard states. The temperature considered for the calculation of ρliq and ΔHvap is also indicated (Tl+g) when differing from 298.15 K. The calculation of ΔGwat and ΔGche was always performed at 298.15 K. Compounds marked with a cross (calibration) were used in the calibration. Only 57 of the 62 compounds considered are listed here, the five nucleotide bases being listed separately in Table 4. A version of this table including additional information as well as a list including chemical structures can be found in SI B (Tables B1 and B2). bRef 116. cRef 117. dRefs 118−120. eRef 121. fRef 122. gRef 123. hRef 124. i Interpolated between values at 243 and 263 K. jRef 125. k293 K. lRef 126. mRef 127. n273 K. o303 K. pRef 128. q384 K. r494 K. sRef 129. tRef 130. u Ref 131. v410 K. w294 K. x318 K. yref 132. z313 K.

(or 53A6OXY+D for the diether dimethoxyethane), 1 N-monosubstituted amide from 53A6OXY+A (N-methylacetamide), 21 compounds of the newly considered chemical functions, and 5 nucleic-acid bases. For the 21 new organic compounds, the nonbonded interaction parameters were systematically and simultaneously reoptimized against experimental data for ρliq, ΔHvap, ΔGwat, and ΔGche. For the five nucleic-acid bases, the Lennard-Jones interaction parameters were transferred by analogy from the reoptimized organic compounds, using pyridine as a template for the aromatic ring atoms. The charges were transferred analogously but subsequently readjusted to reproduce the experimental water-to-chloroform transfer free energies ΔGtrn of the bases. Based on the final 2016H66 parameters, agreement with experimental data in the context of the 57 organic molecules considered (excluding the nucleic-acid bases) is checked not only in terms of the calibration targets ρliq, ΔHvap, ΔGwat, and ΔGche but also in terms of 7 other pure-liquid thermodynamic, transport, and dielectric properties, namely the molar isobaric the heat capacity cP, the isothermal compressibility κT, the isobaric thermal-expansion coefficient αP, the surface tension coefficient γ, the self-diffusion coefficient D, the shear viscosity η and the static relative dielectric permittivity ε. For the sake of conciseness, this article does not detail the optimization process that has led to the 2016H66 parameter set, i.e., only the final set is described. For the same reason, the main article only discusses the validation of the set in relatively broad terms, mainly considering overall deviations for compounds representative of a specific chemical function. However, all the details can be found in two Supporting Information documents A (SI A) and B (SI B). SI A provides a specification of the 2016H66 parameter set. This information corresponds primarily to the content of the GROMOS interaction function parameter (ifp) file. It must be complemented by that contained in the GROMOS molecular topology building block (mtb) files, providing the topology of the 62 organic molecules considered plus cyclohexane. The 2016H66 files for GROMOS (ifp and mtb), as well as corresponding files for GROMACS (itp and rtp), can be downloaded free of charge from the Web site of the last author.134 These files also encompass the GROMOS solvents74,135−142 (see Table A10), alkali and halide ions,90 and 2016H66compatible biomolecular building blocks for peptides, nucleic acids, carbohydrates (including cyclodextrins143), lipids (including mono- and diglycerides 97,101,107−110,144), and cofactors (see details in the Appendix of SI A; note that the GROMACS carbohydrate files built using pdb2gmx must be postprocessed using a script provided in the distribution). It should be stressed, however, that these biomolecular building blocks have undergone only preliminary tests and not yet a thorough validation. SI B provides more information on the organic molecules considered, including the experimental values of their properties

as well as the detailed simulation results concerning these properties.

2. COMPUTATIONAL DETAILS 2.1. Force-Field Representation. 2.1.1. Covalent Interaction Parameters. Except for the few changes made in77 53A6OXY (change of assignment for some reference bond lengths and angles within oxygen-containing chemical functions) and in96 53A6OXY+D (addition of six new torsional-angle types for vicinal diether functions), the covalent interaction parameters of 2016H66 relevant for the 62 organic molecules considered as well as cyclohexane are the same as those one would select for the corresponding molecules in the 53A6 force field.63 The list of bond, bond-angle, torsional-dihedral, and improper-dihedral types of 2016H66 can be found in SI A (Tables A1−A4). 2.1.2. Atom Types and Combination Rules. In GROMOS, the van der Waals interactions are calculated according to a Lennard-Jones function71 VLJ(r) =

⎛C

∑ ∑′ ⎜⎜ i

j>i

⎝

12, ij rij12

−

C6, ij ⎞ ⎟ rij6 ⎟⎠

(1)

where r is the Cartesian coordinate vector of the entire system, rij is the (minimum-image) distance between two interaction sites i and j in configuration r, and the prime at the second summation sign indicates that excluded neighbors are omitted. Normally, the parameters C6,ij and C12,ij are determined by the atom types (integer atom code or IAC) I(i) and J(j) of atoms i and j, respectively, according to a geometric-mean combination rule72,73 1/2 1/2 C12, ij = C12, m(I ) · C12, n(J )

and

C6, ij = C61/2(I ) ·C61/2(J ) (2)

C1/2 12,m

For each atom type, up to three different parameters are 1/2 1/2 defined (labeled C1/2 12,I, C12,II, and C12,III), which are associated with increasing repulsiveness. They are to be used in the cases of nonhydrogen-bonding, uncharged hydrogen-bonding, or oppositely charged (C1/2 12,III used for the anion only) interactions, respectively. The selection in eq 2 is made according to a combination matrix M in the form m = MIJ

and

n = MJI

(3)

There are two types of exceptions to the application of eqs 2 and 3: (1) third covalent neighbors rely on a special set of parameters, 1/2 labeled C1/2 6,nei and C12,nei; (2) exceptions can be explicitly included in the topology file either for specific atom-type pairs (currently used only for the GROMOS chloroform model136) or for specific atom pairs in a molecule (currently used only for carbohydrates in the 56A6CARBO and 56A6CARBO_R parameter sets93,94). The 2016H66 set encompasses 66 atom types (listed in Table 2). The first 53 types have the same name and application range as in set63 53A6, although the Lennard-Jones interaction parameters associated with a number of these types have 3828


Article

Journal of Chemical Theory and Computation Table 2. Atom Types of 2016H66a IAC

type

usage

1 2

O OM

carbonyl oxygen oxygen bearing a negative charge (e.g., carboxylate) oxygen bound to hydrogen (e.g., alcohol, carboxylic acid) ether or ester oxygen SPC or SPC/E water oxygen amide nitrogen amine nitrogen

3

OA

4 5 6 7

OE OW N NT

8

NL

9

NR

10 11 12

NZ NE C

13 14 15 16 17 18

CH0 CH1 CH2 CH3 CH4 CH2r

19 20 21

CR1 HC H

22 23

DUM S

primary guanidinium nitrogen secondary guanidinium nitrogen sp2 bare carbon (e.g., carbonyl, isobutylene, benzene) nonring sp3 bare carbon (e.g., neopentane) aliphatic nonring united-atom CH group aliphatic nonring united-atom CH2 group aliphatic nonring united-atom CH3 group united-atom methane molecule aliphatic united-atom CH2 group in ring (e.g., cyclohexane, hexopyranoses) united-atom aromatic CH group hydrogen bound to carbon hydrogen bound to element other than carbon dummy atom sulfur

24 25 26 27 28 29 30 31

CU1+ CU2+ FE ZN2+ MG2+ CA2+ P,SI AR

copper (charge +1) copper (charge +2) iron (heme group) zinc (charge +2) magnesium (charge +2) calcium (charge +2) phosphor or silicon argon

nitrogen bearing a positive charge (e.g., ammonium) aromatic nitrogen

parameters from

IAC

53A6OXY 53A6 53A6OXY 53A6OXY 53A6 53A6OXY+A present work 53A6 present work 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 56A6CARBO 53A6 53A6 53A6 53A6 present work 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6

type

32 33 34 35

F CL BR CMet

36 37 38 39 40 41 42 43

OMet NA+ CL− CChl CLChl HChl SDmso CDmso

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

ODmso CCl4 CLCl4 FTFE CTFE CHTFE OTFE CUrea OUrea NUrea CH3p LLI+ LNA+ LK+ LRB+ LCS+ LF− LCL− LBR− LI− Or CH0r

66

CH1r

usage fluorine (nonionic) chlorine (nonionic) bromine (nonionic) united-atom methyl group in methanol (solvent) oxygen in methanol (solvent) sodium (charge +1) chloride (charge −1) carbon in chloroform (solvent) chloride in chloroform (solvent) hydrogen in chloroform (solvent) sulfur in DMSO (solvent) united-atom methyl group in DMSO (solvent) oxygen in DMSO (solvent) carbon in carbontetrachloride (solvent) chloride in carbontetrachloride (solvent) fluor in trifluoroethanol carbon in trifluoroethanol united-atom CH2 group in trifluoroethanol oxygen in trifluoroethanol carbon in urea oxygen in urea nitrogen in urea united-atom methyl group in choline lithium (charge +1) sodium (charge +1) potassium (charge +1) rubidium (charge +1) cesium (charge +1) fluoride (charge −1) chloride (charge −1) bromide (charge −1) iodide (charge −1) oxygen in ring (e.g., hexopyranoses) sp3 bare carbon in ring (e.g., hexopyranoses) aliphatic united-atom CH1 group in ring (e.g., hexopyranoses)

parameters from 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 53A6 54A7 ref 90, set L ref 90, set L ref 90, set L ref 90, set L ref 90, set L ref 90, set L ref 90, set L ref 90, set L ref 90, set L 56A6CARBO 56A6CARBO 56A6CARBO

a

For the 66 types, the integer atom code (IAC) and corresponding atom-type name (type) are provided along with the typical usage (usage) and the 1/2 origin of the associated Lennard-Jones interaction parameters. The corresponding single-atom Lennard-Jones interaction parameters C1/2 6 and C12 entering into eq 2, the mass-type codes, and the combination matrix of eq 3 are provided in SI A (Tables A5−A7; see also Appendix therein).

It was applied here to 21 other compounds involving different chemical functions. A list of the 57 organic molecules considered, along with associated three- or four-letter codes used to refer to them in this article, is provided in Table 1. Additional information can be found in SI B in the form of a more detailed table including chemical abstract service (CAS) registry numbers and chemical formulae as well as melting (Tm) and boiling (Tb) temperatures (Table B1) and a corresponding list including drawings of the chemical structures (Table B2). The thermodynamic observables against which the parametrization was performed are: (1) the density ρliq of the pure liquid; (2) the vaporization enthalpy ΔHvap of the pure liquid; (3) the hydration free energy ΔGwat of the compound; (4) the solvation free energy ΔGche of the compound in cyclohexane. The experimental standard values for these four quantities at 1 bar and, unless otherwise specified, 298.15 K are reported in Table 1 based on refs 116−133, 145. For the solvation free energies, the standard-state definition involves identical reference

been reoptimized. The last 13 types are irrelevant for the organic molecules considered here, and included only for the extension of 2016H66 to biomolecules (see Appendix of SI A for more information). The Lennard-Jones interaction parameters associated with the 66 atom types of 2016H66 as well as the corresponding masses are provided in SI A (Tables A5 and A6). The 2016H66 parameter set uses the same combination matrix M as 53A6 for their 53 common atom types. Extension to the 13 additional types leads to the matrix provided in SI A (Table A7). Note that the change of selection rule that was introduced67 from 53A6 to 54A7 for the interaction of type O with type N, namely the use of C1/2 12,I instead of 1/2 C1/2 12,II for type O (along with C12,II for type N) has not been retained. 2.2. Parametrization Strategy. 2.2.1. Organic Liquids. The calibration strategy for the nonbonded interaction parameters was essentially the same as that used in77 53A6OXY for 35 oxygen compounds and in95 53A6OXY+A for N-methylacetamide. 3829


Article

Journal of Chemical Theory and Computation concentrations (e.g., 1 mol dm−3) for the gas-phase and solute standard states (rather than a reference pressure of 1 bar in the gas phase). Exceptions regarding the temperature concern five compounds that are not liquid under ambient conditions (see Tm and Tb in Table B1), and for which a slightly different temperature was selected for the pure-liquid and gas-phase simulations required in the determination of ρliq and ΔHvap (see Tl+g in Table 1). The parameter calibration focused primarily on 26 compounds (indicated with a cross in the table), whereas the 31 others were used to test the transferability of the reoptimized parameters. In this sense, the latter compounds belong to a validation set rather than to the calibration set. Considering the 57 molecules, experimental values were available either for all compounds (ρliq and ΔHvap), for all but two compounds (ΔGwat; no value for 3HN and EDAN), or for 33 compounds (ΔGche). For the 26 calibration molecules, only a few experimental values were missing (10 values for ΔGche). The Lennard-Jones interaction parameters of the aliphatic bare carbon, of the united-atom aliphatic groups, and of methane (atom types CH0, CH1, CH2, CH2r, CH3, and CH4 in Table 2) were not altered relative to those of 53A6, with one exception (CH2r; see below). These parameters have been previously optimized56,60 in the context of linear, branched, and cyclic alkanes against experimental data for ρliq, ΔHvap, and ΔGwat. They were subsequently shown to reproduce very well other experimental properties, including solvation free energies in different polar81 and nonpolar78,81 solvents, as well as conformational properties of hydrocarbons.59 These parameters have also been validated in the context of numerous biomolecular simulations and especially of lipids.89,97,101,107−110,144 The only change performed in 2016H66 concerns the ring methylene atom type CH2r, for which the third-neighbor interaction parameter C1/2 12,nei was reduced to the value adopted in the 56A6CARBO parameter set,93 as this modification was shown to improve the conformational properties of cyclohexane without altering significantly its pure-liquid and solvation properties. Additional aliphatic ring atoms types (CH0r and CH1r) were also imported from 56A6CARBO to handle substituted cyclohexanes. In the 53A6OXY parameter set,77 the calibration procedure involved the refinement of the Lennard-Jones interaction para1/2 1/2 meters C1/2 6 , C12,I and C12,II associated with the atom types O (carbonyl oxygen), OA (alcohol or carboxylic acid oxygen), and OE (ether or ester oxygen), as well as of the partial charges of the atoms involved in all oxygen-containing chemical functions. This led to distinct charge sets for the alcohol, ether, aldehyde, ketone, carboxylic acid, and ester functional groups. The adjustment of the C1/2 6 parameters was kept minimal, considering the qualitative connection that exists between this parameter and the electronic polarizability of the corresponding atom.146 In the 53A6OXY+A parameter set,95 the same procedure was used to refine the Lennard-Jones interaction parameters associated with atom type N (amide nitrogen), as well as the partial charges for the N-monosubstituted amide group. In the present work, the same procedure was applied to calibrate the Lennard-Jones interaction parameters associated with the atom types NT (amine nitrogen), NR (aromatic nitrogen), and S (thiol, sulfide, and disulfide sulfur) and the partial charges of the functional groups amine, unsubstituted amide, N-disubstituted amide, thiol, sulfide, and disulfide, as well as of the aromatic compounds. Additional features of the calibration were that: (1) the reparametrization did not involve any change in the third-neighbor Lennard-Jones interaction parameters of the concerned atom

types; (2) the repulsion parameters C1/2 12,III corresponding to oppositely charged interactions are irrelevant (undefined) for these atom types because they always belong to neutral functional groups; (3) the 2016H66 parameters (LennardJones and charge sets) for oxygen-containing functions and N-monosubstituted amides are identical to those of 53A6OXY and 53A6OXY+A, respectively.77,95 The nonbonded interaction parameters of 2016H66 are provided in SI A, Table A5 (Lennard-Jones; see also the atom types in Table 2, the mass types in Table A6, and the combination matrix in Table A7) and in Table 3 (atomic partial charges in the functional groups relevant for the 57 organic molecules). The details concerning the simulations required for the evaluation of ρliq, ΔHvap, ΔGwat, and ΔGche are provided in sections 2.3.1−2.3.6 and those for the more extensive analysis of pure-liquid properties are provided in section 2.4. 2.2.2. Nucleic-Acid Bases. The five methylated nucleicacid bases relevant for DNA and RNA are 9-methyladenine, 9-methylguanine, 1-methylcytosine, 1-methylthymine, and 1-methyluracil (Table 4). The corresponding Lennard-Jones interaction parameters were transferred directly from the previously reoptimized aromatic compounds benzene, toluene, and pyridine (atom types HC, C, NR, and CH3) along with the functional groups amine and ketone (atom types H, NT, and O) via the choice of identical atom types. For the optimization of the partial charges, the water-tochloroform transfer free energies ΔGtrn = ΔGchl − ΔGwat, where ΔGchl is the solvation free energy in chloroform, were taken as target. The corresponding experimental values for the five methylated bases were derived from the measured147 chloroform− water partition coefficients P (concentration in chloroform divided by concentration in water) at 1 bar and 298.15 K as ΔGtrn = −2.303RT log P, where R is the ideal-gas constant. The experimental values of P and ΔGtrn can be found in SI B (Table B7). The influence of the partial-charge set on the calculated ΔGtrn value was explored for each base using the thermodynamic cycle illustrated in Figure 1. First, ΔGtrn was calculated accurately using thermodynamic integration148 (TI) considering the original 53A6 charge set64 along with the 2016H66 Lennard-Jones interaction parameters. This was done by decoupling the solute−solvent interactions in two distinct calculations, one in water and one in chloroform. Second, two long plain MD simulations were performed using the same force-field combination for the base, one in chloroform and one in water. Based on these two simulations, the effect of various charge changes on the relative free energy in the two environments was estimated using the one-step perturbation (OSP) method.149,150 Combined with the above TI value, these perturbation estimates permitted a prediction of ΔGtrn values for arbitrary charge sets. The charge adjustments leading to a predicted ΔGtrn value in agreement with experiment, while remaining as limited as possible (and consistent across bases as well as compatible with the charge-group definitions of GROMOS), were retained to define the charge set of 2016H66. The final atom-type assignments and charge sets for the five bases are reported in Table 4. For these final parameters, the values of ΔGtrn were also recalculated accurately using TI. For comparison purposes, they were calculated as well considering the parameters (Lennard-Jones and charge set) of the original 53A6 force field.63 2.3. Simulation Protocols. 2.3.1. Simulation Parameters. All MD and stochastic dynamics151−155 (SD) simulations were performed using the latest version of the GROMOS package of programs.69,84−88 The parameter sets considered are the 2016H66 3830


Article


Table 3. Atom-Type Assignments and Atomic Partial Charges of the Different Organic Chemical Functions in 2016H66a

a

For each functional group (group) and for each atom within the group (or set of equivalent atoms for the aromatic compounds), the integer atom code (IAC), atom-type name (type), and atomic partial charge (q) are reported. This table includes only the chemical functions relevant for the 57 organic compounds of Table 1 and is repeated in SI A (Table A8). The corresponding information for the nucleic-acid bases can be found in Table 4. A complete topology specification (including charge-group definitions) for the 57 organic molecules in 2016H66 can be found in the mtb files.134

set in its final form and, for comparison purposes, the 53A5 and 53A6 sets.63 The following choices were made in the simulations employing 53A5 or 53A6: (1) the ethers have identical parameters in 53A5 and 53A6 (see, however, ref 66 for a 53A6 variant); (2) the aldehydes are available neither in 53A5 nor in 53A6; (3) the ketones have identical parameters in 53A5 and 53A6; (4) the charge set retained for the esters in 53A6 is the one of ref 62 (entry q1 of Table 2 therein); (5) the charge set retained for the amines in 53A6 is the one of ref 156 (Table 1 therein); (6) the amide DAMD and the aromatic compound

PYRI are available neither in 53A5 nor in 53A6. For the condensed-phase simulations, the solvent models employed are the SPC water model,74 the 2016H66 cyclohexane model (53A6 model60,63 with the modification93 of C1/2 12,nei for atom type CH2r), and the 53A6 chloroform model.63,136 The details of the SD simulations in the gas phase are provided in section 2.3.3, and the protocols described in this section pertain only to the condensed-phase MD simulations. All condensed-phase MD simulations were carried out under periodic boundary conditions based on cubic computational boxes. 3831


Article

Journal of Chemical Theory and Computation Table 4. Atom-Type Assignments and Atomic Partial Charges of the Nucleic-Acid Bases in 2016H66a

a

For each methylated base (base) and for each atom within the base, the integer atom code (IAC), atom-type name (type), and atomic partial charge (q) are reported. This table is repeated in SI A (Table A9). The corresponding information for the organic chemical functions can be found in Table 3. A complete topology specification (including charge-group definitions) for the 5 methylated bases can be found in the mtb files.134

of 298.15 K. Exceptions are the pure-liquid simulations of five compounds performed at a temperature Tl+g differing from 298.15 K (Table 1) for determining ρliq and ΔHvap as explained in section 2.2.1. The temperature was maintained close to its reference value by weak coupling to external heat baths157 using a coupling time of 0.1 ps. The details of the coupling (number of baths and coupled degrees of freedom) vary between the different types of simulations, and are specified in sections 2.3.2, 2.3.4, and 2.3.7. The pressure was maintained close to its reference value by isotropic weak coupling of the atomic coordinates and box dimensions to a pressure bath157 using a group-based virial and a coupling time of 0.5 ps. The isothermal compressibility involved in the pressure-coupling scheme was set to 4.575 × 10−4 kJ−1 mol nm3 for the simulations in water (as appropriate for solvated biomolecules55,158), to 11.2 × 10−4 kJ−1 mol nm3 for the simulations in cyclohexane (based on the experimental value159), and to 1.654 × 10−3 kJ−1 mol nm3 for the simulations in chloroform (based on the experimental value117). For the pureliquid simulations of the 57 organic molecules, a value of 4.575 × 10−4 kJ−1 mol nm3 was used. This choice was made for compatibility with the standard GROMOS setup, although the experimental compressibilities of the organic liquids considered typically range from 8 × 10−3 to 13 × 10−3 kJ−1 mol nm3 (see Table 2 in ref 77). Although the choice of a specific compressibility for the weak-coupling barostat algorithm affects the pressure fluctuations (magnitude and timescale), it has no influence on the averages of the calculated thermodynamic properties. Newton’s equations of motion were integrated using the leapfrog scheme160 with a timestep of 2 fs. Constraints on all bond lengths as well as the full rigidity of the water and

Figure 1. Illustration of the thermodynamic cycle employed for calibrating the charge sets of the nucleic-acid bases. For each of the five methylated bases (shown here is 9-methyladenine), a reference value ΔGref trn for the water-to-chloroform transfer free energy is calculated using TI (solute−solvent decoupling calculations in each of the two solvents) based on the reoptimized 2016H66 Lennard-Jones interaction parameters along with an initial charge set q (53A6 charges). Long plain MD simulations of the compound in the two solvents are then 2O 3 and ΔGCHCl used to estimate via OSP the free-energy changes ΔGHq→q′ q→q′ with respect to using an altered charge set q′. This results in a CHCl3 H2O predicted value ΔGtrn = ΔGref trn + ΔGq→q′ − ΔGq→q′ for the transfer free energy when employing the altered charge set q′.

They were performed in the isothermal−isobaric ensemble at a reference pressure of 1 bar and, in general, a reference temperature 3832


Article


out to evaluate the solvation free energies ΔGwat and ΔGche, respectively, of the 57 compounds (Table 1) at 1 bar and 298.15 K. The initial coordinates for these simulations were generated by randomly placing one solute molecule and 1000 or 300 molecules of water or cyclohexane, respectively, in cubic computational boxes of dimensions appropriate for the experimental densities of the pure solvents (997 kg m−3 for water and 779 kg m−3 for cyclohexane117). After energy minimization, the systems were equilibrated at constant temperature and pressure for 0.3 ns. For the free-energy calculations, the solute−solvent interactions were perturbed using a coupling parameter λ based on a soft-core scheme,168 where λ = 0 corresponds to full and λ = 1 corresponds to vanishing solute−solvent interactions. The electrostatic and Lennard-Jones soft-core parameters were set to 0.5 nm2 and 0.5, respectively. The solvation free energies were calculated using the TI method148 based on 21 equidistant λ values and trapezoidal integration. The configuration obtained after 50 ps of simulation at a given λ point was used as the initial configuration for the simulation at the next λ point. For each λ value, the system was equilibrated for 0.1 ns followed by a production simulation of 0.6 ns. The estimates of ΔGwat and ΔGche calculated in this way, expressed on a per mole basis, were compared directly to experimental data according to a standardstate definition involving identical reference concentrations for the gas-phase and solute standard states (section 2.2.1). For these calculations, the translational and internal-plus-rotational degrees of freedom of the molecules were jointly coupled to a heat bath using two separate baths, one for the solute molecule and one for the solvent molecules. 2.3.5. Organic Liquids: Error Estimates. Estimates for the statistical error on the calculated quantities were obtained by block averaging.1 For the pure-liquid properties ρliq and ΔHvap, the errors were always below 1% and are not reported. For the solvation free energies ΔGwat and ΔGche, the total error was estimated by weighted summation of error contributions evaluated separately for the simulations at the successive λ points. 2.3.6. Organic Liquids: Compatibility with Previous Work. Besides statistical errors, the simulation results reported in the present article for ρliq, ΔHvap, ΔGwat, and ΔGche in the 2016H66 force field may somewhat differ from those reported in the previous articles describing the 53A6OXY set,77 the 53A6OXY+A set,95 and the 53A6OXY+D set.96 The corresponding maximum differences over the 36 common molecules are on the order of 1 kg m−3, 1 kJ mol−1, 2 kJ mol−1, and 2−3 kJ mol−1, respectively. This is in part due to the small difference in the parameters for atom type CH2r, which affects the solvent cyclohexane (and thus slightly alters all ΔGche values) as well as the compound CHXL. However, it is also due to slight changes in the simulation protocol, namely: (1) the pressure coupling in the pure-liquid simulations relies on a common compressibility (section 2.3.1; vs experimental value for the specific liquid in the previous articles); (2) the temperature coupling in the pure-liquid simulations relies on a joint coupling of the translational and internal-plus-rotational degrees of freedom of the molecules (section 2.3.2; vs separate coupling to two heat baths in the previous articles); (3) the gas-phase simulations rely on single-molecule SD simulations (section 2.3.3; vs MD simulations of 512 molecules placed at noninteracting distances with separate coupling of the translational and internal-plus-rotational degrees of freedom to two heat baths in the previous articles). The first two changes have a negligible impact on the results. However, the third change leads to noticeable differences in the calculated ΔHvap (see below).

chloroform molecules were enforced by application of the SHAKE procedure161 with a relative geometric tolerance of 10−4. The center of mass motion was removed every 100 ps. The nonbonded interactions were calculated using a twinrange scheme,162 with short- and long-range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of 5 timesteps for the short-range pairlist and intermediate-range interactions. A reaction-field correction75,76 was applied to account for the mean effect of electrostatic interactions beyond the long-range cutoff distance. The corresponding relative dielectric permittivities were set to 61 for the aqueous systems (as appropriate for the SPC model163), to 6 for the simulations in cyclohexane (the same value as used in ref 63, which is somewhat higher than the experimental value164 of 2.05 to avoid cutoff artifacts), to 4.7 for the simulations in chloroform (based on the experimental value117), or to the experimental permittivity ε (Table 1) for the pure-liquid simulations of the 57 organic molecules. The reaction-field self-term and excludedatom-term contributions165 to the energy, forces, and virial were included as described in ref 166. Configurations were saved to file every 5 ps for analysis, except for the simulations involved in the viscosity calculations (every timestep, see section 2.4.6). 2.3.2. Organic Liquids: Pure-liquid Simulations. The pureliquid simulations of the 57 compounds considered (Table 1) were performed at 298.15 K or, in five cases, at the indicated temperature Tl+g. The initial coordinates for these simulations were generated by randomly placing 512 molecules in cubic computational boxes of dimensions appropriate for the experimental liquid density ρliq. After energy minimization, the systems were equilibrated at constant volume (0.25 ns) and then at constant pressure (1 ns). The production simulations at constant pressure used for the calculation of ρliq and ΔHvap lasted 5 ns after equilibration. For these simulations, the translational and internal-plus-rotational degrees of freedom of all molecules were jointly coupled to a heat bath, using a single bath for the entire system. 2.3.3. Organic Liquids: Gas-Phase Simulations. The estimation of the gas-phase intramolecular energies is required for the calculation of ΔHvap. For compatibility with the corresponding pure-liquid simulations, the gas-phase simulations were performed at 298.15 K or, in five cases, at the temperature Tl+g (Table 1). These simulations were carried out using SD for a single molecule in vacuum with a friction coefficient γ = 91 ps−1 mimicking water154 and lasted 500 ns. Note that the choice of a specific value for γ affects the dynamics of the system, but not the averages (and fluctuations) of the calculated thermodynamic properties. All the other simulation parameters including the nonbonded interaction scheme (with εRF = 1) were the same as for the liquid simulations (see section 2.3.1). The vaporization enthalpy ΔHvap was calculated as the average potential energy per molecule in the gas phase (SD simulation) minus the corresponding value in the pure liquid (MD simulation; see section 2.3.2), expressed on a per mole basis and increased by RT, where T is 298.15 K or Tl+g. This corresponds to a standard-state definition involving a reference pressure of 1 bar for both the (ideal) gas phase and the (real) liquid phase as recommended167 by IUPAC. The term Mmolρ−1 liq accounting for the pressure−volume enthalpy contribution of the pure liquid, where Mmol is the molar mass of the compound, is in all cases very small and was neglected. 2.3.4. Organic Liquids: Solvation Free-Energy Calculations. Simulations in water and in cyclohexane were carried 3833


Article


This permitted to predict ΔGtrn values for the bases described with the same Lennard-Jones interaction parameters along with arbitrary charge sets. In this way, the partial charges could be optimized to reproduce the experimental transfer free energies ΔGtrn. Additional TI calculations were performed using the 53A6 parameters (Lennard-Jones and charges) for comparison purposes, and using the final 2016H66 parameters so as to validate the estimates obtained using OSP. 2.4. Analysis of Other Liquid Properties. Additional simulations were performed to calculate seven other thermodynamic, transport, and dielectric properties of the pure liquids (besides ρliq and ΔHvap). These properties are the molar isobaric heat capacity cP, the isothermal compressibility κT, the isobaric thermal-expansion coefficient αP, the surface tension coefficient γ, the self-diffusion coefficient D, the shear viscosity η, and the static relative dielectric permittivity ε. Unless otherwise specified, the setup of these simulations (box shape, number of molecules, pressure and temperature coupling, and other simulation parameters) was the same as those described in sections 2.3.1 and 2.3.2, and the production simulations for these calculations lasted 5 ns after equilibration (2 ns for η, see section 2.4.6; 60 ns for ε, see section 2.4.7). Constant-pressure simulations were performed at two different temperatures so as to calculate cP and αP by finite difference, namely at T1 = 288.15 K and T2 = 308.15 K or, in five cases (Table 1), at T1 = Tl+g − 10 K and T2 = Tl+g + 10 K. Constant-volume simulations were performed at the equilibrium density ρliq calculated based on the constant-pressure simulations. In addition, constant-volume simulations were performed at two different densities so as to calculate κT by finite difference, namely at ρ1 = ρliq − 25 kg m−3 and ρ2 = ρliq + 25 kg m−3. All the calculated quantities are reported along with an estimate of the statistical error, calculated as the standard deviation over five independent simulation repeats (for D and ε, see sections 2.4.5 and 2.4.7), or using block averaging1 over a single simulation (for all other quantities). For each property, experimental results were usually only available for a subset of the 57 molecules. The corresponding experimental values are listed in Table B3 based on refs 116−133, 145, 175−187. 2.4.1. Molar Isobaric Heat Capacity. The classical188−191 (i.e., without quantum corrections) molar isobaric heat capacity cP at 298.15 K and 1 bar was determined based on two constant-pressure simulations carried out at temperatures T1 and T2 via finite difference as192

Similar considerations apply to the results reported here using the 53A5 and 53A6 force fields for comparison purposes. These were recalculated as well using the same SD protocol as that employed for 2016H66, and the results may also somewhat differ from those reported in the original article describing these two parameter sets63 (where the gas-phase calculations relied on MD simulations of 512 molecules placed at noninteracting distances with a joint coupling of the translational and internalplus-rotational degrees of freedom to a single heat bath). Independent investigations (data not shown) revealed that the MD protocol results in a progressive flow of kinetic energy from the intramolecular degrees of freedom to the translational (only when using a single thermostat) and rotational (when using either one or two thermostats) ones. Similar equipartition violations induced by velocity-scaling thermostats have been reported previously.169−174 This effect results in lower intramolecular temperatures (vibrations and internal rotations) and in an underestimation of the potential energy. It should be stressed, however, that the errors incurred in the parametrizations of refs 63, 77, 95, and 96 against ΔHvap remain limited (typical error of about 3 kJ mol−1 or less), because these studies relied on gas-phase MD simulations of very short effective durations (50 ps between successive velocity reassignments from a Maxwell−Boltzmann distribution), i.e., the decrease of the intramolecular temperature remained limited. In any case, the SD protocol employed systematically in this article for the three force fields guarantees an equipartition of the kinetic energy across all molecular degrees of freedom, and the reported results reflect the true accuracy of these parameter sets in terms of ΔHvap. 2.3.7. Nucleic-Acid Bases. The calculation of the water-tochloroform transfer free energies ΔGtrn for the methylated nucleic-acid bases relied on the thermodynamic cycle illustrated in Figure 1. The initial coordinates for these simulations were generated by randomly placing one solute molecule and about 1300 water or 300 chloroform molecules (exact numbers depending on the system), respectively, in cubic computational boxes of dimensions (edge length of about 3.5 nm) appropriate for the experimental densities of the pure solvents (997 kg m−3 for water and 1492 kg m−3 for chloroform117). After energy minimization, the systems were equilibrated at constant temperature and pressure for 1 ns. For the calculations relying on the TI method,148 the solute−solvent interactions were perturbed using a coupling parameter λ based on a soft-core scheme,168 where λ = 0 corresponds to full and λ = 1 to vanishing solute−solvent interactions. The electrostatic and Lennard-Jones soft-core parameters were set to 0.5 nm2 and 0.5, respectively. The TI simulations were carried out at 21 equidistant λ values, and the integration relied on the trapezoidal rule. The configuration obtained after 0.1 ns of simulation at a given λ point was used as the initial configuration for the simulation at the next λ point. For each λ point, the system was equilibrated for 0.1 ns followed by a production simulation of 1 ns. For these calculations, the translational and internal-plus-rotational degrees of freedom of the molecules were jointly coupled to a heat bath using two separate baths, one for the solute molecule and one for the solvent molecules. Such a TI calculation was performed to determine ΔGref trn considering the 2016H66 Lennard-Jones interaction parameters along with the 53A6 charge set. The simulations of the bases in water and in chloroform at λ = 0 were then prolonged to 10 ns and served as reference simulations for the calculaH2O CHCl and ΔGq→q′ 3 using the OSP method.149,150 tion of ΔGq→q′

cP =

⟨U ⟩T2 − ⟨U ⟩T1 ⎛ ∂H ⎞ ⎛ ∂U ⎞ ⎜ ⎟ ≈ ⎜ ⎟ ≈ ⎝ ∂T ⎠ P ⎝ ∂T ⎠ P T2 − T1

where H = U + PV is the enthalpy, U is the total (kinetic plus potential) energy of the system expressed on a per molecule basis, and ⟨ ··· ⟩T stands for ensemble averaging at temperature T. The contribution of the pressure−volume term PV is on the order of 10−3 J K−1 mol−1 for the liquids considered and was neglected. 2.4.2. Isothermal Compressibility. The isothermal compressibility κT at 298.15 K and 1 bar was determined based on two constant-volume simulations carried out at densities ρ1 and ρ2 via finite difference as136 κT = −

ln ρ2 − ln ρ1 1 ⎛⎜ ∂V ⎞⎟ ≈ ⎝ ⎠ V ∂P T ⟨P⟩ρ2 − ⟨P⟩ρ1

where P is the instantaneous pressure and ⟨ ··· ⟩ρ stands for ensemble averaging at density ρ. 3834


Article

Journal of Chemical Theory and Computation 2.4.3. Isobaric Thermal-Expansion Coefficient. The isobaric thermal-expansion coefficient αP at 298.15 K and 1 bar was determined based on two constant-pressure simulations carried out at temperatures T1 and T2 via finite difference as136

ε=

where ρ is the instantaneous density and ⟨ ··· ⟩T stands for ensemble averaging at temperature T. 2.4.4. Surface Tension Coefficient. The surface tension coefficient γ at 298.15 K and 1 bar was determined based on a constant-volume simulation of a rectangular computational box (edge lengths Lx = Ly = L and Lz = 5L) presenting a liquid− vapor interface in the xy plane as193

3. RESULTS AND DISCUSSION 3.1. Organic Liquids: Reproduction of the Target Properties. The results obtained using the final 2016H66 parameter set in terms of the four target properties ρliq, ΔHvap, ΔGwat, and ΔGche for the 57 organic liquids considered are compared to experimental data in Figure 2. Experimental values were available for all compounds, except for ΔGwat (2 values missing) and ΔGche (24 values missing). The corresponding numerical values can be found in SI B (experimental values in Table B3, calculated values and deviations from experiment in Table B4). The results obtained using the previous versions63 53A5 and 53A6 of the GROMOS force field are also compared to experimental data in Figure 2 (gray symbols), and the corresponding numerical values are reported in SI B (Tables B5 and B6). Considering the three parameter sets, a simplified comparison in terms of the root-mean-square deviations (RMSD) and average deviations (AVED) from experiment, either sorted by chemical function or considering the entire set of compounds (calibration as well as validation compounds in Table 1), is also provided in Table 5. In the original article63 reporting the 53A5 and 53A6 parameter sets, which considered a smaller collection of compounds and relied on slightly different simulation protocols (see section 2.3.6), set 53A5 was found to reproduce ρliq and ΔHvap very well, but to systematically overestimate ΔGwat. Set 53A6 was found to improve the agreement with experiment in terms of ΔGwat, but at the cost of systematically overestimating ΔHvap. These observations still essentially hold based on the values recalculated for these two sets in the present work. However, compared to both 53A5 and 53A6, the 2016H66 set improves the agreement with experiment in terms of RMSD and AVED over the entire set of compounds and in terms of all four target properties. Considering ρliq, the 2016H66 set leads to an overall RMSD of 32.4 kg m−3 (vs 38.9 and 49.8 for 53A5 and 53A6, respectively) and an AVED of 1.0 kg m−3 (vs 16.2 and 26.6 for 53A5 and 53A6, respectively). The largely reduced AVED suggests that a systematic error (overestimation of ρliq) has been corrected. However, the residual errors in 2016H66 are not distributed entirely homogeneously across the chemical functions. In particular, the density is always slightly underestimated for the alcohols (except MTL), ethers, and esters, whereas it is always slightly overestimated for the amines, thiols, and sulfides. These deviations nevertheless remain below 10% in magnitude, except for the diamine EDAN (+13%) and the thiol ETSH (+10%). Considering ΔHvap, the 2016H66 set leads to an overall RMSD of 3.5 kJ mol−1 (vs 4.3 and 7.5 for 53A5 and 53A6, respectively) and an AVED of 0.2 kJ mol−1 (vs −0.4 and 4.2 for 53A5 and 53A6, respectively). The value of the AVED is again close to zero for 2016H66, suggesting that the residual error is

Lz 1 Pzz − (Pxx + Pyy) 2 2

where Pαα (α = x, y, z) is a diagonal element of the instantaneous pressure tensor. In practice, the system was first equilibrated in a cubic box of volume V = L3. The box was then extended in the z direction, thereby creating a liquid−vapor interface in the xy plane. 2.4.5. Self-Diffusion Coefficient. The self-diffusion coefficient D at 298.15 K and 1 bar was determined based on a constant-pressure simulation from the long-time limit of the mean-square displacement of the molecules using the Einstein relation4,194 D = lim

⟨(ri(τ + t ) − ri(τ ))2 ⟩i , τ 6t

t →∞

where ri is the instantaneous position of the center of geometry of molecule i (following molecules across periodic boundaries) and ⟨ ··· ⟩i,τ stands for averaging over all molecules i and time origins τ. In practice, a least-squares fitting over the interval 0−2.5 ns was performed to estimate the value of D. The resulting correlation coefficients R2 were at least 0.99 in all cases. The result was averaged over 5 independent repeats (different initial configurations and velocities) and is reported along with the corresponding standard deviation as an error estimate. 2.4.6. Shear Viscosity. The shear viscosity η at 298.15 K and 1 bar was determined based on a constant-volume simulation as described in ref 195 from the displacement of the offdiagonal components Pαβ (α, β = x, y, z) of the pressure tensor ΔPαβ(t ) =

∫0

t

dt ′Pαβ(t ′)

via the Einstein-type relation η=

3εo(2εRF + 1)RT ⟨V ⟩ − [⟨M2⟩ − ⟨M⟩2 ]

where εo is the permittivity of vacuum, V is the box volume, and εRF is the reaction-field permittivity (set equal to the experimental value for the liquid; Table 1). These simulations were carried out for 60 ns. The result was averaged over 5 independent repeats (different initial configurations and velocities) and is reported along with the corresponding standard deviation as an error estimate.

ln⟨ρ⟩T2 − ln⟨ρ⟩T1 1 ⎛ ∂V ⎞ αP = ⎜ ⎟ ≈ − ⎝ ⎠ V ∂T P T2 − T1

γ=

3εo(2εRF + 1)RT ⟨V ⟩ + 2εRF[⟨M2⟩ − ⟨M⟩2 ]

V d 2 lim ⟨ΔPαβ (t )⟩ 2RT t →∞ dt

where V is the box volume. In practice, the infinitesimal dt was approximated by the timestep Δt, and a least-squares fitting over the interval 5−20 ps was performed to estimate the value of η. The resulting correlation coefficients R2 were at least 0.99 in all cases. These simulations were carried out for a duration of 2 ns, over which trajectory configurations were saved to file every single timestep. 2.4.7. Static Relative Dielectric Permittivity. The static relative dielectric permittivity ε at 298.15 K and 1 bar was determined based on a constant-pressure simulation from the fluctuations of box dipole-moment vector M using the Neumann relation196 3835


Article


Figure 2. Comparison between experimental and simulated properties calculated using the parameter sets 53A5, 53A6, and 2016H66 for 57 organic molecules. The properties considered are the liquid density (ρliq), vaporization enthalpy (ΔHvap), hydration free energy (ΔGwat), and solvation free energy in cyclohexane (ΔGche). The compounds considered are listed in Table 1. The crosses indicate compounds used for calibration, and the circles, squares, and diamonds indicate compounds considered for validation. The straight line marks the ideal case of a perfect agreement between simulation and experiment. Note that 2 and 24 compounds are omitted in the graphs for ΔGwat and ΔGche, respectively, owing to the lack of experimental data. The values are reported numerically in SI B (experimental values in Table B3, simulated values and deviations from experiment in Tables B4−B6 for 2016H66, 53A5, and 53A6, respectively).

represents an improvement in terms of AVED, but still shows a large RMSD value. The RMSD is reduced in 2016H66 compared to 53A6, although the analysis per chemical function reveals again that the residual errors are not equally distributed. In particular, the hydration free energy is still systematically overestimated for the amides and, to a lesser extent, for the amines. The most significant deviations in terms of ΔGwat (5 kJ mol−1, i.e. about 2RT, or larger) are observed for the alcohols CHXL (+6.4), 2M2P (+6.9), and 2M2B (+6.5), the amines DEAN (+5.6) and TEAN (+7.3), the amides AMD (+5.8), LNMA (+13.4), and DAMD (+14.2), and the disulfide DMDS (+5.0). Considering ΔGche, the three parameter sets provide results that are in very good agreement with experimental data. The 2016H66 set leads to an overall RMSD of 2.1 kJ mol−1 (vs 3.1 and 3.2 for 53A5 and 53A6, respectively) and an AVED of 1.0 kJ mol−1 (vs −1.1 and −1.2 for 53A5 and 53A6, respectively). The 2016H66 set mainly improves over 53A5 and 53A6 for the ketones and esters. For the remaining chemical functions, the three parameter sets perform similarly well.

essentially statistical and no longer systematic. However, here also, the errors are not distributed entirely homogeneously across the chemical functions. The vaporization enthalpy is always slightly underestimated for the secondary and tertiary alcohols, and slightly overestimated for MTL and the carboxylic acids. The most significant deviations in terms of ΔHvap (5 kJ mol−1, i.e. about 2RT, or larger) are observed for the alcohols 2M2P and 2M2B (−8.3 and −6.4, respectively), the carboxylic acid BTA (+6.2), the diamine EDAN (+10.5), and the aromatic compounds PHL and PHT (+12.5 and +8.8, respectively). Note that the ΔHvap errors for these compounds are also comparatively large for 53A5 and 53A6 (except for BTA in 53A5, as well as 2M2P and 2M2B in 53A6, where the errors are smaller compared to those of 2016H66). Considering ΔGwat, the 2016H66 set leads to an overall RMSD of 4.1 kJ mol−1 (vs 8.0 and 6.3 for 53A5 and 53A6, respectively) and an AVED of 2.6 kJ mol−1 (vs 3.8 and −0.3 for 53A5 and 53A6, respectively). The positive AVED for 53A5 suggests that the organic molecules considered are not sufficiently polar within this set. Compared to 53A5, 53A6 3836


Article


Table 5. Comparison between experimental and simulated target properties calculated using the parameter sets 53A5, 53A6, and 2016H66 for 57 organic moleculesa chemical function alcohols

N 15

ethers

4

aldehydes ketones

3 6

acids

3

esters

4

amines

6

amides

3

thiols

3

sulfides/disulfides

4

aromatic compounds

6

entire set of compounds

57

parameter set 53A5 53A6 2016H66 53A5/6 2016H66 2016H66 53A5/6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66

ρliq (kg m−3) 28.8 31.4 34.0 13.8 23.3 10.0 18.4 13.7 78.9 69.1 19.6 16.3 48.3 25.8 42.5 73.3 50.2 55.1 100.2 12.0 49.3 61.8 54.4 35.2 55.6 43.0 54.0 47.6 23.8 38.9 49.8 32.4

[−19.6] [−11.0] [−24.1] [−0.4] [−21.0] [−1.9] [11.2] [1.2] [70.6] [61.1] [5.2] [15.7] [47.6] [−24.9] [33.1] [64.0] [32.7] [55.1] [93.1] [−0.6] [46.3] [58.1] [48.0] [24.7] [50.3] [37.1] [49.7] [32.8] [15.8] [16.2] [26.6] [1.0]

ΔHvap (kJ mol−1) (15) (15) (15) (4) (4) (3) (6) (6) (3) (3) (3) (4) (4) (4) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (52) (52) (57)

3.3 4.2 3.3 1.2 0.7 2.6 1.8 1.3 5.6 11.6 4.4 0.7 9.9 1.9 3.8 11.1 4.4 2.7 15.6 2.9 2.5 1.0 1.4 3.8 2.4 0.8 10.1 10.4 6.6 4.3 7.5 3.5

[−1.3] [2.8] [−1.4] [−1.2] [0.3] [2.6] [−1.7] [−1.0] [−5.3] [11.5] [4.1] [0.2] [9.9] [−1.8] [2.0] [10.0] [2.5] [−0.1] [12.2] [−1.2] [−2.5] [−0.9] [−1.3] [−3.8] [0.4] [0.0] [8.2] [6.9] [2.3] [−0.4] [4.2] [0.2]

(15) (15) (15) (4) (4) (3) (6) (6) (3) (3) (3) (4) (4) (4) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (52) (52) (57)

ΔGwat (kJ mol−1) 3.4 2.9 3.8 10.4 1.3 1.6 4.8 2.1 7.7 5.5 3.0 8.7 8.6 2.1 12.7 5.9 4.8 16.2 13.3 11.8 0.8 1.9 1.8 4.6 6.3 2.8 10.9 6.4 1.7 8.0 6.3 4.1

[2.3] [−1.4] [3.0] [8.9] [−0.1] [1.3] [4.7] [1.9] [0.4] [−5.5] [2.9] [8.7] [−8.5] [2.0] [11.5] [4.8] [4.5] [16.1] [7.2] [11.1] [−0.8] [−1.8] [1.5] [4.5] [−2.6] [1.1] [−9.2] [−4.9] [−0.2] [3.8] [−0.3] [2.6]

(15) (15) (15) (4) (4) (3) (5) (5) (3) (3) (3) (4) (4) (4) (5) (5) (5) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (50) (50) (55)

ΔGche (kJ mol−1) 1.8 1.9 2.1 3.4 0.6 2.4 1.1 4.9 5.4 2.7 6.0 6.1 1.0 2.0 1.8 2.3 2.7 2.8 4.1 2.1 2.5 2.7 0.5 1.1 1.9 2.0 2.3 2.6 3.1 3.2 2.1

[0.2] [0.0] [1.3] [−3.4] [−0.6] (0) [−2.3] [0.9] [−4.0] [−4.5] [−0.3] [−6.0] [−6.0] [−0.0] [1.4] [1.5] [2.0] [2.7] [2.8] [4.1] [−2.1] [−2.5] [−2.7] [−0.5] [−1.1] [−1.9] [1.3] [1.3] [1.9] [−1.1] [−1.2] [1.0]

(9) (9) (9) (1) (1) (5) (5) (2) (2) (2) (4) (4) (4) (3) (3) (3) (1) (1) (1) (1) (1) (1) (1) (1) (1) (5) (5) (6) (32) (32) (33)

a The properties considered are the pure-liquid density (ρliq), vaporization enthalpy (ΔHvap), hydration free energy (ΔGwat), and solvation free energy in cyclohexane (ΔGche). The compounds considered are listed in Table 1. For each property and chemical function, the RMSD between simulation and experiment as well as the corresponding AVED (simulation minus experiment, between square brackets) are reported. The number of compounds for each chemical function is also indicated (N) as well as the number of compounds for which experimental values for a given property were available (between parentheses). Note that 53A5 and 53A6 are equivalent for the ethers and the ketones, and do not encompass parameters for the aldehydes as well as the compounds DAMD and PYRI. The charge set for the esters in 53A6 is the one of ref 62 (entry q1 of Table 2 therein), and the charge set for the amines in 53A6 is the one of ref 156 (Table 1 therein). The detailed results for the 57 compounds can be found in SI B (experimental values in Table B3, simulated values and deviations from experiment in Tables B4−B6 for 2016H66, 53A5, and 53A6, respectively).

The deviations in terms of ΔGche in 2016H66 are below 5 kJ mol−1 (i.e., about 2RT) for all compounds considered. 3.2. Organic Liquids: Analysis of Additional Nontarget Properties. The 2016H66 force field was tested in terms of seven additional pure-liquid properties besides ρliq and ΔHvap, namely, the (classical) isobaric heat capacity cP, the isothermal compressibility κT, the isobaric thermal-expansion coefficient αP, the surface tension coefficient γ, the self-diffusion coefficient D, the shear viscosity η, and the static relative dielectric permittivity ε. The results for the 57 organic compounds considered are compared to experimental data in Figure 3. The corresponding numerical values can be found in SI B (experimental values in Table B3, calculated values and deviations from experiment in Table B4). Experimental values for γ and ε were available for all compounds, whereas 1, 21, 12, 29, and 1 values were missing for cP, κT, αP, D, and η, respectively. The results obtained using the previous versions63 53A5 and 53A6 of the GROMOS force field are also compared

to experimental data in Figure 3 (gray symbols), and the corresponding numerical values are reported in SI B (Tables B5 and B6). Considering the three parameter sets, an overall assessment in terms of RMSD and AVED, either sorted by chemical function or considering the entire set of compounds (calibration as well as validation compounds in Table 1), is provided in Table 6. Because these seven properties were not considered as target for the calibration of 2016H66, their comparison with experimental data belongs to the validation of the force field. In terms of the thermodynamic properties cP, κT, αP, and γ, the results obtained with 2016H66 can be considered satisfactory. The largest deviations are observed for cP, and are particularly large for the amines, the amides, and the aromatic compounds. For these compounds, cP is systematically overestimated, sometimes by a very large amount (e.g., nearly a factor of 2 for benzene), and in a rather nonregular fashion across the different molecules. However, if these three classes 3837


Article


Figure 3. Comparison between experimental and simulated properties calculated using the parameter sets 53A5, 53A6, and 2016H66 for 57 organic molecules. The properties considered are the (classical) isobaric heat capacity (cP), isothermal compressibility (κT), isobaric thermal-expansion coefficient (αP), surface tension coefficient (γ), self-diffusion coefficient (D), shear viscosity (η), and static relative dielectric permittivity (ε). The compounds considered are listed in Table 1. The crosses indicate compounds used for calibration, and the circles, squares, and diamonds indicate compounds considered for validation. The straight line marks the ideal case of a perfect agreement between simulation and experiment. Note that some compounds are omitted from specific graphs due to the lack of experimental data. In addition, CHXL is omitted from the η graph and AMD, LNMA and DAMD are omitted from the ε graph owing to very high experimental values. LNMA is also omitted from the αP graph owing to a very high simulated value. The values are reported numerically in SI B (experimental values in Table B3, simulated values and deviations from experiment in Tables B4−B6 for 2016H66, 53A5, and 53A6, respectively). 3838


3839

53A5 53A6 2016H66 53A5/6 2016H66 2016H66 53A5/6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66 53A5 53A6 2016H66

parameter set

21.3 32.5 22.3 46.8 33.5 11.9 37.7 34.8 12.0 8.1 10.8 29.5 22.9 27.4 45.9 52.6 62.8 70.3 66.3 54.2 13.0 8.1 10.5 44.4 40.7 41.1 103.0 99.5 103.8 46.5 47.6 47.4

[−14.6] [−25.8] [−14.8] [−45.3] [−32.3] [−7.3] [−36.4] [−33.1] [11.0] [−0.3] [1.6] [−28.7] [−22.0] [−26.6] [17.9] [22.5] [38.3] [47.3] [34.3] [19.6] [−12.3] [−6.6] [−7.4] [−44.1] [−40.4] [−40.9] [93.9] [94.9] [98.7] [−4.0] [−6.7] [1.0]

cP (J mol−1 K−1) (15) (15) (15) (4) (4) (3) (5) (5) (3) (3) (3) (4) (4) (4) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (51) (51) (56) 17.0 8.4 2.3 8.4 15.4 3.6

3.0 1.9 12.9 22.4 1.3 0.3 2.0 2.8 13.1 52.6 10.0 17.1 17.4 6.6

2.3 5.4 2.4 1.5 3.5

[−5.9] [−2.6] [1.8] [−1.3] [0.6] [−0.6]

[1.5] [0.4] [−4.9] [−11.0] [0.5] [−0.3] [1.8] [−2.7] [−7.7] [14.6] [−5.0] [17.1] [17.4] [4.8]

[−0.8] [0.5] [−1.2] [−1.4] [−3.4]

κT (10−5 bar−1) (14) (14) (14) (2) (2) (0) (5) (5) (3) (3) (3) (3) (3) (3) (2) (2) (2) (1) (1) (2) (0) (0) (0) (0) (0) (0) (4) (4) (5) (34) (34) (36) 1.5 0.8 0.5 2.1 2.0 2.2 3.8 4.0 4.4

1.4 1.8 1.8 1.3 0.7 0.9 0.5 0.5 1.2 1.4 2.0 2.1 3.5 0.9 2.1 3.9 1.9 21.9 21.1 19.6

[1.5] [−0.8] [−0.5] [0.6] [1.2] [1.7] [0.7] [0.3] [1.2]

[0.7] [0.5] [1.0] [−1.2] [−0.3] [−0.8] [−0.4] [0.1] [−1.1] [−0.7] [0.9] [−2.1] [−3.4] [−0.9] [1.2] [−0.6] [1.1] [21.9] [21.1] [13.8]

αP (10−4 K−1) (15) (15) (15) (2) (2) (3) (4) (4) (3) (3) (3) (4) (4) (4) (5) (5) (5) (1) (1) (2) (0) (0) (0) (1) (1) (1) (5) (5) (6) (40) (40) (45)

2.8 2.8 3.3 1.2 0.6 3.3 1.2 1.2 6.0 5.2 1.5 3.4 7.9 0.9 1.8 8.6 4.4 3.5 0.7 6.3 1.8 0.6 1.0 2.1 3.5 2.4 5.9 3.9 4.8 3.2 4.4 3.3

[−2.2] [−2.0] [−2.8] [1.2] [−0.1] [−2.5] [0.7] [−0.4] [5.0] [4.6] [−0.8] [3.4] [7.7] [−0.9] [0.1] [5.1] [−0.6] [−3.0] [−0.2] [−4.5] [−1.6] [−0.6] [0.2] [−2.1] [1.6] [1.9] [2.0] [−3.2] [−3.6] [−0.0] [0.8] [−1.6]

γ (mN m−1) (15) (15) (15) (4) (4) (3) (6) (6) (3) (3) (3) (4) (4) (4) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (52) (52) (57) 0.6 0.6 0.5 0.6 0.6 0.5

0.5 0.6 0.8 0.8 0.5 1.0 1.0 0.3 0.7 0.7 0.7 0.3 0.3 0.3

0.6 0.6 0.5 0.5 0.3

[−0.5] [−0.5] [−0.4] [−0.4] [−0.4] [−0.2]

[−0.3] [−0.5] [−0.8] [−0.8] [−0.5] [−0.9] [−0.9] [0.3] [−0.7] [−0.7] [−0.6] [−0.3] [−0.3] [0.3]

[−0.3] [−0.3] [−0.0] [0.4] [−0.3]

(10) (10) (10) (2) (2) (0) (3) (3) (1) (1) (1) (2) (2) (2) (2) (2) (2) (1) (1) (2) (0) (0) (0) (0) (0) (0) (5) (5) (6) (26) (26) (28)

D (10−9 m2 s−1) 27.3 26.3 28.1 0.6 0.2 0.9 0.8 0.6 4.0 2.0 3.3 1.5 1.5 1.6 16.0 57.0 28.3 13.6 47.1 12.7 1.3 1.5 0.6 1.5 1.5 0.9 10.6 32.6 20.6 16.1 27.9 18.6

[−23.0] [−21.6] [−23.7] [−0.5] [−0.1] [−0.4] [−0.7] [−0.3] [3.4] [−0.8] [−2.6] [−1.3] [0.5] [−1.4] [6.2] [27.1] [11.2] [−13.6] [24.6] [−9.6] [−1.3] [−1.4] [−0.5] [−1.0] [−1.2] [0.2] [1.9] [1.5] [−5.7] [−6.2] [−2.0] [−6.3]

(14) (14) (14) (4) (4) (3) (6) (6) (3) (3) (3) (3) (3) (3) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (50) (50) (55)

η [10−4 kg.m−1.s−1] 8.5 7.9 6.8 2.9 9.7 1.9 8.7 5.4 3.0 1.8 2.3 3.9 6.9 1.1 4.2 7.5 1.9 73.7 74.3 56.9 1.2 0.6 0.7 4.6 8.4 4.3 2.4 3.0 3.2 15.7 16.0 14.0

ε [−8.1] [−7.3] [−6.4] [−2.5] [8.7] [0.2] [−8.4] [−5.0] [−2.7] [−0.4] [−2.2] [−3.9] [6.2] [−1.0] [3.0] [6.5] [−1.2] [−67.6] [−68.1] [−47.6] [−1.2] [−0.4] [−0.5] [−3.7] [3.9] [1.6] [−1.6] [−0.7] [−2.3] [−6.7] [−4.5] [−4.6]

(15) (15) (15) (4) (4) (3) (6) (6) (3) (3) (3) (4) (4) (4) (6) (6) (6) (2) (2) (3) (3) (3) (3) (4) (4) (4) (5) (5) (6) (52) (52) (57)

The pure-liquid properties considered are the (classical) isobaric heat capacity (cP), isothermal compressibility (κT), isobaric thermal-expansion coefficient (αP), surface tension coefficient (γ), selfdiffusion coefficient (D), shear viscosity (η), and static relative dielectric permittivity (ε). The compounds considered are listed in Table 1. For each property and chemical function, the RMSD between simulation and experiment as well as the corresponding AVED (simulation minus experiment, between square brackets) are reported. Note that CHXL was omitted from the statistics for η owing to its anomalously high experimental value. The number of compounds for each chemical function is also indicated (N) as well as the number of compounds for which experimental values for a given property were available (between parentheses). Note that 53A5 and 53A6 are equivalent for the ethers and the ketones, and do not encompass parameters for the aldehydes as well as the compounds DAMD and PYRI. The charge set for the esters in 53A6 is the one of ref 62 (entry q1 of Table 2 therein), and the charge set for the amines in 53A6 is the one of ref 156 (Table 1 therein). The detailed results for the 57 compounds can be found in SI B (experimental values in Table B3, simulated values and deviations from experiment in Tables B4−B6 for 2016H66, 53A5, and 53A6, respectively).

a

57

entire set of compounds

3

amides

6

6

amines

aromatic compounds

4

esters

4

3

acids

sulfides/ disulfides

3 6

aldehydes ketones

3

4

ethers

thiols

15

N

alcohols

chemical function

Table 6. Comparison between Experimental and Simulated Validation Properties Calculated using the parameter sets 2016H66, 53A5, and 53A6 for 57 Organic Moleculesa

Journal of Chemical Theory and Computation Article


Article

Journal of Chemical Theory and Computation are disregarded, the extent of correlation between simulation and experiment becomes very good, albeit with a systematic offset of about −20 J mol−1 K−1 in the calculated values. These discrepancies can at least in part be attributed to the absence of quantum corrections to the results calculated at the classical level.188−191 Concerning this issue, it is interesting to compare our observations concerning 2016H66 to those of ref 191 (see Figure 7 therein), where the classical molar heatcapacities calculated using the OPLS/AA force field37,41,43,44 and the GAFF force field22 for 130 liquids were found to systematically overestimate the experimental values by about a factor of 2. The difference probably comes from the fact that the two force fields considered are all-atom force fields whereas 2016H66 is a united-atom force field. This implicitly eliminates the contribution of high-frequency vibrations involving the aliphatic hydrogen atoms, leading to more realistic classical cP estimates. This interpretation is in line with the overestimation observed in 2016H66 for the aromatic compounds, for which the hydrogen atoms are explicit. Overestimation for the amines and amides may be due as well to the presence of explicit hydrogen atoms, along with peculiar vibrational properties of these groups (amide vibration bands, amine pyramidal inversion). Although classical heat capacities may be of moderate relevance in terms of a detailed comparison with experimental data in the absence of quantum corrections, they nevertheless determine the temperature derivatives of the thermodynamic properties within the force field. In this sense, the qualitative agreement with experimental data observed for most compounds suggests a good transferability of the 2016H66 force field to temperatures differing from 298.15 K (within reasonable bounds), which is actually an advantage of the united-atom over an all-atom representation. The calculated values of κT present in general an excellent agreement with the experimental ones except for two noticeable outliers, the amide LNMA and the aromatic amine ANLN, where κT is overestimated by about a factor of 2. The excellent agreement between simulation and experiment in terms of κT, in line with the results of ref 191 (see Figure 6 therein) for the OPLS/AA force field37,41 and the GAFF force field,22 suggests that the accurate representation of both ρliq and κT is possible given the employed force-field functional form and, in particular, the inverse-twelfth power assigned to the repulsion component of the Lennard-Jones function.71 This is an interesting observation, considering that the latter choice is ad hoc and has sometimes been suggested to be too steep, early force fields commonly relying on an inverse-ninth power instead.72,73,197−200 Considering αP, the overall agreement between calculated and experimental values is also good, except for the compounds MTL, ACA, EDAN, BZN, and ANLN, where αP is overestimated by 30−50%, and the amide LNMA, where αP is overestimated by about a factor of 4, owing to which the representative point for this compound is not displayed in Figure 3. Finally, considering γ, the overall agreement between calculated and experimental values is excellent. The largest deviations of about 25−30% are observed for the compounds 2M2P, EAN, and AMD. It has been recently suggested201−203 that the application of a long-range Lennard-Jones correction or the use of a lattice-sum Lennard-Jones representation is required to accurately reproduce γ. However, the present agreement with experiment in terms of γ suggests that it is not necessarily the case. With the Lennard-Jones cutoff of 1.4 nm employed here, a careful parametrization against ΔHvap leads to effective C61/2

dispersion coefficients that appropriately compensate for the omitted long-range Lennard-Jones interactions. In terms of the transport properties D and η, the results are also satisfactory. Experimental values for D were only available for 28 compounds. For this subset, 2016H66 provides good results. The largest deviations are found for the compounds MTL, ACA, EDAN, PHL, and PHT, where D is underestimated by about a factor of 2, and the compound 2M2P, where D is overestimated by a factor of 4. Experimental values for η were available for all but one compound, MPE. The comparison with experiment is overall satisfactory. However, the viscosity is systematically and significantly underestimated for the alcohols (except MTL). The case of CHXL is clearly also special, with an experimental value122 of 575 × 10−4 kg m−1 s−1 (i.e., anomalously high), owing to which this compound was omitted from the statistics in Table 6 and its representative point is not displayed in Figure 3. This anomaly is most likely due to the formation of a glassy liquid phase under ambient conditions.204−208 The value of 29 × 10−4 kg m−1 s−1 calculated using 2016H66 underestimates the experimental one by a factor of 20, which suggests that the compound still behaves as a normal liquid within the force-field representation. However, the properties of this compound at 298.15 K may be difficult to capture accurately because this temperature is almost equal to the melting temperature117 of 298.3 K. In terms of the dielectric permittivity ε, a good agreement between calculated and experimental values is observed for the aldehydes, the esters, the amines DEAN, TEAN, and EDAN, the thiols, and the aromatic compounds PHL and PHT. The calculated values for the alcohols, the ketones, the carboxylic acids, the other amines, the amides, and the disulfides are systematically underestimated compared to the experimental ones. The deviation is particularly large for the amides AMD, LNMA, and DAMD, characterized by anomalously high experimental permittivities of116 65.0, 139.5, and 37.8, respectively, owing to which their representative points are not displayed in Figure 3. The corresponding values calculated using 2016H66 are significantly lower, namely 30.3, 49.0, and 20.2, respectively. Note that the representation of the amides and in particular of their dielectric properties within classical force fields is a challenging problem that has attracted considerable attention over the years.95,209−213 The origin of the amide peculiarity has been attributed by different authors to the formation of hydrogenbonded microclusters or chains214,215 stabilized by strong and cooperative216 NH ··· O hydrogen bonds (except for N-disubstituted amides) as well as weaker CH ··· O interactions211 (improper hydrogen bonds210), or to the large magnitude of electronic polarizability effects for this chemical function.212,213 In contrast, the calculated ε values for the ethers and the sulfides are systematically overestimated. For the aromatic compounds, ε is underestimated for BZN, TOL, ANLN, and PYRI. It is overestimated for PHT and, to a lesser extent, PHL. 3.3. Nucleic-Acid Bases. Figure 4 displays the correlation between the experimental and calculated water-to-chloroform transfer free energies ΔGtrn for the five methylated bases MADE, MGUA, MTHY, MCYT, and MURA, considering both the 2016H66 and the 53A6 parameter sets. The corresponding numerical values can be found in SI B (Table B7). Note that for the 2016H66 parameter set, the values predicted using OSP are very close to the ones calculated more accurately using TI. The OSP approach is thus a very efficient and reasonably accurate method to refine the charge set of a solute molecule during force-field parametrization. The values calculated using 2016H66 3840


Article


force field, namely sets63 53A5 and 53A6 (note that set67 54A7 and set68 54A8 are identical to 53A6 for the molecules considered here). The agreement with experiment in terms of all four target properties over the 57 organic molecules is significantly improved by 2016H66. The improvement is observed not only in terms of RMSD but also in terms of AVED, which suggests that systematic errors have been reduced. It should be stressed, however, that this improvement results mainly from an increase in the flexibility of the parametrization procedure. More precisely, the present calibration involves: (1) a slight readjustment of the Lennard-Jones C1/2 6 dispersion parameters (vs adjustment of the C1/2 12 repulsion parameters only in 53A5 and 53A6); (2) a simultaneous adjustment of the partial charges and Lennard-Jones interaction parameters (vs separate adjustments in 53A5 and 53A6). Owing to this increased flexibility, it becomes possible to reproduce simultaneously the four target properties with a higher fitting precision and within a single parameter set (vs two sets for 53A5 and 53A6, reproducing accurately only the pure-liquid or the solvation properties, respectively). As an initial validation of the 2016H66 parameter set, seven additional thermodynamic (cP, κT, αP, and γ), transport (D, η), and dielectric (ε) properties are calculated for the 57 organic molecules in the liquid phase and compared with experimental data. When considering this comparison and its significance in terms of validation, one should keep in mind that it characterizes three features of the force-field representation simultaneously, namely: (1) the accuracy of the underlying physical model (classical description, empirical force-field functional form, and other simulation approximations), i.e., its intrinsic ability to capture the physics of a liquid in all of its relevant aspects; (2) the extent of degeneracy of the parametrization procedure given the selected physical model and target properties, i.e., to which extent the fitting of the target properties against experiment determines a unique solution for the parameter set; (3) the precision of the parametrization procedure, i.e., the accuracy with which the target properties are reproduced by the parameter set. These three aspects are clearly difficult to disentangle in practice. The precision of the parametrization procedure in terms of the target properties can in principle be evaluated by monitoring RMSD and AVED values (see above), but the extent to which residual errors affect the accuracy of the physical model and its level of parametrization degeneracy is difficult to assess. And even assuming a nearly optimal parametrization against the target properties, the two other aspects remain intermingled. More precisely, agreement with experiment in terms of a large number of nontarget properties suggests (but does not prove) that the model physics is correct and that the model parameters are uniquely determined by the selected target properties. Conversely, disagreements in terms of nontarget properties may arise from an intrinsic inaccuracy of the physical model, from a degeneracy of the parametrization, or from both. Here, degeneracy of the parametrization means that the selected target properties do not encompass sufficient information for determining a unique parameter set and that the solution retained, although nearly optimal in terms of the target properties, remains suboptimal with respect to nontarget ones. In practice, the agreement between the 2016H66 results and experimental data for the seven additional properties over the 57 organic molecules appears to be very reasonable. This makes it likely (but not guaranteed) that the physical model employed is appropriate, and that the extent of parametrization degeneracy is limited. Of particular interest is the observation that, for each of the seven properties, the agreement with

Figure 4. Correlation between the experimental and calculated waterto-chloroform transfer free energies (ΔGtrn) for the five methylated nucleic-acid bases. The abbreviations for the bases are given in Table 4. The straight line marks the ideal case of a perfect agreement between simulation and experiment. The experimental values are taken from ref 147. The simulation estimates are obtained using OSP for the 2016H66 force field according to the thermodynamic cycle illustrated in Figure 1 (the reference state involving 2016H66 Lennard-Jones interaction parameters along with the 53A6 charge set), or using TI for either the 2016H66 force field or the 53A6 force field. The values are reported numerically in SI B (Table B7).

are in excellent agreement with the target experimental data. The RMSD and AVED values relative to experiment considering the five bases are 1.7 and 0.8 kJ mol−1, respectively. For comparison purposes, the values calculated using the 53A6 set63 are also shown in Figure 4. These are systematically underestimated, in line with the results reported in a recent study.217

4. CONCLUSIONS In this article, we report on the calibration and validation of a new GROMOS-compatible parameter set 2016H66 for small organic molecules in the condensed phase. The calibration follows the principles and procedures previously adopted in the development of the set77 53A6OXY, the set95 53A6OXY+A, and the set96 53A6OXY+D, and these sets are encompassed within 2016H66. This calibration is based on 62 organic molecules spanning the chemical functions alcohol, ether, aldehyde, ketone, carboxylic acid, ester, amine, amide, thiol, sulfide, and disulfide, as well as aromatic compounds and nucleic-acid bases. For 57 organic compounds, the calibration targets are the experimental pure-liquid properties ρliq and ΔHvap as well as the solvation properties ΔGwat and ΔGche. The final RMSD values for these four quantities over the entire set of compounds are 32.4 kg m−3, 3.5 kJ mol−1, 4.1 kJ mol−1, and 2.1 kJ mol−1, respectively, and the corresponding AVED values are 1.0 kg m−3, 0.2 kJ mol−1, 2.6 kJ mol−1, and 1.0 kJ mol−1, respectively. The low values for the AVED (compared to those of the RMSD) suggest that the residual deviations are predominantly statistical, i.e., exempt of significant systematic errors. For the five nucleic-acid bases, the parametrization is mainly performed by transferring the final 2016H66 parameters from analogous organic compounds. However, the partial charges are slightly readjusted to reproduce the experimental water-to-chloroform transfer free energies ΔGtrn. The final RMSD value for this property over the five bases is 1.7 kJ mol−1, and the corresponding AVED is 0.8 kJ mol−1. These deviations are compared to those observed when using the two most recent major reparametrizations of the GROMOS 3841


Article


the dynamics (time-dependence) of the system, whereas D and η are transport parameters that characterize precisely this dynamics. The dynamical information is thus implicitly added by the physical model. At first sight, the calibration targets ΔGwat and ΔGche may appear irrelevant for defining the component of the model governing pure-liquid properties. In practice, however, this is not the case as the atomic partial charges and Lennard-Jones interaction parameters are jointly determining both the pureliquid and solvation properties. As a result, the calibration targets ΔGwat and ΔGche add physics into the model, which is relevant even for the pure-liquid properties. We believe that the inclusion of this information in 2016H66 was essential to lift the parametrization degeneracy, and that the calibration of a liquid force field solely on the basis of ρliq and ΔHvap would be insufficient even in the context of the sole pureliquid properties. Although it has long been realized that force-field calibration against solvation properties is essential,12,29,43,56,63,128,258−263 the idea that this optimization must be performed simultaneously with the calibration against pure-liquid properties has emerged only recently in the literature.77,264−266 The observation of such a good general agreement in terms of nontarget observables suggests that the choice of parametrization targets is appropriate, which is very encouraging from the point of view of force-field development. The calibration targets ρliq and ΔHvap are easy to parametrize against (relatively short and inexpensive simulations) and available experimentally for many organic compounds.191 The calibration targets ΔGwat and ΔGnps (solvation free energy in nonpolar solvents, e.g., ΔGche) are only slightly more difficult to parametrize against (free-energy calculations) and still available experimentally for many organic compounds.118−120,128,267−269 These four properties are therefore not only appropriate but also relatively easy parametrization targets. At present, the 2016H66 force field is mainly intended for the simulation of small organic molecules, either as pure liquids or as solutes in aqueous or nonpolar solutions, although preliminary 2016H66-compatible parameters for biomolecules are also available134 (see Appendix of SI A). In the long term, the goal of this work is to provide a database of carefully parametrized organic molecules to serve as fragments for: (1) a rigorously fragment-based topology builder (no structure- or QM-based parameter inference as commonly used in available builders270−282) for molecules relevant in organic chemistry, biochemistry, and material sciences; (2) the parametrization of a new GROMOS biomolecular force field; (3) the parametrization of ligand molecules for applications in drug design relying on GROMOS. Work is currently in progress toward these three goals, which are in part similar to those pursued in the development of the all-atom force field37,41,43,44 (OPLS/ AA) for OPLS, of the general AMBER force field22 (GAFF) for AMBER, or of the CHARMM general force field11 (CGenFF) for CHARMM, as well as earlier in the development of the universal force field283,284 (UFF) and the Merck molecular force field285−287 (MMFF). In this respect, it is interesting to compare the overall RMSD values in Tables 5 and 6 to the corresponding values reported in ref 203 (see Table 1 therein, entries without LJPME) for OPLS/AA, GAFF, and CGenFF. The 2016H66 parameter set presents the lowest RMSD values for ρliq, ΔHvap, and γ, comparable RMSD values for αP and κT, and somewhat higher RMSD values for ε. Although RMSD values should not be overinterpreted (as they may be quite sensitive to the molecule

experiment is good for most compounds, with a few exceptions (outliers) typically concerning either a specific chemical function or a specific molecule. Clearly, these special cases should be investigated in detail, as the deviations suggest either a defect in the parametrization procedure, an inaccurate experimental value, or a special experimental property that is not captured by the model physics. In the latter situation, the basic assumptions of the physical model may have to be reconsidered. There exist many different opinions concerning the presumably weakest component of the classical force-field representation, i.e., the element of the model one should reconsider with highest priority. One may cite in particular: (1) the use of simple and uncoupled covalent interaction terms218 (vs inclusion of anharmonicities or/and cross terms); (2) the neglect of longrange van der Waals interactions;201−203,219,220 (3) the isotropic united-atom representation of aliphatic groups156,221−226 (vs anisotropic united-atom or all-atom representations); (4) the application of ad hoc combination rules;12,29,31,219,220,227−232 (5) the neglect of long-range electrostatic interactions;233,234 (6) the implicit (via effective charges) representation of electronic polarizability;12,192,235−247 (7) the atom-centered point-charge representation of electrostatic properties44,242,248−250 (vs inclusion of off-atom charges or/and higher-order multipoles); (8) the use of conformation-independent charges218,250−252 (vs inclusion of conformation-dependent charge fluxes); (9) the decoupling of the van der Waals from the charge model253−256 (i.e., use of QM-derived charges along with a standardized van der Waals atom-type set). One of the main difficulties in assessing the relative importance of these possible shortcomings is that comparisons at a suboptimal level of parametrization can be very misleading,156,203 i.e., two fully optimized parameter sets (one with and one without the change of representation) should always be compared. All the above points except (5) and (9) may potentially affect 2016H66, and we believe that points (2) and (4) are the ones to be reinvestigated with highest priority. Point (5) is addressed here using a reaction-field approximation. Although lattice-sum electrostatics is nowadays very popular, both approaches consider the long-range tail of electrostatic interactions and large differences should not be expected, at least in the context of pure organic liquids. Importantly, 2016H66 also addresses point (9). In many current small-molecule force fields,11,22,43,44,253,254,257 QM-derived (thus also in principle dependent or the choice of a reference structure and environment) atomic charges are used along within standardized van der Waals atom-type sets, a situation that is somewhat reminiscent of the separate optimization of van der Waals parameters and charges in the design of the 53A6 set63, and may result in suboptimal combinations.77 In contrast, in 2016H66, the van der Waals and charge models are refined simultaneously and consistently for all atoms within the different functional groups. The general agreement observed for most compounds within 2016H66 in terms of the seven nontarget properties is in fact truly remarkable. It suggests that, given a physical model based on classical mechanics with an empirical (and relatively simple) interaction function, agreement with experiment solely in terms of the liquid properties ρliq and ΔHvap and the solvation properties ΔGwat and ΔGche is already sufficient to reach an appropriate description of the liquid in terms of a wide spectrum of its physical properties. As a simple illustration of this fact, we note that all four target properties are purely thermodynamic parameters that do not encompass any information on 3842


Article


Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (8) MacKerell, A. D.; Banavali, N.; Foloppe, N. Development and current status of the CHARMM force field for nucleic acids. Biopolymers 2000, 56, 257−265. (9) Brooks, B. R.; Brooks, C. L.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (10) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; McKerell, A. D.; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B 2010, 114, 7830− 7843. (11) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; McKerell, A. D., Jr CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671−690. (12) Baker, C. M.; Lopes, P. E. M.; Zhu, X.; Roux, B.; McKerell, A. D. Accurate calculation of hydration free energies using pair-specific Lennard-Jones parameters in the CHARMM Drude polarizable force field. J. Chem. Theory Comput. 2010, 6, 1181−1198. (13) Guvench, P.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W.; MacKerell, A. D., Jr CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate−protein modeling. J. Chem. Theory Comput. 2011, 7, 3162−3180. (14) Zhu, Z.; Lopes, P. E. M.; McKerell, A. D., Jr Recent developments and applications of the CHARMM force fields. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 167−185. (15) Vanommeslaeghe, K.; MacKerell, A. D., Jr CHARMM additive and polarizable force fields for biophysics and computer-aided drug design. Biochim. Biophys. Acta 2015, 1850, 861−871. (16) Weiner, P. K.; Kollman, P. A. AMBERassisted model building with energy refinementa general program for modeling molecules and their interactions. J. Comput. Chem. 1981, 2, 287−303. (17) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E.; DeBolt, S.; Fergusson, D.; Seibel, G.; Kollman, P. AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 1995, 91, 1−41. (18) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (19) Senderowitz, H.; Parish, C.; Still, W. C. Carbohydrates: united atom AMBER* parametrization of pyranoses and simulations yelding anomeric free energies. J. Am. Chem. Soc. 1996, 118, 2078−2086. (20) Senderowitz, H.; Still, W. C. A quantum mechanically derived all-atom force field for pyranose oligosaccharides. AMBER* parameters and free energy simulations. J. Org. Chem. 1997, 62, 1427−1438. (21) Glennon, T. M.; Merz, K. M., Jr A carbohydrate force field for AMBER and its application to the study of saccharide to surface adsorption. J. Mol. Struct.: THEOCHEM 1997, 395, 157−171. (22) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general AMBER force field. J. Comput. Chem. 2004, 25, 1157−1174.

set considered and to the presence of outliers), the united-atom 2016H66 parameter set certainly stands the comparison against the three all-atom force fields. Work is in progress to further extend the scope of 2016H66 so as to encompass new organic molecules involving more diverse chemical functions. However, because the goals stated above are somewhat distinct from those pursued in the development of the GROMOS biomolecular force field, a new naming convention (2016H66) has been adopted to label the first parameter set of this new force-field line (H-series).

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00187. The specification of the 2016H66 parameter set, i.e., the content of the GROMOS ifp and (partially only) mtb files (PDF) More information on the organic molecules considered, including the experimental values of their properties as well as the detailed simulation results concerning these properties (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. +55 21 3938 7132. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Wilfred van Gunsteren for numerous discussions and fruitful suggestions concerning this work. Thanks are also due to Andreas Eichenberger for developing preliminary torsional parameters for polypeptides (see Appendix in SI A), to Victor Holanda Rusu for helping with the GROMACS files, to Niels Hansen for helping with the forcefield files for cyclodextrin and DMSO2, and to Wojtek Plazinski for helping with the carbohydrate. Financial support from the Swiss National Science Foundation (grant nos. 21-132739 and 21-138020) is also gratefully acknowledged.

■

REFERENCES

(1) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (2) van Gunsteren, W. F.; Berendsen, H. J. C. Computer simulation of molecular dynamics: methodology, applications and perspectives in chemistry. Angew. Chem., Int. Ed. 1990, 29, 992−1023. (3) van Gunsteren, W. F.; Bakowies, D.; Baron, R.; Chandrasekhar, I; Christen, M.; Daura, X.; Gee, P.; Geerke, D. P.; Glättli, A.; Hünenberger, P. H.; Kastenholz, M. A.; Oostenbrink, C.; Schenk, M.; Trzesniak, D.; van der Vegt, N. F. A.; Yu, H. B. Biomolecular modelling: goals, problems, perspectives. Angew. Chem., Int. Ed. 2006, 45, 4064−4092. (4) Berendsen, H. J. C. Simulating the Physical World; Cambridge University Press: Cambridge, U.K., 2007. (5) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: a program for macromolecular energy, minimization and dynamics calculations. J. Comput. Chem. 1983, 4, 187−217. (6) MacKerell, A. D.; Wiorkiewicz-kuczera, J.; Karplus, M. An allatom empirical energy function for the simulation of nucleic-acids. J. Am. Chem. Soc. 1995, 117, 11946−11975. (7) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; 3843


Article


energies using the next generation OPLS force field. J. Chem. Theory Comput. 2012, 8, 2553−2558. (44) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A. OPLS3: a force field providing broad coverage of drug-like small molecules and proteins. J. Chem. Theory Comput. 2016, 12, 281−296. (45) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (46) Chen, B.; Potoff, J. J.; Siepmann, J. I. Monte Carlo calculations for alcohols and their mixtures with alkanes. Transferable potentials for phase equilibria. 5. United-atom description of primary, secondary, and tertiary alcohols. J. Phys. Chem. B 2001, 105, 3093−3104. (47) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable potentials for phase equilibria. 6. United-atom description for ethers, glycols, ketones, and aldehydes. J. Phys. Chem. B 2004, 108, 17596−17605. (48) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 7. Primary, secondary, and tertiary amines, nitroalkanes and nitrobenzene, nitriles, amides, pyridine, and pyrimidine. J. Phys. Chem. B 2005, 109, 18974−18982. (49) Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 8. United-atom description for thiols, sulfides, disulfides, and thiophene. J. Phys. Chem. B 2005, 109, 24100−24107. (50) Kamath, G.; Robinson, J.; Potoff, J. J. Application of TraPPE-UA force field for determination of vapor−liquid equilibria of carboxylate esters. Fluid Phase Equilib. 2006, 240, 46−55. (51) Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 10. Explicit-hydrogen description of substituted benzenes and polycyclic aromatic compounds. J. Phys. Chem. B 2012, 117, 273− 288. (52) Eggimann, B. L.; Sunnarborg, A. J.; Stern, H. D.; Bliss, A. P.; Siepmann, J. I. An online parameter and property database for the TraPPE force field. Mol. Simul. 2014, 40, 101−105. (53) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J. P. M. A consistent empirical potential for water−protein interactions. Biopolymers 1984, 23, 1513−1518. (54) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation (GROMOS) Library Manual; BIOMOS: Groningen, the Netherlands, 1987. (55) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hünenberger, P. H.; Krüger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; Verlag der Fachvereine: Zürich, Switzerland, 1996. (56) Daura, X.; Mark, A. E.; van Gunsteren, W. F. Parametrization of aliphatic CHn united atoms of GROMOS96 force field. J. Comput. Chem. 1998, 19, 535−547. (57) van Gunsteren, W. F.; Daura, X.; Mark, A. E. In Encyclopedia of Computational Chemistry; Schleyer, P., Ed.; John Wiley & Sons: Chichester, U.K., 1998; Vol. 2, pp 1211−1216. (58) Scott, W. R. P.; Hünenberger, P. H.; Tironi, I. G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Krüger, P.; van Gunsteren, W. F. The GROMOS biomolecular simulation program package. J. Phys. Chem. A 1999, 103, 3596−3607. (59) Schuler, L. D.; van Gunsteren, W. F. On the choice of dihedral angle potential energy functions for n-alkanes. Mol. Simul. 2000, 25, 301−319. (60) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase. J. Comput. Chem. 2001, 22, 1205−1218. (61) Chandrasekhar, I.; Kastenholz, M. A.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman, D. P.; van Gunsteren, W. F. A consistent potential energy parameter set for lipids: dipalmitoylphosphatidylcholine as a benchmark of the GROMOS96 45A3 force field. Eur. Biophys. J. 2003, 32, 67−77. (62) Chandrasekhar, I.; Oostenbrink, C.; van Gunsteren, W. F. Simulating the physiological phase of hydrated DPPC bilayers: the ester moiety. Soft Mater. 2004, 2, 27−45.

(23) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The AMBER biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668−1688. (24) Hornak, V.; Abel, R.; Okur, A.; Stockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple AMBER force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Bioinf. 2006, 65, 712−725. (25) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzálezOuteriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: a generalizable biomolecular force field. Carbohydrates. J. Comput. Chem. 2008, 29, 622−655. (26) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the AMBER ff99SB protein force field. Proteins 2010, 78, 1950− 1958. (27) Chapman, D. E.; Steck, J. K.; Nernberg, P. S. Optimizing protein-protein van der Waals interactions for the AMBER ff9x/ff12 force field. J. Chem. Theory Comput. 2014, 10, 273−281. (28) Dickson, C. J.; Madej, B. D.; Skjevik, A. A.; Betz, R. M.; Teigen, K.; Gould, I. R.; Walker, R. C. Lipid14: the AMBER lipid force field. J. Chem. Theory Comput. 2014, 10, 865−879. (29) Jämbeck, J. P. M.; Lyubartsev, A. P. Update of the general AMBER force field for small solutes with an emphasis on free energies of hydration. J. Phys. Chem. B 2014, 118, 3793−3804. (30) Nikitin, A. M.; Milchevskiy, Y. V.; Lyubartsev, A. P. A new AMBER-compatible force field parameter set for alkanes. J. Mol. Model. 2014, 20, 2143/1−2143/10. (31) Nikitin, A.; Milchevskiy, Y.; Lyubartsev, A. AMBER-II: new combining rules and force field for perfluoroalkanes. J. Phys. Chem. B 2015, 119, 14563−14573. (32) Jorgensen, W. L. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J. Am. Chem. Soc. 1981, 103, 335−340. (33) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized intermolecular potential functions for liquid hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (34) Jorgensen, W. L. Optimized intermolecular potential functions for liquid alcohols. J. Phys. Chem. 1986, 90, 1276−1284. (35) Jorgensen, W. L.; Tirado-Rives, J. The OPLS potential functions for proteins. Energy minimizations for crystals of cyclic peptides and crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (36) Pranata, J.; Wierschke, S. G.; Jorgensen, W. L. OPLS potential functions for nucleotide bases. Relative association constants of hydrogen-bonded base pairs in chloroform. J. Am. Chem. Soc. 1991, 113, 2810−2819. (37) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (38) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. OPLS all-atom force field for carbohydrates. J. Comput. Chem. 1997, 18, 1955−1970. (39) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B 2001, 105, 6474−6487. (40) Kony, D.; Damm, W.; Stoll, S.; van Gunsteren, W. F. An improved OPLS-AA force field for carbohydrates. J. Comput. Chem. 2002, 23, 1416−1429. (41) Jorgensen, W. L.; Tirado-Rives, J. Potential energy functions for atomic-level simulations of water and organic and biomolecular systems. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6665−6670. (42) Jorgensen, W. L.; Schyman, P. Treatment of halogen bonding in the OPLS-AA force field: application to potent anti-HIV agents. J. Chem. Theory Comput. 2012, 8, 3895−3901. (43) Shivakumar, D.; Harder, E.; Damm, W.; Friesner, R. A.; Sherman, W. Improving the prediction of absolute solvation free 3844


Article


electrolyte solutions, lipid bilayers, and proteins. J. Chem. Theory Comput. 2013, 9, 1247−1264. (84) Eichenberger, A. P.; Allison, J. R.; Dolenc, J.; Geerke, D. P.; Horta, B. A. C.; Meier, K.; Oostenbrink, C.; Schmid, N.; Steiner, D.; Wang, D.; van Gunsteren, W. F. The GROMOS++ software for the analysis of biomolecular simulation trajectories. J. Chem. Theory Comput. 2011, 7, 3379−3390. (85) Schmid, N.; Allison, J. R.; Dolenc, J.; Eichenberger, A. P.; Kunz, A. P. E.; van Gunsteren, W. F. Biomolecular structure refinement using the GROMOS simulation software. J. Biomol. NMR 2011, 51, 265− 281. (86) Riniker, S.; Christ, C. D.; Hansen, H. S.; Hünenberger, P. H.; Oostenbrink, C.; Steiner, D.; van Gunsteren, W. F. Calculation of relative free energies for ligand−protein binding, solvation and conformational transitions using the GROMOS biomolecular simulation software. J. Phys. Chem. B 2011, 115, 13570−13577. (87) Kunz, A.-P. E.; Allison, J. R.; Geerke, D. P.; Horta, B. A. C.; Hünenberger, P. H.; Riniker, S.; Schmid, N.; van Gunsteren, W. F. New functionalities in the GROMOS biomolecular simulation software. J. Comput. Chem. 2012, 33, 340−353. (88) Schmid, N.; Christ, C. D.; Christen, M.; Eichenberger, A. P.; van Gunsteren, W. F. Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput. Phys. Commun. 2012, 183, 890−903. (89) Poger, D.; van Gunsteren, W. F.; Mark, A. E. A new force field for simulating phosphatidylcholine bilayers. J. Comput. Chem. 2010, 31, 1117−1125. (90) Reif, M. M.; Hünenberger, P. H. Computation of methodologyindependent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones parameter sets for the alkali and halide ions in water. J. Chem. Phys. 2011, 134, 144104/1−144104/25. (91) Pol-Fachin, L.; Rusu, V. H.; Verli, H.; Lins, R. D. GROMOS 53A6GLY C, an improved force field for hexopyranose-based carbohydrates. J. Chem. Theory Comput. 2012, 8, 4681−4690. (92) Pol-Fachin, L.; Verli, H.; Lins, R. D. Extension and validation of the GROMOS 53A6GLYC parameter set for glycoproteins. J. Comput. Chem. 2014, 35, 2087−2095. (93) Hansen, H. S.; Hünenberger, P. H. A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers and glycosidic linkage conformers. J. Comput. Chem. 2011, 32, 998−1032. (94) Plazinski, W.; Lonardi, A.; Hünenberger, P. H. Revision of the GROMOS 56A6CARBO force field: improving the description of ringconformational equilibria in hexopyranose-based carbohydrates chains. J. Comput. Chem. 2016, 37, 354−365. (95) Horta, B. A. C.; Lin, Z.; Huang, W.; Riniker, S.; van Gunsteren, W. F.; Hünenberger, P. H. Reoptimized interaction parameters for the peptide-backbone model compound N-methylacetamide in the GROMOS force field: influence on the folding properties of two beta-peptides in methanol. J. Comput. Chem. 2012, 33, 1907−1917. (96) Fuchs, P. F. J.; Hansen, H. S.; Hünenberger, P. H.; Horta, B. A. C. A GROMOS parameter set for vicinal diether functions: properties polyethyleneoxide and polyethyleneglycol. J. Chem. Theory Comput. 2012, 8, 3943−3963. (97) Horta, B. A. C.; Hünenberger, P. H. Enantiomeric segregation in the gel phase of lipid bilayers. J. Am. Chem. Soc. 2011, 133, 8464− 8466. (98) Riniker, S.; Horta, B. A. C.; Thijssen, B.; Gupta, S.; van Gunsteren, W. F.; Hünenberger, P. H. Temperature dependence of the dielectric permittivity of acetic acid, propionic acid and their methyl esters: a molecular dynamics simulation study. ChemPhysChem 2012, 13, 1182−1190. (99) Rossi, G.; Fuchs, P. F. J.; Barnoud, J.; Monticelli, L. A coarsegrained MARTINI model of polyethylene glycol and of polyethylene alkyl ether surfactants. J. Phys. Chem. B 2012, 116, 14353−14362. (100) Hansen, N.; Hünenberger, P. H.; van Gunsteren, W. F. Efficient combination of environment change and alchemical perturbation within the enveloping distribution sampling (EDS)

(63) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (64) Soares, T. A.; Hünenberger, P. H.; Kastenholz, M. A.; Kräutler, V.; Lenz, T.; Lins, R. D.; Oostenbrink, C.; van Gunsteren, W. F. An improved nucleic-acid parameter set for the GROMOS force field. J. Comput. Chem. 2005, 26, 725−737. (65) Lins, R. D.; Hünenberger, P. H. A new GROMOS force field for hexopyranose-based carbohydrates. J. Comput. Chem. 2005, 26, 1400− 1412. (66) Winger, M.; de Vries, A. H.; van Gunsteren, W. F. Force-field dependence of the conformational properties of α,ω-dimethoxypolyethylene glycol. Mol. Phys. 2009, 107, 1313−1321. (67) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F. Definition and testing of the GROMOS force-field versions 54A7 and 54B7. Eur. Biophys. J. 2011, 40, 843−856. (68) Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. New interaction parameters for charged amino acid side chains in the GROMOS force field. J. Chem. Theory Comput. 2012, 8, 3705−3723. (69) van Gunsteren, W. F. The GROMOS Sof tware for Biomolecular Simulation. http://www.gromos.net (accessed May 5, 2011). (70) Levitt, M.; Lifson, S. Refinement of protein conformations using a macromolecular energy minimization procedure. J. Mol. Biol. 1969, 46, 269−279. (71) Lennard-Jones, J. E. The equation of state of gases and critical phenomena. Physica 1937, 4, 941−956. (72) Hagler, A. T.; Huler, E.; Lifson, S. Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals. J. Am. Chem. Soc. 1974, 96, 5319− 5327. (73) Lifson, S.; Hagler, A. T.; Dauber, P. Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 1. Carboxylic acids, amides, and the C=O···H hydrogen bonds. J. Am. Chem. Soc. 1979, 101, 5111−5121. (74) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, the Netherlands, 1981; pp 331−342. (75) Barker, J. A.; Watts, R. O. Monte Carlo studies of the dielectric properties of water-like models. Mol. Phys. 1973, 26, 789−792. (76) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 1995, 102, 5451−5459. (77) Horta, B. A. C.; Fuchs, P. F. J.; van Gunsteren, W. F.; Hü nenberger, P. H. New interaction parameters for oxygen compounds in the GROMOS force field: improved pure-liquid and solvation properties for alcohols, ethers, aldehydes, ketones, carboxylic acids and esters. J. Chem. Theory Comput. 2011, 7, 1016−1031. (78) Villa, A.; Mark, A. E. Calculation of the free energy of solvation for neutral analogues of amino acid side chains. J. Comput. Chem. 2002, 23, 548−553. (79) Soares, T. A.; Daura, X.; Oostenbrink, C.; Smith, L. J.; van Gunsteren, W. F. Validation of the GROMOS force-field parameter set 45A3 against nuclear magnetic resonance data of hen egg lysozyme. J. Biomol. NMR 2004, 30, 407−422. (80) Oostenbrink, C.; Soares, T. A.; van der Vegt, N. F. A.; van Gunsteren, W. F. Validation of the 53A6 GROMOS force field. Eur. Biophys. J. 2005, 34, 273−284. (81) Geerke, D. P.; van Gunsteren, W. F. Force field evaluation for biomolecular simulation: free enthalpies of solvation of polar and apolar compounds in various solvents. ChemPhysChem 2006, 7, 671− 678. (82) Huang, W.; Lin, Z.; van Gunsteren, W. F. Validation of the GROMOS 54A7 force field with respect to β-peptide folding. J. Chem. Theory Comput. 2011, 7, 1237−1243. (83) Reif, M. M.; Winger, M.; Oostenbrink, C. Testing of the GROMOS force-field parameter set 54A8: structural properties of 3845


Article

Journal of Chemical Theory and Computation scheme: twin-system EDS and application to the determination of octanol−water partition coefficients. J. Chem. Theory Comput. 2013, 9, 1334−1346. (101) Laner, M.; Horta, B. A. C.; Hünenberger, P. H. Phasetransition properties of glycerol-monopalmitate lipid bilayers investigated by molecular dynamics simulation: influence of the system size and force-field parameters. Mol. Simul. 2013, 39, 563−583. (102) Pechlaner, M.; Sigel, R. K. O.; van Gunsteren, W. F.; Dolenc, J. Structure and conformational dynamics of the domain 5 RNA hairpin of a bacterial group II intron revealed by solution nuclear magnetic resonance and molecular dynamics simulations. Biochemistry 2013, 52, 7099−7113. (103) Huang, W.; Lin, Z.; van Gunsteren, W. F. Use of enveloping distribution sampling to evaluate important characteristics of biomolecular force fields. J. Phys. Chem. B 2014, 118, 6424−6430. (104) Benson, S. P.; Pleiss, J. Molecular dynamics simulations of selfemulsifying drug-delivery systems (SEDDS): influence of excipients on droplet nanostructure and drug localization. Langmuir 2014, 30, 8471−8480. (105) Tang, X.; Koenig, P. H.; Larson, R. G. Molecular dynamics simulations of sodium dodecyl sulfate micelles in water. The effect of the force field. J. Phys. Chem. B 2014, 118, 3864−3880. (106) Tang, X.; Huston, K. J.; Larson, R. G. Molecular dynamics simulations of structure−property relationships of Tween 80 surfactants in water and at interfaces. J. Phys. Chem. B 2014, 118, 12907−12918. (107) Laner, M.; Horta, B. A. C.; Hünenberger, P. H. Effect of the cosolutes trehalose and methanol on the equilibrium and phasetransition properties of glycerol-monopalmitate lipid bilayers investigated using molecular dynamics simulations. Eur. Biophys. J. 2014, 43, 517−544. (108) Laner, M.; Horta, B. A. C.; Hünenberger, P. H. Long-timescale motions in glycerol-monopalmitate lipid bilayers investigated using molecular dynamics simulation. J. Mol. Graphics Modell. 2015, 55, 48− 64. (109) Laner, M.; Hünenberger, P. H. Effect of methanol on the phase-transition properties of glycerol-monopalmitate lipid bilayers investigated using molecular dynamics simulations: in quest of the biphasic effect. J. Mol. Graphics Modell. 2015, 55, 85−104. (110) Laner, M.; Hünenberger, P. H. Phase-transition properties of gycerol-dipalmitate lipid bilayers investigated using molecular dynamics simulations. J. Mol. Graphics Modell. 2015, 59, 136−147. (111) Huston, K. J.; Larson, R. G. Reversible and irreversible adsorption energetics of poly(ethyleneglycol) and sorbitan poly(ethoxylate) at a water/alkane interface. Langmuir 2015, 31, 7503− 7511. (112) Bieler, N. S.; Tschopp, J. P.; Hünenberger, P. H. Multistate λ-local-elevation umbrella-sampling (MS-λ-LEUS): method and application to the complexation of cations by crown ethers. J. Chem. Theory Comput. 2015, 11, 2575−2588. (113) Stanzione, F.; Jayaraman, A. Hybrid atomistic and coarsegrained molecular dynamics simulations of polyethylene glycol (PEG) in explicit water. J. Phys. Chem. B 2015, 120, 4160−4173. (114) Reif, M. M.; Hünenberger, P. H. Origin of asymmetric solvation effects for ions in water and organic solvents investigated using molecular dynamics simulations: the Swain acity-basity scale revisited. J. Phys. Chem. B 2016, in press, DOI: 10.1021/ acs.jpcb.6b02156. (115) IUPAC. Compendium of Chemical Terminology (the “Gold Book”), 2nd ed.; Blackwell Scientific Publications: Oxford, U.K., 1997. (116) Wohlfahrt, C. Landolt-Börnstein. Numerical Data and Functional Relationships in Science and Technology; Madelung, O., Eds.; Springer: Berlin, Germany, 1991; Vol. 6, pp 5−228. (117) Riddick, J. A.; Bunger, W. B.; Sakano, T. K. Organic Solvents, Physical Properties and Methods of Purification; John Wiley & Sons: New York, 1986. (118) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Universal quantum mechanical model for solvation free energies based on gasphase geometries. J. Phys. Chem. B 1998, 102, 3257−3271.

(119) Li, J.; Zhu, T.; Hawkins, G. D.; Winget, P.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G. Extension of the platform of applicability of the SM5.42R universal solvation model. Theor. Chem. Acc. 1999, 103, 9−63. (120) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Generalized Born solvation model SM12. J. Chem. Theory Comput. 2013, 9, 609− 620. (121) Yaws, C. L. The Yaws Handbook of Physical Properties for Hydrocarbons and Chemicals; Gulf Professional Publishing: Oxford, U.K., 2015. (122) Haynes, W. M. CRC Handbook of Chemistry and Physics, 96th ed.; CRC Press/Taylor and Francis: Boca Raton, FL, 2015. (123) Gallicchio, E.; Zhang, L. Y.; Levy, R. M. The SGB/NP hydration free energy model based on the surface generalized born solvent reaction field and novel nonpolar hydration free energy estimators. J. Comput. Chem. 2001, 23, 517−529. (124) Wu, J.; Liu, Z.; Bi, S.; Meng, X. Viscosity of saturated liquid dimethyl ether from (227 to 343) K. J. Chem. Eng. Data 2003, 48, 426−429. (125) InfoTherm release 3.3, 2015.11. http://www.infotherm.com (accessed 2015). (126) Pedley, J. B.; Naylor, R. D.; Kirby, S. P. Thermodynamical Data of Organic Compounds; Chapman and Hall: London, 1986. (127) Majer, V.; Svoboda, V. Enthalpies of Vaporization of Organic Compounds: A Critical Review and Data Compilation; Blackwell Scientific Publications: Oxford, U.K., 1985. (128) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. Small molecule hydration free energies in explicit solvent: an extensive test of fixed-charge atomistic simulations. J. Chem. Theory Comput. 2009, 5, 350−358. (129) Radzicka, A.; Wolfenden, R. Comparing the polarities of the amino acids: side-chain distribution coefficients between the vapor phase, cyclohexane, 1-octanol, and neutral aqueous solution. Biochemistry 1988, 27, 1664−1670. (130) Lemire, R. J.; Sears, P. G. In Topics in Current Chemistry; Boschke, F. L., Ed.; Springer: Berlin, Germany, 1978; Vol. 74, pp 45−91. (131) MacKerell, A. D.; Shim, J. H.; Anisimov, V. M. Re-evaluation of the reported experimental values of the heat of vaporization of N-methylacetamide. J. Chem. Theory Comput. 2008, 4, 1307−1312. (132) Gagarin, S. G.; Gyul’maliev, A. M. Heat of vaporization of methyl-substituted phenol derivatives. Coke Chem. 2007, 50, 232−236. (133) Liessmann, G.; Schmidt, W.; Reiffarth, S. Data Compilation of the Sächsische Olefinwerke Boehlen, DETHERM database, Germany; 1995; p 48. (134) Hünenberger, P. H. Computer Simulation of Molecular Systems/ GROMOS. http://www.csms.ethz.ch/files_and_links/GROMOS.html (accessed June 1, 2016). (135) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (136) Tironi, I. G.; van Gunsteren, W. F. A molecular dynamics study of chloroform. Mol. Phys. 1994, 83, 381−403. (137) Geerke, D. P.; Oostenbrink, C.; van der Vegt, N. F. A.; van Gunsteren, W. F. An effective force field for molecular dynamics simulations of dimethyl sulfoxide and dimethyl sulfoxide−water mixtures. J. Phys. Chem. B 2004, 108, 1436−1445. (138) Walser, R.; Mark, A. E.; van Gunsteren, W. F.; Lauterbach, M.; Wipff, G. The effect of force-field parameters on properties of liquids: parametrization of a simple three-site model for methanol. J. Chem. Phys. 2000, 112, 10450−10459. (139) Tironi, I. G.; Fontana, P.; van Gunsteren, W. F. A molecular dynamics simulation study of liquid carbon tetrachloride. Mol. Simul. 1996, 18, 1−11. (140) Fioroni, M.; Burger, K.; Mark, A. E.; Roccatano, D. A new 2,2,2-trifluoroethanol model for molecular dynamics simulations. J. Phys. Chem. B 2000, 104, 12347−12354. (141) Smith, L. J.; Berendsen, H. J. C.; van Gunsteren, W. F. Computer simulation of urea−water mixtures: a test of force-field 3846


Article

Journal of Chemical Theory and Computation parameters for use in biomolecular simulation. J. Phys. Chem. B 2004, 108, 1065−1071. (142) Hansen, N.; Kraus, P.; Sassmannhausen, H.; Timmerscheidt, T.; van Gunsteren, W. F. An effective force field for molecular dynamics simulations of dimethyl sulfone. Mol. Phys. 2011, 109, 2593−2605. (143) Gebhardt, J.; Hansen, N. Calculation of binding affinities for linear alcohols to α-cyclodextrin by twin-system enveloping distribution sampling simulations. Fluid Phase Equilib. 2016, 422, 1−17. (144) Horta, B. A. C.; de Vries, A. H.; Hünenberger, P. H. Simulating the transition between gel and liquid-crystal phases of lipid bilayers: dependence of the transition temperature on the hydration level. J. Chem. Theory Comput. 2010, 6, 2488−2500. (145) Marcus, Y. The Properties of Solvents, 1st ed.; John Wiley & Sons: Chichester, U.K., 1998. (146) Slater, J. C.; Kirkwood, J. G. The van der Waals forces in gases. Phys. Rev. 1931, 37, 682−697. (147) Cullis, P. M.; Wolfenden, R. Affinities of nucleic acid bases for solvent water. Biochemistry 1981, 20, 3024−3028. (148) Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300−313. (149) Zwanzig, R. W. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J. Chem. Phys. 1954, 22, 1420−1426. (150) Liu, H.; Mark, A. E.; van Gunsteren, W. F. Estimating the relative free energy of different molecular states with respect to a single reference state. J. Phys. Chem. 1996, 100, 9485−9494. (151) van Gunsteren, W. F.; Berendsen, H. J. C.; Rullmann, J. A. C. Stochastic dynamics for molecules with constraints. Brownian dynamics of n-alkanes. Mol. Phys. 1981, 44, 69−95. (152) Guàrdia, E.; Padró, J. A. Generalized Langevin dynamics simulation of interacting particles. J. Chem. Phys. 1985, 83, 1917−1920. (153) van Gunsteren, W. F.; Berendsen, H. J. C. A leap-frog algorithm for stochastic dynamics. Mol. Simul. 1988, 1, 173−185. (154) Yun-yu, S.; Lu, W.; van Gunsteren, W. F. On the approximation of solvent effects on the conformation and dynamics of cyclosporin A by stochastic dynamics simulation techniques. Mol. Simul. 1988, 1, 369−383. (155) van Gunsteren, W. F. In Computer Simulation of Biomolecular Systems, Theoretical and Experimental Applications; van Gunsteren, W. F., Weiner, P. K., Wilkinson, A. J., Eds.; ESCOM Science Publishers, B.V.: Leiden, the Netherlands, 1993; Vol. 2, pp 3−36. (156) Oostenbrink, C.; Juchli, D.; van Gunsteren, W. F. Amine hydration: a united-atom force-field solution. ChemPhysChem 2005, 6, 1800−1804. (157) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; di Nola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (158) Gavish, B.; Gratton, E.; Hardy, C. J. Adiabatic compressibility of globular proteins. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 750−754. (159) Aicart, E.; Tardajos, G.; Diaz Pena, M. Isothermal compressibility of cyclohexane + n-hexane, cyclohexane + n-heptane, cyclohexane + n-octane, and cyclohexane + n-nonane. J. Chem. Eng. Data 1980, 25, 140−145. (160) Hockney, R. W. The potential calculation and some applications. In Methods in Computational Physics: Academic Press, New York, 1970; Vol. 9, pp 135−211. (161) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (162) Berendsen, H. J. C.; van Gunsteren, W. F.; Zwinderman, H. R. J.; Geurtsen, R. G. Simulations of proteins in water. Ann. N.Y. Acad. Sci. 1986, 482, 269−286. (163) Heinz, T. N.; van Gunsteren, W. F.; Hünenberger, P. H. Comparison of four methods to compute the dielectric permittivity of liquids from molecular dynamics simulations. J. Chem. Phys. 2001, 115, 1125−1136.

(164) Staudhammer, P.; Seyer, W. F. The dielectric constant of cisand trans-decahydronaphthalene and cyclohexane as a function of temperature and frequency. J. Am. Chem. Soc. 1958, 80, 6491−6494. (165) Heinz, T. N.; Hünenberger, P. H. Combining the lattice-sum and reaction-field approaches for evaluating long-range electrostatic interactions in molecular simulations. J. Chem. Phys. 2005, 123, 034107/1−034107/19. (166) Christen, M.; Hünenberger, P. H.; Bakowies, D.; Baron, R.; Bürgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenholz, M. A.; Kräutler, V.; Oostenbrink, C.; Peter, C.; Trzesniak, D.; van Gunsteren, W. F. The GROMOS software for biomolecular simulation: GROMOS05. J. Comput. Chem. 2005, 26, 1719−1751. (167) IUPAC. Quantities, Units and Symbols in Physical Chemistry (the “Green Book”), 2nd ed.; Blackwell Scientific Publications: Oxford, U.K., 1993. (168) Beutler, T. C.; Mark, A. E.; van Schaik, R.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529−539. (169) Harvey, S. C.; Tan, R. K.-Z.; Cheatham, T. E., III The flying ice cube: velocity rescaling in molecular dynamics leads to violation of energy equipartition. J. Comput. Chem. 1998, 19, 726−740. (170) Chiu, S.-W.; Clark, M.; Subramaniam, S.; Jakobsson, E. Collective motion artifacts arising in long-duration molecular dynamics simulations. J. Comput. Chem. 2000, 21, 121−131. (171) Hünenberger, P. H. Thermostat algorithms for molecular dynamics simulations. Adv. Polym. Sci. 2005, 173, 105−149. (172) Lingenheil, M.; Denschlag, R.; Reichold, R.; Tavan, P. The “hot-solvent/cold-solute” problem revisited. J. Chem. Theory Comput. 2008, 4, 1293−1306. (173) Eastwood, M. P.; Stafford, K. A.; Lippert, R. A.; Jensen, M.Ø.; Maragakis, P.; Predescu, C.; Dror, R. O.; Shaw, D. E. Equipartition and the calculation of temperature in biomolecular simulations. J. Chem. Theory Comput. 2010, 6, 2045−2058. (174) Basconi, J. E.; Shirts, M. R. Effects of temperature control algorithms on transport properties and kinetics in molecular dynamics simulations. J. Chem. Theory Comput. 2013, 9, 2887−2899. (175) Uosaki, Y.; Motoki, T.; Hamaguchi, T.; Moriyoshi, T. Excess molar volumes of binary mixtures of 1,3-dimethylimidazolidin-2-one with an alkan-1-ol at the temperatures 283.15 K, 298.15 K, and 313.15 K. J. Chem. Thermodyn. 2007, 39, 810−816. (176) Oswal, S. L.; Oswal, P.; Gardas, R. L.; Patel, S. G.; Shinde, R. G. Acoustic, volumetric, compressibility and refractivity properties and reduction parameters for the ERAS and Flory models of some homologous series of amines from 298.15 to 328.15 K. Fluid Phase Equilib. 2004, 216, 33−45. (177) Matsuo, S.; Makita, T. Volumetric properties of 1-alkanols at temperatures in the range 298−348 K and pressures up to 40 MPa. Int. J. Thermophys. 1989, 10, 885−897. (178) Sahli, B. P.; Gager, H.; Richard, A. J. Ultracentrifugal studies of the isothermal compressibilities of organic alcohols and alkanes. Correlation with surface tension. J. Chem. Thermodyn. 1976, 8, 179− 188. (179) Hamad, E. Z.; Mansoori, G. A.; Matteoli, E.; Lepori, L. Relations among concentration fluctuation integrals in mixtures (theory and experiments). Z. Phys. Chem. (Muenchen, Ger.) 1989, 162, 27−45. (180) Smith, T. E.; Bonner, R. F. Acetaldehyde, propionaldehyde, and n-butyraldehyde: some physical properties. Ind. Eng. Chem. 1951, 43, 1169−1173. (181) Mohammad, A. A.; Alkhaldi, K. H. A. E.; AlTuwaim, M. S.; AlJimaz, A. S. Viscosity and surface tension of binary systems of N, Ndimethylformamide with alkan-1-ols at different temperatures. J. Chem. Thermodyn. 2013, 56, 106−113. (182) Hopfe, D. Data Compilation of FIZCHEMIE, Fiz Chemie Berlin/Wiley - VCH, Germany; 1990; p 20. (183) Coto, B.; Cabañas, A.; Pando, C.; Menduiña, C.; Rubio, R. G.; Renuncio, J. A. R. Bulk and surface properties of the highly non-ideal 3847


Article

Journal of Chemical Theory and Computation associated mixtures formed by methanol and propanal. J. Chem. Soc., Faraday Trans. 1995, 91, 2779−2787. (184) Owen, K.; Quayle, O. R.; Clegg, W. J. A study of organic parachors. V. Constitutive variations of the parachors of a series of normal ketones. J. Am. Chem. Soc. 1942, 64, 1294−1296. (185) González, B.; Domínguez, A.; Tojo, J. Physical properties of the binary systems methylcyclopentane with ketones (acetone, butanone and 2-pentanone) at T = (293.15, 298.15, and 303.15) K. New UNIFAC-VISCO interaction parameters. J. Chem. Thermodyn. 2006, 38, 707−716. (186) Thole, F. B. CLXXXVIthe effect of ring formation on viscosity. J. Chem. Soc., Trans. 1914, 105, 2004−2012. (187) Suárez-Iglesias, O.; Medina, I.; Sanz, M. A.; Pizarro, C.; Bueno, J. L. Self-diffusion in molecular fluids and noble gases: available data. J. Chem. Eng. Data 2015, 60, 2757−2817. (188) Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water. J. Chem. Phys. 1983, 79, 2375−2389. (189) Postma, J. P. M. MD of H2O, a Molecular Dynamics Study of Water. Thesis, University of Groningen, 1985. (190) Waheed, Q.; Edholm, O. Quantum corrections to classical molecular dynamics simulations of water and ice. J. Chem. Theory Comput. 2011, 7, 2903−2909. (191) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force field benchmark of organic liquids: density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory Comput. 2012, 8, 61−74. (192) Kunz, A.-P. E.; van Gunsteren, W. F. Development of a nonlinear classical polarization model for liquid water and aqueous solutions: COS/D. J. Phys. Chem. A 2009, 113, 11570−11579. (193) Zhang, Y.; Feller, S. E.; Brooks, B. R.; Pastor, R. W. Computer simulation of liquid/liquid interfaces. I. Theory and application to octane/water. J. Chem. Phys. 1995, 103, 10252−10266. (194) Einstein, A. Ü ber die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann. Phys. (Berlin, Ger.) 1905, 322, 549−560. (195) Smith, P. E.; van Gunsteren, W. F. The viscosity of SPC and SPC/E water at 277 and 300 K. Chem. Phys. Lett. 1993, 215, 315−318. (196) Neumann, M. Dipole moment fluctuation formulas in computer simulations of polar systems. Mol. Phys. 1983, 50, 841−858. (197) Hagler, A. T.; Lifson, S. Energy functions for peptides and proteins. II. The amide hydrogen bond and calculation of amide crystal properties. J. Am. Chem. Soc. 1974, 96, 5327−5335. (198) Hagler, A. T.; Lifson, S. A procedure for obtaining energy parameters from crystal packing. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1974, 30, 1336−1341. (199) Hagler, A.T.; Lifson, S.; Dauber, P. Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 2. A benchmark for the objective comparison of alternative force fields. J. Am. Chem. Soc. 1979, 101, 5122−5130. (200) Hagler, A. T.; Dauber, P.; Lifson, S. Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 3. The C=O···H hydrogen bonds and the analysis of the energetics and packing of carboxylic acids. J. Am. Chem. Soc. 1979, 101, 5131−5141. (201) Mı ́guez, J. M.; Piñeiro, M. M.; Blas, F. J. Influence of the longrange corrections on the interfacial properties of molecular models using Monte Carlo simulation. J. Chem. Phys. 2013, 138, 034707/1− 034707/10. (202) Zubillaga, R. A.; Labastida, A.; Cruz, B.; Martı ́ez, J. C.; Sánchez, E.; Alejandre, J. Surface tension of organic liquids using the OPLS/AA force field. J. Chem. Theory Comput. 2013, 9, 1611−1615. (203) Fischer, N. M.; van Maaren, P. J.; Ditz, J. C.; Yildirim, A.; van der Spoel, D. Properties of organic liquids when simulated with longrange Lennard-Jones interactions. J. Chem. Theory Comput. 2015, 11, 2938−2944. (204) Adachi, K.; Suga, H.; Seki, S. Phase changes in crystalline and glassy-crystalline cyclohexanol. Bull. Chem. Soc. Jpn. 1968, 41, 1073− 1087.

(205) Suga, H.; Seki, S. Thermodynamic investigation on glassy states of pure simple compounds. J. Non-Cryst. Solids 1974, 16, 171− 194. (206) Pelous, J.; Vacher, R.; Bonnet, J. P.; Ribet, J. L. Brillouin scattering study of the orientationally disordered phase of cyclohexanol. Phys. Lett. A 1980, 78, 195−197. (207) Kuhns, P. L.; Conradi, M. S. NMR study of molecular motions in cyclohexanol, a glass-forming rotor crystal. J. Chem. Phys. 1984, 80, 5851−5858. (208) Suzuki, H.; Hoshina, H.; Otani, C. Kinetics of polymorphic transitions of cyclohexanol investigated by terahertz absorption spectroscopy. Cryst. Growth Des. 2014, 14, 4087−4093. (209) Kang, M.; Smith, P. E. A Kirkwood-Buff derived force field for amides. J. Comput. Chem. 2006, 27, 1477−1485. (210) Whitfield, T. W.; Martyna, G. J.; Allison, S.; Bates, S. P.; Vass, H.; Crain, J. Structure and hydrogen bonding in neat N-methylacetamide: classical molecular dynamics and raman spectroscopy studies of a liquid of peptide fragments. J. Phys. Chem. B 2006, 110, 3624−3637. (211) Zhang, R.; Zheng, D.; Pan, Y.; Luo, S.; Wu, W.; Li, H. All-atom simulation and excess properties study on intermolecular interactions of amide−water system. J. Mol. Struct. 2008, 875, 96−100. (212) Harder, E.; Anisimov, V. M.; Whitfield, T.; MacKerell, A. D.; Roux, B. Understanding the dielectric properties of liquid amides from a polarizable force field. J. Phys. Chem. B 2008, 112, 3509−3521. (213) Lin, B.; Lopes, P. E. M.; Roux, B.; MacKerell, A. D., Jr Kirkwood-Buff analysis of aqueous N-methylacetamide and acetamide solutions modeled by the CHARMM additive and Drude polarizable force fields. J. Chem. Phys. 2013, 139, 084509/1−084509/11. (214) Trabelsi, S.; Bahri, M.; Nasr, S. X-ray scattering and densityfunctional theory calculations to study the presence of hydrogenbonded clusters in liquid N-methylacetamide. J. Chem. Phys. 2005, 122, 024502/1−024502/8. (215) Gainaru, C.; Bauer, S.; Vynokur, E.; Wittkamp, H.; Hiller, W.; Richert, R.; Böhmer, R. Dynamics in supercooled secondary amide mixtures: dielectric and hydrogen bond specific spectroscopies. J. Phys. Chem. B 2015, 119, 15769−15779. (216) Jiang, X.-N.; Wang, C.-S. Rapid prediction of hydrogen bond cooperativity in N-nethylacetamide chains. ChemPhysChem 2009, 10, 3330−3336. (217) Wolf, M. G.; Groenhof, G. Evaluating nonpolarizable nucleic acid force fields: a systematic comparison of the nucleobases hydration free energies and chloroform-to-water partition coefficients. J. Comput. Chem. 2012, 33, 2225−2232. (218) Hagler, A. T. Quantum derivative fitting and biomolecular force fields: functional form, coupling terms, charge flux, nonbond anharmonicity, and individual dihedral potentials. J. Chem. Theory Comput. 2015, 11, 5555−5572. (219) Peña, M. D.; Pando, C.; Renuncio, J. A. R. Combination rules for intermolecular potential parameters. I. Rules based on approximations for the long-range dispersion energy. J. Chem. Phys. 1982, 76, 325−332. (220) Wennberg, C. L.; Murtola, T.; Páll, S.; Abraham, M. J.; Hess, B.; Lindahl, E. Direct-space corrections enable fast and accurate Lorentz-Berthelot combination rule Lennard-Jones lattice summation. J. Chem. Theory Comput. 2015, 11, 5737−5746. (221) Chen, C.; Depa, P.; Sakai, V. G.; Maranas, J. K.; Lynn, J. W.; Peral, I.; Copley, J. R. D. A comparison of united atom, explicit atom, and coarse-grained simulation models for poly(ethylene oxide). J. Chem. Phys. 2006, 124, 234901/1−234901/11. (222) Chen, C.; Depa, P.; Maranas, J. K.; Sakai, V. G. Comparison of explicit atom, united atom, and coarse-grained simulations of poly(methyl methacrylate). J. Chem. Phys. 2008, 128, 124906/1− 124906/12. (223) Perez-Pellitero, J.; Bourasseau, E.; Demachy, I.; Ridard, J.; Ungerer, P.; Mackie, A. D. Anisotropic united-atoms (AUA) potential for alcohols. J. Phys. Chem. B 2008, 112, 9853−9863. 3848


Article

Journal of Chemical Theory and Computation (224) Li, C.; Choi, P.; Sundararajan, P. R. Simulation of chain folding in polyethylene: a comparison of united atom and explicit hydrogen atom models. Polymer 2010, 51, 2803−2808. (225) Liu, X.; Zhang, X.; Zhou, G.; Yao, X.; Zhang, S. All-atom and united-atom simulations of guanidinium-based ionic liquids. Sci. China: Chem. 2012, 55, 1573−1579. (226) Hemmen, A.; Gross, J. Transferable anisotropic united-atom force field based on the Mie potential for phase equilibrium calculations: n-alkanes and n-olefins. J. Phys. Chem. B 2015, 119, 11695−11707. (227) Halgren, T. A. Representation of van der Waals (vdW) interactions in molecular mechanics force-fields: Potential form, combination rules, and vdW parameters. J. Am. Chem. Soc. 1992, 114, 7827−7843. (228) Waldman, M.; Hagler, A. T. New combining rules for rare gas van der Waals parameters. J. Comput. Chem. 1993, 14, 1077−1084. (229) Song, W.; Rossky, P. J.; Maroncelli, M. Modeling alkane− perfluoroalkane interactions using all-atom potentials: failure of the usual combining rules. J. Chem. Phys. 2003, 119, 9145−9162. (230) Al-Matar, A. K.; Rockstraw, D. A. A generating equation for mixing rules and two new mixing rules for interatomic potential energy parameters. J. Comput. Chem. 2004, 25, 660−668. (231) Goodwin, A. R. H.; Sandler, S. I. In Applied Thermodynamics of Fluids; Goodwin, A. R. H., Sengers, J. V., Peters, C. J., Eds.; Royal Society of Chemistry: Cambridge, U.K., 2010; pp 84−134. (232) Fyta, M.; Netz, R. R. Ionic force field optimization based on single-ion and ion-pair solvation properties: going beyond mixing rules. J. Chem. Phys. 2012, 136, 124103/1−124103/11. (233) Kastenholz, M. A.; Hünenberger, P. H. Influence of artificial periodicity and ionic strength in molecular dynamics simulations of charged biomolecules employing lattice-sum methods. J. Phys. Chem. B 2004, 108, 774−788. (234) Reif, M. M.; Kräutler, V.; Kastenholz, M. A.; Daura, X.; Hünenberger, P. H. Explicit-solvent molecular dynamics simulations of a reversibly-folding β-heptapeptide in methanol: influence of the treatment of long-range electrostatic interactions. J. Phys. Chem. B 2009, 113, 3112−3128. (235) Halgren, T. A.; Damm, W. Polarizable force fields. Curr. Opin. Struct. Biol. 2001, 11, 236−242. (236) Yu, H.; van Gunsteren, W. F. Accounting for polarization in molecular simulation. Comput. Phys. Commun. 2005, 172, 69−85. (237) Vorobyov, I. V.; Anisimov, V. M.; MacKerell, A. D., Jr Polarizable force field for alkanes based on the classical Drude oscillator model. J. Phys. Chem. B 2005, 109, 18988−18999. (238) Warshel, A.; Kato, M.; Pisliakov, A. V. Polarizable force fields: history, test cases, and prospects. J. Chem. Theory Comput. 2007, 3, 2034−2045. (239) Lopes, P. E.; Lamoureux, G.; Roux, B.; MacKerell, A. D., Jr Polarizable empirical force field for aromatic compounds based on the classical Drude oscillator. J. Phys. Chem. B 2007, 111, 2873−2885. (240) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; HeadGordon, T. Current status of the AMOEBA polarizable force field. J. Phys. Chem. B 2010, 114, 2549−2564. (241) Ponomarev, S. Y.; Kaminski, G. A. Polarizable simulations with second-order interaction model (POSSIM) force field: developing parameters for alanine peptides and protein backbone. J. Chem. Theory Comput. 2011, 7, 1415−1427. (242) Ren, P.; Wu, C.; Ponder, J. W. Polarizable atomic multipolebased molecular mechanics for organic molecules. J. Chem. Theory Comput. 2011, 7, 3143−3161. (243) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D., Jr Polarizable force field for peptides and proteins based on the classical Drude oscillator. J. Chem. Theory Comput. 2013, 9, 5430−5449. (244) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable atomic multipole-based AMOEBA force field for proteins. J. Chem. Theory Comput. 2013, 9, 4046−4063.

(245) Bachmann, S. J.; van Gunsteren, W. F. Polarizable model for DMSO and DMSO−water mixtures. J. Phys. Chem. B 2014, 118, 10175−10186. (246) Szklarczyk, O. M.; Bachmann, S. J.; van Gunsteren, W. F. A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons. J. Comput. Chem. 2014, 35, 789−801. (247) Baker, C. M. Polarizable force fields for molecular dynamics simulations of biomolecules. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2015, 5, 241−254. (248) Jakobsen, S.; Jensen, F. Systematic improvement of potentialderived atomic multipoles and redundancy of the electrostatic parameter space. J. Chem. Theory Comput. 2014, 10, 5493−5504. (249) Cardamone, S.; Hughes, T. J.; Popelier, P. L. A. Multipolar electrostatics. Phys. Chem. Chem. Phys. 2014, 16, 10367−10387. (250) Cardamone, S.; Popelier, P. L. A. Prediction of conformationally dependent atomic multipole moments in carbohydrates. J. Comput. Chem. 2015, 36, 2361−2373. (251) Palmo, K.; Mannfors, B.; Krimm, S. Balanced charge treatment of intramolecular electrostatic interactions in molecular mechanics energy functions. Chem. Phys. Lett. 2003, 369, 367−373. (252) Palmo, K.; Mannfors, B.; Mirkin, N. G.; Krimm, S. Inclusion of charge and polarizability fluxes provides needed physical accuracy in molecular mechanics force fields. Chem. Phys. Lett. 2006, 429, 628−632. (253) Stachowicz, A.; Styrcz, A.; Korchowiec, J. Charge sensitivity analysis in force-field-atom resolution. J. Mol. Model. 2011, 17, 2217− 2226. (254) Stachowicz, A.; Korchowiec, J. Generalized charge sensitivity analysis. Struct. Chem. 2012, 23, 1449−1458. (255) Ahmed, A.; Sandler, S. I. Hydration free energies of multifunctional nitroaromatic compounds. J. Chem. Theory Comput. 2013, 9, 2774−2785. (256) Jämbeck, J. P. M.; Mocci, F.; Lyubartsev, A. P.; Laaksonen, A. Partial atomic charges and their impact on the free energy of solvation. J. Comput. Chem. 2013, 34, 187−197. (257) Rai, B. K.; Bakken, G. A. Fast and accurate generation of ab initio quality atomic charges using nonparametric statistical regression. J. Comput. Chem. 2013, 34, 1661−1671. (258) Shivakumar, S.; Deng, Y.; Roux, B. Computations of absolute solvation free energies of small molecules using explicit and implicit solvent model. J. Chem. Theory Comput. 2009, 5, 919−930. (259) Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the OPLS force field. J. Chem. Theory Comput. 2010, 6, 1509−1519. (260) Liu, Y.; Fu, J.; Wu, J. High-throughput prediction of the hydration free energies of small molecules from a classical density functional theory. J. Phys. Chem. Lett. 2013, 4, 3687−3691. (261) Knight, J. L.; Yesselman, J. D.; Brooks, C. L., III Assessing the quality of absolute hydration free energies among CHARMMcompatible ligand parameterization schemes. J. Comput. Chem. 2013, 34, 893−903. (262) Martins, S. A.; Sousa, S. F.; Ramos, M. J.; Fernandes, P. A. Prediction of solvation free energies with thermodynamic integration using the general AMBER force field. J. Chem. Theory Comput. 2014, 10, 3570−3577. (263) Misin, M.; Fedorov, M. V.; Palmer, D. S. Hydration free energies of molecular ions from theory and simulation. J. Phys. Chem. B 2016, 120, 975−983. (264) Nerenberg, P. S.; Jo, B.; So, C.; Tripathy, A.; Head-Gordon, T. Optimizing solute−water van der Waals interactions to reproduce solvation free energies. J. Phys. Chem. B 2012, 116, 4524−4534. (265) Zhang, J.; Tuguldur, B.; van der Spoel, D. Force field benchmark of organic liquids. 2. Gibbs energy of solvation. J. Chem. Inf. Model. 2015, 55, 1192−1201. (266) Beauchamp, K. A.; Behr, J. M.; Rustenburg, A. S.; Bayly, C. I.; Kroenlein, K.; Chodera, J. D. Toward automated benchmarking of atomistic force fields: neat liquid densities and static dielectric constants from the ThermoML data archive. J. Phys. Chem. B 2015, 119, 12912−12920. 3849


Article

Journal of Chemical Theory and Computation (267) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on the generalized Born approximation with asymmetric descreening. J. Chem. Theory Comput. 2009, 5, 2447− 2464. (268) Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J. Comput.-Aided Mol. Des. 2014, 28, 711−720. (269) Fu, J.; Wu, J. Toward high-throughput predictions of the hydration free energies of small organic molecules from first principles. Fluid Phase Equilib. 2016, 407, 304−313. (270) Krieger, E.; Koraimann, G.; Vriend, G. Increasing the precision of comparative models with YASARA NOVA. A self-parameterizing force field. Proteins: Struct., Funct., Genet. 2002, 47, 393−402. (271) Schüttelkopf, A. W.; van Aalten, D. M. F. PRODRGa tool for high-throughput crystallography of protein-ligand complexes. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2004, 60, 1355−1363. (272) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (273) Ribeiro, A. A. S. T.; Horta, B. A. C.; de Alencastro, R. B. MKTOP: a program for automatic construction of molecular topologies. J. Braz. Chem. Soc. 2008, 19, 1433−1435. (274) Zoete, V.; Cuendet, M. A.; Grosdidier, A.; Michielin, O. SwissParam: a fast force field generation tool for small organic molecules. J. Comput. Chem. 2011, 32, 2359−2368. (275) Vanommeslaeghe, K.; MacKerell, A. D., Jr Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. J. Chem. Inf. Model. 2012, 52, 3144−3154. (276) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D., Jr Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J. Chem. Inf. Model. 2010, 52, 3155−3168. (277) Malde, A. K.; Zuo, L.; Breeze, M.; Stroet, M.; Poger, D.; Nair, P. C.; Oostenbrink, C.; Mark, A. E. An automated force field topology builder (ATB) and repository: version 1.0. J. Chem. Theory Comput. 2011, 7, 4026−4037. (278) Zhang, Q.; Zhang, W.; Li, Y.; Wang, J.; Zhang, L.; Hou, T. A rule-based algorithm for automatic bond type perception. J. Cheminf. 2012, 4, 26. (279) Yesselman, J. D.; Price, D. J.; Knight, J. L.; Brooks, C. L., III MATCH: an atom-typing toolset for molecular mechanics force fields. J. Comput. Chem. 2012, 33, 189−202. (280) Koziara, K. B.; Stroet, M.; Malde, A. K.; Mark, A. E. Testing and validation of the automated topology builder (ATB) version 2.0: prediction of hydration free enthalpies. J. Comput.-Aided Mol. Des. 2014, 28, 221−233. (281) Gapsys, V.; Michielssens, S.; Seeliger, D.; de Groot, B. L. pmx: automated protein structure and topology generation for alchemical perturbations. J. Comput. Chem. 2015, 36, 348−354. (282) Loeffer, H. H.; Michel, J.; Woods, C. FESetup: automating setup for alchemical free energy simulations. J. Chem. Inf. Model. 2015, 55, 2485−2490. (283) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a full periodic-table force-field for molecular mechanics and molecular-dynamics simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (284) Casewit, C. J.; Colwell, K. S.; Rappé, A. K. Application of a universal force field to organic molecules. J. Am. Chem. Soc. 1992, 114, 10035−10046. (285) Halgren, T. A. Merck molecular force field. 1. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 1996, 17, 490−519. (286) Halgren, T. A. Merck molecular force field. 2. MMFF94 van der waals and electrostatic parameters for intermolecular interactions. J. Comput. Chem. 1996, 17, 520−552. (287) Halgren, T. A. MMFF VII. Characterization of MMFF94, MMFF94s, and other widely available force fields for conformational energies and for intermolecular-interaction energies and geometries. J. Comput. Chem. 1999, 20, 730−748. 3850


A GROMOS-Compatible Force Field for Small Organic Molecules in

Recommend Documents