Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant.

The chemical composition of small organic molecules is often very similar to amino acid side chains or the bases in nucleic acids, and hence there is no a priori reason why a molecular mechanics force field could not describe both organic liquids and biomolecules with a single parameter set. Here, we devise a benchmark for force fields in order to test the ability of existing force fields to reproduce some key properties of organic liquids, namely, the density, enthalpy of vaporization, the surface tension, the heat capacity at constant volume and pressure, the isothermal compressibility, the volumetric expansion coefficient, and the static dielectric constant. Well over 1200 experimental measurements were used for comparison to the simulations of 146 organic liquids. Novel polynomial interpolations of the dielectric constant (32 molecules), heat capacity at constant pressure (three molecules), and the isothermal compressibility (53 molecules) as a function of the temperature have been made, based on experimental data, in order to be able to compare simulation results to them. To compute the heat capacities, we applied the two phase thermodynamics method (Lin et al. J. Chem. Phys.2003, 119, 11792), which allows one to compute thermodynamic properties on the basis of the density of states as derived from the velocity autocorrelation function. The method is implemented in a new utility within the GROMACS molecular simulation package, named g_dos, and a detailed exposé of the underlying equations is presented. The purpose of this work is to establish the state of the art of two popular force fields, OPLS/AA (all-atom optimized potential for liquid simulation) and GAFF (generalized Amber force field), to find common bottlenecks, i.e., particularly difficult molecules, and to serve as a reference point for future force field development. To make for a fair playing field, all molecules were evaluated with the same parameter settings, such as thermostats and barostats, treatment of electrostatic interactions, and system size (1000 molecules). The densities and enthalpy of vaporization from an independent data set based on simulations using the CHARMM General Force Field (CGenFF) presented by Vanommeslaeghe et al. (J. Comput. Chem.2010, 31, 671) are included for comparison. We find that, overall, the OPLS/AA force field performs somewhat better than GAFF, but there are significant issues with reproduction of the surface tension and dielectric constants for both force fields.

Introduction

Parameters in most force fields have been derived incrementally, that is, building on previous work by adding support for different chemical moieties in a sequential fashion. While the focus of many force fields is on biomolecules, the chemical basis lies in organic molecules. Of the major force fields available today OPLS/AA (optimized parameters for liquid simulations, all atoms) is one of the few that “specializes” in simple liquids.[1] The generalized Amber force field (GAFF) was introduced recently[2] (together with the Antechamber set of programs[3]) to aid in the derivation of force field parameters for small molecules that are often involved in binding to biomolecules. Accurate parameters are crucial for predicting, for instance, the Gibbs energy of ligand binding, a key property in drug design.[4] The GAFF parameters for small molecules are intended to be combined with the Amber force field[5] although there are studies of proteins using GAFF parameters as well.[6]

A critical component in force field development is generation of partial charges. The method for deriving partial charges by optimizing their values to reproduce the electrostatic potential (ESP) was introduced in the 1980s by Kollman et al.[7,8] The electron density taken from a quantum chemistry calculation, together with the nuclear charges, generates an electrostatic potential around the molecule. Typically, the set of partial charges for a molecule, for use in force field calculations, is determined by minimizing the (square) difference between the ESP generated by the partial charges and the ESP generated by the quantum chemistry calculation. A set of partial charges (or indeed any atom-centered set of spherically distributed charges) can never completely reproduce the ESP due to the fact that electron density is not completely spherically symmetric around the nuclei (for instance, due to p and higher orbitals). A further issue is due to the fact that the fitting points are highly correlated, and hence atoms far from the ESP data points (e.g., the buried carbon in isobutanol) may end up being a sink for the fit[9,10] and get arbitrary values. An ad hoc refinement of the ESP method to overcome this problem is the restrained ESP (RESP) method.[11] The RESP method does the same fit, however with an added penalty on the absolute value of the charge. The RESP method is an integral part of the Antechamber package,[2,3] which relies on either quantum calculations or empirical methods, such as AM1-BCC,[12,13] to provide the partial charges.

Mobley et al. tested the performance of GAFF parameters for Gibbs energies of hydration using two different water models.[14,15] They paid particular attention to the way the partial charges were determined and found that the final results are related to the level of theory used, something that was corroborated by Wallin et al., who did a similar study of charge schemes for ligand binding.[16] The CM1 charge model for OPLS/AA,[17] used in the study of Wallin et al.,[16] performs well,[18,19] although some degradation for conformational energetics is expected. The differences are generally considered to be minor.[1] There are some drawbacks with these studies however. First, they involve complex systems, where a subset of the parameters was changed and the “quality” of the charges evaluated on the basis of a single number, the free energy, hereby ignoring the interdependency between Lennard-Jones parameters and point charges. Second, free energy calculations depend critically on the amount of the sampling that was used, although it is possible to ascertain that the errors due to sampling are small.[20] In order to test the validity of force field, it would be good to take one step back and evaluate the performance for simple systems first, in order to avoid systematic errors due to water model and/or protein force fields. A recent review by Jorgensen and Tirado-Rives provides further background information on the topic of force field development.[1]

To assess the state of the art of GAFF and OPLS/AA force fields, we provide a comprehensive benchmark of the liquid properties of molecules in each of the GAFF and OPLS/AA force fields. Previous simulations of mixtures of alcohol and water[21,22] using the OPLS/AA force field showed that many properties of the pure liquids are reproduced faithfully, but the heat of mixing and the density of mixing are slightly, but significantly, off. Similar comparisons of force fields for water models are numerous in the literature (see for example, refs (23−29)), while for organic liquids there are some papers by Kaminski and Jorgensen,[30,31] and a recent paper by Wang and Tingjun,[32] which we discuss in the Discussion section.

Liquid properties are usually known experimentally with high accuracy, and their calculation is most often straightforward. Rather, the time goes into the preparation and equilibration of the systems. A total of 146 molecular liquids was prepared and simulated using these force fields in the GROMACS molecular simulation package,[33−36] and from these molecular dynamics simulations, we extract the density ρ (from constant pressure simulations), the enthalpy of vaporization ΔHvap, the heat capacities at constant pressure cP and volume cV, the volumetric expansion coefficient αP, the isothermal compressibility κT, the surface tension γ, and the static dielectric constant ε(0). Although, in principle, more observables could be computed, this set includes the most important thermodynamic properties of the liquids, including temperature derivatives of energy and volume. The intention of this work is to supply a large number of tests for further force field development. To this end, the topologies and structures have been made available on a dedicated Web site at http://virtualchemistry.org, while the simulation parameters are available as Supporting Information to this paper. These topologies and structure files may be useful for simulations of biomolecules in organic liquids as well. The recently presented all atom CHARMM general force field (CGenFF)[37] would be an equally well suited candidate for inclusion in this comparison, but we have chosen to limit our simulations to two force fields only. However, to allow the reader to compare OPLS/AA and GAFF to a similar study based on CGenFF, we have included results on density and enthalpy of vaporization from that paper.[37]

Methods

Energy Function

Most force fields use the same functional form for the intermolecular part of the interaction function, based on the Coulomb potential and the Lennard-Jones potential:where r is the distance between two atoms i and j, q and q are the partial charges on the atoms, ε0 is the permittivity of vacuum, σ is the van der Waals radius, and ε is the well-depth for this atom pair. In most force fields, the parameters σ and ε are derived from the atomic values σ and ε using a simple equation (the combination rule). Suffice to say that we have applied the standard combination rules for GAFF (Lorentz–Berthelot[38]) and for OPLS/AA (σ = (σσ)1/2 and ε = (εε)1/2[39]) in this work.

Molecule Selection and Preparation

A set of organic molecules was selected for which both enthalpy of vaporization and density are known at room temperature. Models for these molecules were built using either PRODRG[40] or Molden.[41] These molecules were optimized using the Gaussian 03 suite of programs[42] at the Hartree–Fock level with the 6-311G** basis set.[43−47]

OPLS/AA Topologies

The OpenBabel (http://openbabel.org) code was used to extract a coordinate file including connectivity information from the Gaussian output files, and this file was used to generate an initial topology using the GROMACS tools[35] for the OPLS/AA force field.[1,39] The topologies were checked manually for correctness before using them, making sure that the total charge of the molecule is zero, and also that the atom types were correct. For molecules containing linear groups (e.g., nitriles), a virtual site construction was added to the topologies preserving the moment of inertia and the total mass, in order to keep the groups perfectly linear.[48]

GAFF Topologies

For the simulations where GAFF[2] was used, the Antechamber software[2,3] was employed to generate the topologies from the coordinate files (which were generated as explained above). Gaussian 03[42] at the Hartee–Fock level with the 6-311G** basis set[43−47] (as provided by the Basis Set Exchange Web site[49,50]) and Merz–Singh–Kollman (MK) scheme[7] were used the to determine the partial charges in Gaussian. This particular basis set was used because it is very similar to the 6/31G* basis set,[51] which is the default for GAFF, while simultaneously supporting a larger number of elements (e.g., I). The MK radius for I is not implemented in Antechamber, we used RI = 2.15 Å. The amb2gmx.pl script[52] was used to convert the AMBER topologies into the GROMACS format (this script is available online at http://ffamber.cnsm.csulb.edu/). The final partial charges were calculated using the RESP method[11] as implemented in Antechamber, and we manually checked that the charges were sane. Note that RESP can be used with any QM method producing electrostatics, not just with HF/6-311G**. No modifications for linear group were made for the GAFF topologies, where the Antechamber software[3] generates a near-linear angle term instead.

Liquid Simulation Box Preparation

To generate liquid simulation boxes, we first made a 2 × 2 × 2 nm3 box containing a single molecule. From 125 such single molecule boxes, we built up a 10 × 10 × 10 nm3 box. These boxes were simulated under high pressure (100 bar) to force the molecules into the liquid phase, and finally we let the systems relax under normal pressure (1 bar) to reach an equilibrated system. For the equilibration simulations, we used Berendsen’s coupling algorithm[53] because of its efficient relaxation properties.[34] To generate our final simulation boxes, we stacked 2 × 2 × 2 of the 125 molecule boxes and ran an additional equilibration simulation. The absolute drift in total energy was automatically checked in the equilibration and production simulations, and the simulations were continued until the drift was below 0.5 J/mol/ns per degree of freedom, which is a very strict criterion but which is necessary to accurately compute fluctuation properties.

Simulation Parameters

The GROMACS suite of programs was used for all simulations.[33−36] Following previous simulations of alcohol water mixtures[21,22] using the OPLS/AA force field,[1,39] we employed a 1.1 nm cutoff for Lennard-Jones interactions and the same distance as the switching distance for the particle mesh Ewald (PME) algorithm for computing Coulomb interactions.[54,55] Although the OPLS/AA force field was not developed for use with PME, extensive studies on water models[56] and proteins in water[57] have shown that correspondence of simulation results with experimental data improves considerably when long-range interactions are taken into account explicitly—irrespective of the force field used. Analytic corrections to pressure and potential energies were made to compensate for the truncation of the Lennard-Jones interactions.[38] In the production simulations, we used the Nosé–Hoover algorithm for temperature coupling,[58,59] in order to provide correct fluctuations, which is necessary to compute fluctuation properties. A time constant for coupling of 1 ps (corresponding to a mass parameter Q of 7.6 ps at room temperature) was used, which is in the range of time scales for intermolecular collisions, as recommended by Holian et al.[60] For production simulations at constant pressure, the Parrinello–Rahman pressure coupling[61] algorithm was used with compressibility set to 5 × 10–5 bar–1 and a time constant of 5 ps. The temperatures of the simulations were selected to fit the experimental data available. In most simulations, the bonds were constrained using the LINCS algorithm[62,63] for all molecules, applying two iterations in order to obtain good energy conservation. Periodic boundary conditions were used in all liquid phase simulations.

Four types of production run simulations were performed according to Table 1. The density of states (DOS) production simulations were performed under constant volume conditions, but they were preceded by equilibration simulations under NPT (without constraints) in order to obtain the equilibrium density at P = 1 bar for the subsequent DOS simulations. In the DOS simulations, slightly stricter energy conservation parameters were used: a neighbor list buffer of 0.3 nm, combined with a switched Lennard-Jones and short-range electrostatics term (1.0–1.1 nm), see reference (56) for a description of the functional form.

Table 1

Simulation Characteristics for the Different Simulation Types

name	length	# molecules	ensemble	constraints	electrostatics
LIQ	10 ns	1000	NPT	all bonds	PME
GAS	100 ns	1	NVT	all bonds	all interactions
SURF	10 ns	1000	NVT	all bonds	PME
DOS	100 ps	1000	NVT	none	PME

The GAS simulations were done using a stochastic dynamics (SD) integrator, which adds a friction and a noise term to Newton’s equation of motion:where m is the mass of atom i, ξ is a friction constant, and ρ(t) is a noise process withwhere kB is Boltzmann’s constant, T is the temperature, δ(s) is the Dirac δ function, and δ is the Kronecker δ function. A leapfrog algorithm adapted for SD simulations[64] was used to integrate eq 2. When 1/ξ is large compared to the time scales present in the system, SD functions like molecular dynamics with stochastic temperature-coupling. One of the benefits with SD as compared to MD is that when simulating a system in a vacuum there is no accumulation of errors for the overall translational and rotational degrees of freedom, making sampling of different configuration states more accurate. SURF and LIQ simulations were done using a conventional MD leapfrog integrator.[65] To enable replication of our simulations and detailed scrutiny of the data, we provide all simulation parameters for each type of run, as well as starting structures and topologies. These files, in GROMACS format, are available for downloading at http://virtualchemistry.org.

To ensure that our liquid systems did not freeze during the simulations, we monitored the changes in diffusion constant ΔD as derived from the mean square displacement during the simulations, defined asThe subscript “begin” means the value is an average over the 1000–1500 ps of the simulation, and “end” means over 8500–9000 ps. |ΔD| is close to zero for most simulations, indicating that D is approximately the same in the beginning and at the end of the simulation. We also verified that D > 0 for all simulations. For the simulations where |ΔD| ≥ 0.5, we ensured that the systems indeed were not frozen, by inspecting the full mean square displacement curve and the trajectory of the simulations. In the Supporting Information (Figure S1), we show ΔD for all of the liquid simulations.

Analysis

The density ρ in a constant pressure simulation follows trivially from the mass M of the system divided by the volume V:The enthalpy of vaporization can be computed fromwhere Eintra is the intramolecular energy in either the gas (g) phase or the liquid (l) phase and Einter represents the intermolecular energy of the system. In practice, we can simply evaluateρ was determined from LIQ simulations and ΔHvap from LIQ and GAS simulations.

The SURF simulations were done using liquid boxes, the size of which in the z direction was extended by a factor of 3, generating a simulation box with two liquid–vacuum interfaces. The surface tension γ then follows fromwhere P is the pressure component in direction n and L is the length of the box in the z direction (perpendicular to the surfaces).

Static dielectric constants ε(0) were computed on the basis of the fluctuations of the total dipole moment M of the simulation box[66,67] in the LIQ simulations:where V is the volume of the simulation box. Errors were estimated by block-averaging over 10 blocks of 1 ns. In order to verify the validity of eq 9, we computed the autocorrelation time τM of the total dipole moment M in the simulation boxes (from the integral of the autocorrelation function). In order for fluctuations to be well-defined, τM should be at least an order of magnitude shorter than the simulation length. Henceforth, we omitted the dielectric constants for those systems where τM was longer than 1 ns. For those systems where this was the case, longer simulations of 50 ns were performed, in most cases without any improvement.

The fluctuation properties αP (the volumetric thermal expansion coefficient) and κT (the isothermal compressibility) are computed from the LIQ simulations according to[38]where H is the enthalpy and δ indicates the fluctuations, andThese two properties can be related to the difference between heat capacities at constant pressure and constant volume throughwhere V is the molecular volume. We can take advantage of this relation in two ways, first by computing αP and κT from our simulations and then computing the constant pressure heat capacity based on the constant volume heat capacity. By using experimental data for αP and κT, we can also establish “experimental” constant volume heat capacities, which are difficult to measure directly. In this work, we have done both, as detailed in the Results and Discussion sections.

The classical—that is, without any quantum corrections—heat capacity cPclass can be obtained from the fluctuations in the enthalpy:[38]Although this is straightforward to calculate, the numbers obtained in this manner are a factor of 2 too high (Table 2). Therefore, we have determined the heat capacities cP and cV on the basis of the two phase thermodynamic method[68−70] (described in the Supporting Information), which is based on the convolution of the density of states with a weighting function based on quantum harmonic oscillators, as introduced originally by Berens et al.[71] The final expression yielding the heat capacity cV isDoSgas and DoSsolid denote the density of states in a gas and a solid, Wgas(ν) and Wsolid(ν) are weighting factors for the same, and cP can be obtained by combining eq 12 and eq 14. For all details and a complete derivation, we refer the reader to the Supporting Information.

Table 2

Statistics of a Linear Fit of Calculated to Experimental Values According to y = ax + ba

force field	N	a	b	RMSD	% dev.	R2
ρ (g/l)
GAFF	235	0.96	58.5	82.9	4	97%
OPLS/AA	235	0.98	20.9	40.4	2	99%
CGenFF[37]	111	1.03	–36.0	26.0	2	99%
OPLS/AA[70]	9	1.01	–24.0	45.3	4	96%
ΔHvap(kJ/mol)
GAFF	231	1.07	0.8	10.6	17	83%
OPLS/AA	231	0.96	3.4	6.5	11	89%
CGenFF[37]	95	0.94	2.4	4.7	7	84%
γ (10–3 N/m)
GAFF	155	0.75	0.9	8.6	23	70%
OPLS/AA	155	0.97	–5.5	7.3	22	89%
ε(0)
GAFF	163	0.27	0.4	15.7	35	55%
OPLS/AA	176	0.16	0.7	15.9	43	55%
αP (10–3/K)
GAFF	221	0.90	0.3	0.3	24	67%
OPLS/AA	221	0.91	0.3	0.3	21	75%
OPLS/AA[70]	9	0.53	0.8	0.7	42	39%
κT (1/GPa)
GAFF	103	0.66	0.0	0.3	27	74%
OPLS/AA	103	0.76	0.1	0.3	19	85%
OPLS/AA[70]	8	0.93	0.0	1.1	59	84%
cP(J/mol K)
GAFF	130	1.08	–30.9	19.8	10	98%
OPLS/AA	132	1.10	–30.2	18.2	10	97%
OPLS/AA[70]	9	0.94	3.5	10.4	7	94%
cV(J/mol K)
GAFF	72	1.02	–17.6	18.8	10	97%
OPLS/AA	72	1.04	–17.9	18.3	9	95%
OPLS/AA[70]	8	1.01	–5.4	10.8	7	95%
cPclass(J/mol K)
GAFF	214	1.77	–21.6	148.3	77	87%
OPLS/AA	214	1.98	–52.8	147.0	73	93%

Uncertainties in the simulation results are used as weights in the fit. The number of (experimental) data points N is given for each property. Root mean square deviation (RMSD) from experimental values, average relative deviation in percent, and the correlation coefficient R2 are given. OPLS/AA results from ref (70) and CGenFF results from ref (37) (using the so called CHARMM generalized force field) are also listed for comparison.

The properties investigated fall into two categories: those that follow directly from the ensemble average of a property (energy, pressure, volume) and those based on fluctuations (heat capacities, compressibility, and expansion coefficient). For the first category, error estimates were based on a block averaging procedure that automatically takes the autocorrelation of the property under investigation into account.[72] Properties like potential energy and density usually have relatively short autocorrelation times. The surface tension fluctuates significantly but also has a short autocorrelation time. For the second category, we have used a different approach when estimating the error. By dividing the entire simulation trajectory into nine, in time, equally long parts, we get nine values for each property, from which we can estimate the total error. In the case of cV, we used five blocks of 20 ps for error estimation instead.

We calculated cP on the basis of eq 12 and estimated the error δcP from the errors in cV (δcV), αP (δαP), and κT (δκT) asor, expressed in Δc (eq 12):

Results

Correlations between experimental data and simulations for observables and derived quantities are plotted in Figures 1–8. The statistics for linear fits to the data (ycalcd = ayexptl + b) are given in Table 2 for each of the observables and the two force fields, plus similar data from refs (37) and (70). To identify which specific molecule generated a certain value in the figures, we refer to Tables S2–S10 in the Supporting Information. An overview of the names of the molecules, their formula, molecular weight, CAS number, and ChemSpider ID is given in Table S1 (Supporting Information). For many molecules, results at different temperatures were generated, and hence the number of data points may be larger than the number of molecules. For densities, heats of vaporization, surface tensions, and dielectric constants, some of the experimental values were generated from analytical functions of temperature based on experimental data, the parameters of which are given in the Handbook of Chemistry and Physics,[73] the Landolt-Bornstein database,[74] and Yaws’ book on Thermophysical Properties of Chemicals and Hydrocarbons.[75] In addition, we parametrized the dielectric constant, heat capacity at constant pressure, and isothermal compressibility as a function of the temperature for some molecules (see below).

Figure 1

Correlation between densities (ρ) calculated by MD simulation using GAFF, OPLS/AA, CGenFF, and experimental results. The CGenFF data were adopted from Vanommeslaeghe et al.[37] and are based on a different (but similar) set of molecules, including 111 molecules. For a full list of the CGenFF data, we refer to the reference and the supplemental files therein.

Figure 8

Correlation between measured heat capacity at constant volume (cV) and computed using the GAFF and the OPLS/AA force fields based on the density of states method, which includes quantum corrections.

Statistics

In the following, we discuss general trends in all properties first; outliers are described separately below. A comparison of the values in Table 2 shows that OPLS/AA is slightly better than GAFF at reproducing experimental data for most observables, with both lower RMSD and higher correlation coefficients R2.

Density

The density ρ (Figure 1, Table S2) of virtually all liquids is reproduced very well, with R2 = 97% (GAFF) and 99% (OPLS/AA) (Table 2). For GAFF, the densities are systematically slightly underestimated (a = 0.96), while for OPLS/AA, a = 0.98, very close to 1, and both have an R2 close to 100%. In a recent publication, Vanommeslaeghe et al. presented the CHARMM general force field (CGenFF).[37] They calculated densities for a set of 111 drug-like molecules, using boxes of 216 molecules. Their reported densities are also very accurate with a = 1.03 and R2 = 99%, see Figure 1 and Table S2.

Enthalpy of Vaporization

ΔHvap (Figure 2, Table S4) correlates very well with experimental data in most cases, with R2 = 83% (GAFF) and 89% (OPLS/AA) (Table 2). The GAFF overestimates ΔHvap with slope a = 1.07, while OPLS/AA underestimates a slightly at 0.96. These deviations cannot just be attributed to a small number of outliers, as may be evident from Figure 2. Vanommeslaeghe et al.[37] calculated enthalpy of vaporization for a set of 95 small molecules. Like for OPLS/AA, ΔHvap is underestimated in CGenFF calculations with a slope of a = 0.94. The correlation between experiments and simulation is similar to the two force fields studied here, R2 = 84%. The CGenFF data set is based on a comparable but different set of molecules than what has been analyzed here (37 molecules overlap between the two studies). To simplify a comparison between OPLS/AA, GAFF, and CGenFF, we have listed the CGenFF ΔHvap values from the study by Vanommeslaeghe et al. next to OPLS/AA and GAFF values in Table S3, and we have plotted them in Figure 2.

Figure 2

Correlation between enthalpy of vaporization (ΔHvap) calculated using GAFF, OPLS/AA, CGenFF, and experimental results. The CGenFF data were adopted from Vanommeslaeghe et al.[37] and are based on a different (but similar) set of molecules, including 95 molecules. For a full list of the CGenFF data, we refer to the reference and the supplemental files therein.

Surface Tension

The surface tension γ (Figure 3, Table S4) seems to be underestimated systematically in both force fields with slope a = 0.75 (GAFF) and 0.97 (OPLS/AA, Table 2). The interactions between molecules on the surface are not sufficiently strong, a well know problem with nonpolarizable force fields.[21,25,76] The values are spread around the diagonal for both GAFF (R2 = 70%) and OPLS/AA (R2 = 89%), and here again OPLS/AA performs slightly better than GAFF.

Figure 3

Correlation between surface tension (γ) calculated using the GAFF and the OPLS/AA force fields and experimental results.

Dielectric Constant

For 32 molecules, a novel parametrization of the temperature dependence of the dielectric constant was made on the basis of experimental values predominantly from the Landolt-Bornstein database.[77] The parametrization is to a polynomial of second or third order (as is used in the Handbook of Chemistry and Physics[73]), and the resulting coefficients are given in Table 3. Interpolations of these polynomials were used in order to compare the simulations to experimental data, and the fits are presented in Figure S3 of the Supporting Information.

Table 3

Parameterization of Temperature Dependence of Dielectric Constants in a Polynomial Form ε(0) = A + BT + CT2 + DT3, Which Is the Same Form Used in the Handbook of Chemistry and Physics[73]a

molecule	N	χ2	Tmin	Tmax	A	B	C	D
bromomethane	12	0.7	194.60	275.70	52.59	–2.812e–01	4.565e–04	0
methanol	92	1.0	175.62	337.75	226.69	–1.319e+00	2.937e–03	–2.359e–06
1,1,1,2,2-pentachloroethane	9	0.0	245.15	338.15	13.81	–5.527e–02	7.186e–05	0
1,1,2,2-tetrachloroethane	14	0.2	231.15	318.15	71.61	–3.630e–01	5.010e–04	0
1,2-dibromoethane	39	0.1	288.15	353.15	10.31	–3.114e–02	4.200e–05	0
1,1-dichloroethane	8	0.2	288.15	323.15	36.77	–1.300e–01	1.361e–04	0
2-chloroethanol	30	3.1	263.15	401.75	105.36	–3.245e–01	3.619e–05	5.019e–07
ethanamide	7	0.3	358.15	448.20	–200.55	1.551e+00	–2.239e–03	0
methyldisulfanylmethane	6	0.0	293.15	323.15	53.55	–2.539e–01	3.571e–04	0
2-aminoethanol	7	0.4	283.65	298.15	166.68	–7.576e–01	1.018e–03	0
1,3-dioxolan-2-one	24	0.5	309.46	364.15	223.34	–4.560e–01	9.143e–05	0
1,3-dioxolane	31	0.2	175.93	303.15	40.61	–2.507e–01	6.323e–04	–5.695e–07
dimethoxymethane	5	0.0	170.65	298.15	2.59	–9.298e–04	3.847e–06	0
ethylsulfanylethane	6	0.1	293.15	323.15	11.68	–1.994e–02	0.000e+00	0
2-methylpropan-2-amine	4	0.0	291.15	303.15	294.70	–1.887e+00	3.060e–03	0
thiophene	14	0.1	252.65	333.15	2.32	5.071e–03	–1.232e–05	0
furan	31	0.2	198.15	303.15	6.69	–2.044e–02	2.644e–05	0
pentane-2,4-dione	9	2.0	291.15	323.15	–532.57	3.658e+00	–5.982e–03	0
3-methylpyridine	6	1.0	293.15	333.00	35.54	–9.303e–02	4.307e–05	0
benzenethiol	6	0.3	293.15	358.15	5.72	–7.033e–03	7.362e–06	0
(E)-hex-2-ene	6	0.0	157.00	295.00	2.43	–1.132e–03	–1.372e–06	0
1-methoxy-2-(2-methoxyethoxy)ethane	5	0.0	298.15	333.15	32.07	–1.359e–01	1.766e–04	0
diethyl propanedioate	7	0.2	293.15	343.15	19.98	–5.034e–02	3.345e–05	0
2,4,6-trimethylpyridine	10	0.1	293.15	358.15	16.67	–3.036e–02	2.361e–06	0
triethyl phosphate	6	0.1	294.15	333.15	–1.59	1.317e–01	–2.780e–04	0
phenylmethanol	26	0.2	278.15	363.15	105.48	–5.130e–01	6.802e–04	0
tetrahydrothiophene 1,1-dioxide	57	0.4	300.75	398.15	488.81	–3.732e+00	1.055e–02	–1.017e–05
2,4,6-trimethylpyridine	10	0.1	293.15	358.15	16.67	–3.036e–02	2.361e–06	0
dimethoxymethane	5	0.0	170.65	298.15	2.92	–4.106e–03	1.126e–05	0
1,3-dichloropropane	5	0.1	298.15	333.15	–61.39	4.818e–01	–8.107e–04	0
methylsulfanylmethane	6	0.2	273.30	310.48	23.41	–8.896e–02	1.076e–04	0
1,2-ethanedithiol	3	0.0	293.15	333.15	11.23	–1.350e–02	0.000e+00	0

Tmin and Tmax (K) indicate the validity range of the parameterization. N indicates the number of points in the fit; χ2 is the root mean square deviation. See the Supporting Information for details.

The dielectric constant ε(0) (Figure 4, Table S5) appears to be the most difficult property to reproduce in our simulations, with slopes a < 0.5 and R2 ≤ 60% for both force fields (Table 2). Apart from lacking explicit polarization, limited sampling (1000 molecules for 10 ns were used in all cases) may be one of the causes; another contributing factor is the high viscosity of molecules containing alcohol or amine groups, further aggravated by the fact that some of these molecules were simulated at temperatures close to the melting temperature.

Figure 4

Correlation between dielectric constant (ε(0)) calculated using the GAFF and the OPLS/AA force fields and experimental results. Note the logarithmic axes.

Some liquids have extremely large dielectric constants, e.g., methanamide (ε(0) = 109) and N-methylformamide (ε(0) = 190). For these molecules, GAFF predicts 41 and 14, respectively, while OPLS/AA predicts 51 and 19. Xie et al. report a simulated dielectric constant of 200 for N-methylformamide, using a polarizable model, with only 256 molecules and 1 ns of simulation, but the authors state that “The dielectric constants have only been averaged for 1 ns of simulation time, and they are almost certain not yet converged.” Indeed, Whitfield et al. had previously concluded that very long times (50 ns) may be needed to obtain converged dielectric constants of molecules like N-methylacetamide because they tend to form long linear chains.[78] Such chains can in periodic simulation systems become “infinite”, which may contribute to the long relaxation time. It should be noted, however, that for most molecules in our study, the values are well converged, as witnessed by small error bars. Deviations from experimental results are therefore due predominantly to a lack of polarization and too low mobility of molecules. Interestingly, GAFF is somewhat better at predicting ε(0) than OPLS/AA (Table 2), most likely because the partial charges are somewhat higher for most molecules.

Volumetric Expansion Coefficient

The volumetric expansion coefficient αP is plotted in Figure 5 and tabulated in Table S6. The slope of the correlation plots is slightly less than 1 for both GAFF (a = 0.9) and OPLS/AA (a = 0.91), and there is a large spread around the y = x line for both OPLS/AA and GAFF with a RMSD of 0.3/GPa in both cases.

Figure 5

Correlation between volumetric expansion coefficient (αP) calculated using the GAFF and the OPLS/AA force fields and experimental results.

Isothermal Compressibility

For 53 molecules, an interpolation of experimental values of the isothermal compressibility κT as a function of the temperature was performed (Table 4 and Figure S3). The simulated κT’s are plotted versus the experimental values in Figure 6 and tabulated in Table S8. Like for αP, the spread in numbers is large, and the slope of the correlation plots is significantly less than 1 (GAFF, 0.66; OPLS/AA, 0.76, Table 2). In general, it seems that fluctuation properties are more difficult to predict than simple linear averages. Although we applied a very strict convergence criterion for the total energy of 0.5 J/mol/ns per degree of freedom, it may be that even longer equilibration times and production simulations are needed.

Table 4

Parameterization of Temperature Dependence of Isothermal Compressibility Constants in a Polynomial Form κT = A + BT + CT2a

molecule	N	χ2	Tmin	Tmax	A	B	C
dichloromethane	3	0.000	293.15	303.15	–1.709e+01	1.144e–01	–1.800e–04
methanamide	5	0.008	288.15	323.15	1.352e–01	9.161e–04	0
nitromethane	4	0.020	298.15	323.15	–1.253e+00	6.666e–03	0
methanol	24	0.014	213.15	333.15	1.004e+00	–6.791e–03	2.557e–05
acetonitrile	5	0.000	298.15	318.15	3.174e+00	–2.209e–02	5.114e–05
1,1,2,2-tetrachloroethane	2	0.000	293.15	303.15	–4.962e–01	3.900e–03	0
1,1,2-trichloroethane	7	0.002	288.15	318.15	–7.213e–01	4.937e–03	0
bromoethane	5	0.010	273.15	323.15	9.748e+00	–6.685e–02	1.287e–04
N-methylformamide	4	0.011	288.15	313.15	6.378e–03	1.968e–03	0
nitroethane	3	0.015	298.15	323.15	–9.873e–01	6.004e–03	0
ethanol	16	0.007	203.15	363.15	1.280e+00	–8.946e–03	2.857e–05
methylsulfinylmethane	7	0.030	293.15	353.15	5.206e–01	–3.136e–03	1.052e–05
2-aminoethanol	6	0.000	278.15	333.15	7.273e–01	–4.276e–03	1.051e–05
1,3-dichloropropane	6	0.000	283.15	323.15	6.932e–01	–4.785e–03	1.678e–05
propan-2-one	10	0.010	293.15	328.15	–3.053e+00	1.468e–02	0
methyl acetate	8	0.012	293.15	328.15	–2.562e+00	1.249e–02	0
1,3-dioxolane	2	0.000	293.15	313.15	–1.317e+00	6.960e–03	0
1-bromopropane	7	0.003	288.15	318.15	–1.264e+00	8.037e–03	0
N,N-dimethylformamide	18	0.018	288.15	333.20	1.748e+00	–1.073e–02	2.367e–05
1-nitropropane	3	0.004	298.15	323.15	–1.111e+00	6.420e–03	0
2-nitropropane	3	0.020	298.15	323.15	–1.060e+00	6.604e–03	0
1,4-dichlorobutane	5	0.004	288.15	318.15	–8.725e–01	5.246e–03	0
propane-1,2,3-triol	19	0.003	293.15	473.15	8.358e–01	–4.323e–03	7.862e–06
propan-1-amine	6	0.036	293.15	323.15	–2.469e+00	1.238e–02	0
N,N-dimethylacetamide	5	0.015	288.15	318.15	–5.890e–01	4.142e–03	0
butan-1-ol	15	0.021	293.15	393.15	1.307e+00	–8.833e–03	2.543e–05
N-ethylethanamine	5	0.002	298.15	318.15	7.548e+00	–5.536e–02	1.188e–04
butan-1-amine	8	0.003	298.15	328.15	2.330e+00	–1.702e–02	4.371e–05
ethyl acetate	9	0.012	298.15	350.30	5.084e+00	–3.567e–02	7.598e–05
oxolane	5	0.001	278.15	323.15	–9.434e–01	4.999e–03	4.886e–06
1-bromobutane	12	0.000	298.15	333.15	2.650e+00	–1.860e–02	4.413e–05
1-chlorobutane	10	0.029	293.15	318.15	–2.399e+00	1.205e–02	0
pentanenitrile	5	0.005	283.15	323.15	8.811e–01	–7.004e–03	2.429e–05
ethyl propanoate	15	0.022	278.15	338.15	6.964e–01	–7.128e–03	2.882e–05
2-methylbutan-2-ol	2	0.000	293.15	298.15	–1.495e+00	8.600e–03	0
pentan-1-ol	8	0.010	293.15	333.15	3.158e+00	–2.044e–02	4.292e–05
pentan-3-ol	10	0.003	293.15	368.15	4.952e+00	–3.315e–02	6.587e–05
nitrobenzene	5	0.009	298.15	323.15	–3.337e–01	2.832e–03	0
cyclohexanone	5	0.021	298.15	308.15	–9.399e–01	5.421e–03	0
hexan-2-one	8	0.022	278.15	338.15	–1.451e+00	8.315e–03	0
1-methoxy-2-(2-methoxyethoxy)ethane	6	0.001	298.15	318.15	–8.794e–01	5.105e–03	0
N,N-diethylethanamine	8	0.006	298.15	328.15	4.400e+00	–3.405e–02	8.064e–05
N-propan-2-ylpropan-2-amine	7	0.001	298.15	328.15	9.459e+00	–6.732e–02	1.357e–04
methoxybenzene	5	0.043	298.15	338.15	–1.520e+00	7.287e–03	0
3-methylphenol	6	0.041	298.15	413.15	1.744e+00	–1.029e–02	2.104e–05
toluene	50	0.006	288.15	333.15	2.342e+00	–1.627e–02	3.853e–05
diethyl propanedioate	7	0.000	298.15	328.15	2.164e+00	–1.397e–02	3.048e–05
heptan-2-one	2	0.000	293.15	298.15	–8.915e–01	6.200e–03	0
ethylbenzene	7	0.008	293.15	333.15	2.524e+00	–1.652e–02	3.683e–05
1,2-dimethylbenzene	10	0.022	273.15	417.50	–2.914e–01	1.846e–03	6.429e–06
octan-1-ol	16	0.033	293.15	413.15	2.242e+00	–1.449e–02	3.206e–05
quinoline	2	0.000	333.15	373.15	–5.477e–01	3.320e–03	0
(1-methylethyl)benzene	3	0.003	293.15	298.15	–6.340e–01	5.110e–03	0

Tmin and Tmax (K) indicate the validity range of the parameterization. N indicates the number of points in the fit; χ2 is the root mean square deviation. See the Supporting Information for details.

Figure 6

Correlation between isothermal compressibility (κT) calculated using the GAFF and the OPLS/AA force fields and experimental results.

Heat Capacities

For three molecules, an interpolation of experimental numbers is presented in Table 5 and Figure S4. The heat capacity is a difficult property to calculate due to significant quantum effects. The simple eq 7 produces numbers (cPclass) that are twice too high (Table 2). Since the energy taken up by vibrations in a classical harmonic oscillator is much higher than for a quantum harmonic oscillator at the same frequency, the cPclass values are too high. Introducing quantum corrections, in the manner proposed by Berens et al.,[71] on which the two phase thermodynamics (2PT) method[68−70] is based, presupposes that the frequencies in the classical simulation are correct: this is often the case since most force constants have been derived from spectroscopic experiments. It should be noted that there is no a priori reason to assume that the intermolecular degrees of freedom behave harmonically, as they are determined by Coulomb and van der Waals interactions. Despite these theoretical shortcomings, the 2PT method produces reasonable results for cP (see Figure 7, Table 2, and Table S8)—much closer to experiment than cPclass on any account. In order to compute cP, it is necessary to add a correction Δc (eq 12) to the heat capacity at constant volume cV that is produced by the density of states analysis. Δc is underestimated by classical force field calculations; however, cP still is estimated reasonably, with a = 1.08 for GAFF and a = 1.02 for OPLS/AA with correlation coefficients R2 = 98% and 97%, respectively. If we compare just cV from our simulations (i.e., without adding in Δc) and subtract the experimental Δc from the measured cP, we find a very good correlation (GAFF, a = 1.02, R2 = 97%; OPLS/AA, a = 1.04, R2 = 95%), see Figure 8 and Table 2. Although correlation between experimental results and calculations can by no means validate the underlying theoretical model, it nevertheless indicates that the results are meaningful, because we have approximately 70 experimental cV values to which to compare. Indeed, although the root-mean-square deviation (RMSD) from experimental results is similar for cV and cP, the fit to experimental results is much better (slope a close to 1) for both OPLS/AA and GAFF. The DOS simulations were performed without constraints, and the heat capacities depend directly on the intra- and intermolecular vibrations. Deviations from the experimental heat capacities could therefore indicate problems with the force constants for intramolecular motions.

Table 5

Parameterization of Temperature Dependence of Heat Capacity at Constant Pressure in a Polynomial Form cP = A + BTa

molecule	N	χ2	Tmin	Tmax	A	B
1,3-dioxolane	9	0.187	288.15	328.15	4.371e+01	2.613e–01
1,2,3,4-tetrafluorobenzene	41	0.145	235.47	319.79	1.158e+02	2.491e–01
1,2,3,5-tetrafluorobenzene	25	0.343	229.32	311.18	1.186e+02	2.400e–01

Tmin and Tmax (K) indicate the validity range of the parameterization. N indicates the number of points in the fit; χ2 is the root mean square deviation. See the Supporting Information for details.

Figure 7

Correlation between measured heat capacity at constant pressure (cP) and computed using the GAFF and the OPLS/AA force fields based on either the density of states (DoS) method, which includes quantum corrections and a Δc correction based on simulations, or based on a purely classical treatment (cPclass, Class.).

Outliers Per Force Field

Table 6 shows how the molecular models of the individual molecules perform relative to the force field as a whole. The average relative deviations in σ and averaged over 1–8 data points (depending on the availability of experimental data) signals how well the force field performs for each molecule. The properties used were density, enthalpy of vaporization, surface tension, dielectric constant, volumetric expansion coefficient, isothermal compressibility, and the heat capacity at constant volume.

Table 6

Average Relative Deviation (σ) from Experimental Values, in Brackets, the Number of Observablesa

name	CGenFF	GAFF	OPLS/AA
1. chloroform		2.1(6)	3.0(7)
2. dichloro(fluoro)methane		1.0(4)	1.3(4)
3. dibromomethane		2.9(6)	1.7(7)
4. dichloromethane		1.7(7)	3.6(7)
5. methanal		0.3(4)	0.3(4)
6. methanoic acid		4.5(6)	2.6(7)
7. bromomethane		1.4(3)	0.4(3)
8. methanamide	0.0(1)	1.2(7)	0.4(6)
9. nitromethane		2.0(7)	0.8(7)
10. methanol	0.0(2)	0.8(7)	0.8(7)
11. 1,1,1,2,2-pentachloroethane		0.5(4)	0.8(4)
12. 1,1,2,2-tetrachloroethane		1.7(7)	1.7(7)
13. 1,1-dichloroethene		1.7(4)	0.8(4)
14. 1,1,2-trichloroethane		1.2(7)	0.9(7)
15. acetonitrile	0.0(1)	1.1(7)	2.2(7)
16. 1,2-dibromoethane		2.6(7)	4.0(7)
17. 1,1-dichloroethane	0.0(1)	0.7(7)	1.7(7)
18. 1,2-dichloroethane		1.6(7)	1.2(7)
19. methyl formate		0.9(4)	0.8(5)
20. bromoethane	0.0(1)	2.2(7)	0.6(7)
21. chloroethane	0.0(1)	0.8(5)	1.3(5)
22. 2-chloroethanol		0.4(4)	0.5(4)
23. ethanamide		0.2(4)	0.8(5)
24. N-methylformamide		1.4(7)	1.4(7)
25. nitroethane		1.5(7)	0.7(7)
26. methoxymethane		0.5(5)	1.3(5)
27. ethanol	0.0(2)	1.0(7)	0.7(6)
28. 1,2-ethanedithiol		0.6(3)	0.1(3)
29. methyldisulfanylmethane	0.1(2)	1.2(5)	1.6(5)
30. methylsulfinylmethane	0.1(1)	1.0(7)	0.6(7)
31. methylsulfanylmethane		1.4(5)	1.2(5)
32. 2-aminoethanol		1.2(5)	1.3(6)
33. ethane-1,2-diamine		1.2(7)	1.9(7)
34. prop-2-enenitrile		1.0(5)	1.2(5)
35. 1,3-dioxolan-2-one		0.5(5)	0.2(4)
36. propanenitrile		1.1(7)	1.9(7)
37. 1,2-dibromopropane		1.1(5)	0.6(4)
38. 1,3-dichloropropane		0.9(7)	1.0(7)
39. (2R)-2-methyloxirane		0.0(2)	0.1(2)
40. propan-2-one	0.0(2)	1.0(7)	0.7(7)
41. methyl acetate	0.0(2)	1.3(7)	0.9(7)
42. 1,3-dioxolane	0.0(1)	1.2(4)	0.6(4)
43. 2-iodopropane		0.7(5)	1.1(5)
44. 1-bromopropane		1.3(7)	0.6(7)
45. N,N-dimethylformamide		0.7(6)	0.5(6)
46. N-methylacetamide	0.0(1)	0.4(4)	0.2(4)
47. 1-nitropropane		1.6(7)	1.2(7)
48. 2-nitropropane		1.6(7)	0.9(7)
49. dimethoxymethane		0.8(5)	0.9(5)
50. propane-1,2,3-triol		1.3(6)	0.8(6)
51. propan-1-amine		1.1(7)	1.5(7)
52. propan-2-amine		0.7(5)	0.6(4)
53. 2-methylpropane	0.0(1)	0.8(5)	1.1(5)
54. ethylsulfanylethane		0.6(5)	0.7(5)
55. butane-1-thiol		0.9(5)	0.5(5)
56. butan-1-ol		1.1(7)	0.9(7)
57. 2-methylpropan-2-ol		0.4(2)	0.1(2)
58. butane-1,4-diol		0.9(6)	0.4(6)
59. (2-hydroxyethoxy)ethan-2-ol		1.2(4)	1.1(5)
60. N-ethylethanamine		1.1(7)	1.2(7)
61. butan-1-amine		1.1(7)	0.9(7)
62. 2-methylpropan-2-amine		1.0(5)	0.8(5)
63. 2-(2-hydroxyethylamino)ethanol		0.5(4)	0.4(4)
64. pyrimidine	0.0(2)	0.7(4)	0.6(4)
65. furan	0.2(2)	1.9(5)	1.9(5)
66. thiophene	0.0(2)	0.7(4)	0.3(5)
67. 1H-pyrrole	0.1(1)	1.3(7)	1.1(7)
68. ethenyl acetate		0.5(4)	0.8(4)
69. oxolan-2-one		0.3(3)	0.3(4)
70. acetyl acetate		1.2(4)	1.2(4)
71. 1,4-dichlorobutane		0.6(7)	0.8(7)
72. oxolane		0.6(6)	1.3(7)
73. ethoxyethene		0.3(3)	0.2(3)
74. ethyl acetate	0.0(2)	1.2(7)	1.1(7)
75. tetrahydrothiophene 1,1-dioxide		0.8(4)	0.9(4)
76. thiolane		0.5(4)	0.4(4)
77. 1-bromobutane		1.1(7)	0.7(7)
78. 1-chlorobutane		1.4(7)	1.8(7)
79. pyrrolidine	0.1(1)	1.3(7)	1.3(7)
80. N,N-dimethylacetamide		1.0(7)	0.9(7)
81. morpholine		0.8(5)	0.9(5)
82. pyridine	0.1(2)	0.6(6)	0.9(7)
83. cyclopentanone		0.8(5)	0.6(5)
84. 1-cyclopropylethanone		0.2(2)	0.1(2)
85. pentane-2,4-dione		0.9(5)	1.3(5)
86. methyl 2-methylprop-2-enoate		0.8(5)	3.6(5)
87. pentanenitrile		0.6(6)	1.6(7)
88. ethyl propanoate		1.3(7)	1.1(7)
89. diethyl carbonate		2.1(7)	0.7(6)
90. pentan-1-ol		1.0(7)	0.9(7)
91. pentan-3-ol		1.0(7)	1.1(7)
92. 2-methylbutan-2-ol		1.1(5)	0.5(5)
93. pentane-1,5-diol		0.8(6)	0.6(6)
94. pentan-3-amine		0.5(4)	0.6(4)
95. 1,2,3,4-tetrafluorobenzene		0.2(2)	0.1(2)
96. 1,2,3,5-tetrafluorobenzene		0.2(2)	0.1(2)
97. 1,3-difluorobenzene	0.2(2)	0.7(4)	1.3(5)
98. 1,2-difluorobenzene		0.7(4)	1.0(5)
99. fluorobenzene	0.1(2)	1.6(7)	0.5(6)
100. nitrobenzene	0.0(2)	1.1(7)	1.1(7)
101. 2-chloroaniline		0.9(4)	0.6(4)
102. phenol		0.8(4)	0.9(5)
103. benzenethiol		1.4(5)	1.3(5)
104. 2-methylpyridine		0.3(4)	0.9(5)
105. 3-methylpyridine	0.1(2)	0.8(5)	0.6(5)
106. 4-methylpyridine	0.0(2)	1.1(7)	0.4(6)
107. cyclohexanone		1.0(7)	0.9(7)
108. (E)-hex-2-ene	0.0(2)	0.0(2)	0.0(2)
109. hexan-2-one		0.8(6)	0.9(7)
110. 2,4,6-trimethyl-1,3,5-trioxane		1.6(4)	1.0(4)
111. cyclohexanamine		0.8(5)	0.7(5)
112. 2-propan-2-yloxypropane		3.3(7)	0.9(7)
113. 1-methoxy-2-(2-methoxyethoxy)ethane		1.5(7)	1.1(7)
114. triethyl phosphate		2.8(6)	2.2(6)
115. N,N-diethylethanamine		1.2(7)	1.0(7)
116. N-propan-2-ylpropan-2-amine		0.8(6)	0.6(6)
117. trifluoromethylbenzene		0.8(5)	0.5(4)
118. benzonitrile		1.0(5)	1.0(5)
119. benzaldehyde	0.2(2)	5.7(7)	3.7(6)
120. toluene	0.1(2)	1.6(7)	1.3(7)
121. methoxybenzene	0.1(2)	1.2(7)	1.1(7)
122. phenylmethanol		1.0(5)	0.8(5)
123. 2-methylphenol		0.9(5)	0.8(5)
124. 3-methylphenol		1.0(5)	0.9(5)
125. 4-methylphenol	0.1(1)	1.2(5)	0.7(5)
126. diethyl propanedioate		1.1(4)	0.8(4)
127. 2,4-dimethylpentan-3-one		0.6(4)	0.4(4)
128. heptan-2-one		1.1(7)	0.7(7)
129. ethenylbenzene		1.2(5)	1.1(5)
130. 1-phenylethanone		1.0(7)	1.1(7)
131. methyl benzoate		0.9(7)	1.0(7)
132. methyl 2-hydroxybenzoate		1.1(5)	0.4(4)
133. ethylbenzene	0.1(2)	1.4(7)	1.1(7)
134. 1,2-dimethylbenzene	0.1(1)	1.7(7)	1.0(7)
135. 1,2-dimethoxybenzene		0.4(4)	0.6(5)
136. 2,4,6-trimethylpyridine		0.9(5)	1.0(5)
137. octan-1-ol		0.8(6)	1.7(7)
138. 1-butoxybutane		0.7(4)	1.0(5)
139. N-butylbutan-1-amine		0.9(7)	0.8(7)
140. isoquinoline	0.0(1)	0.7(4)	1.3(4)
141. quinoline	0.1(2)	1.1(7)	1.2(7)
142. (1-methylethyl)benzene	0.1(2)	1.0(6)	0.7(6)
143. 1,2,4-trimethylbenzene	0.1(1)	1.2(6)	1.0(6)
144. 2,6-dimethylheptan-4-one		1.0(5)	0.9(5)
145. 1-chloronaphthalene		0.5(6)	1.3(7)
146. phenoxybenzene		0.6(4)	1.1(5)

Average relative deviation larger than 1σ is printed in bold, larger than 1.5σ in bold italic.

Some types of molecules are problematic in both of the force fields considered here. Small molecules containing more than one Cl or Br atom generally have both density and enthalpy of vaporization values that deviate significantly from experimental reference. This is not the case for molecules containing only one of these atoms, or molecules where there is a spacer (e.g., a CH2 group) between them. It could therefore be that the differences are caused by overlapping atoms. By introducing a new atom type of Br and Cl for the case where there are two such atoms next to each other on the carbon chain, these problems might be resolved.

The density and enthalpy of vaporization of methanoic acid (formic acid) were particularly hard to reproduce, as was noted previously by Jedlovsky and Turi, who constructed a specific potential for this molecule.[79] The main feature responsible for the improved model in this case was a higher charge (≈ 0.1e) on the C–H atom than is used in either OPLS/AA (0) or GAFF (0.04). Methanoic acid forms very strong linear chains, which are difficult to break. This leads to long correlation times for the system dipoles and to dielectric constants that are far from the experimental values (Table S5).

Benzaldehyde and furan are also problematic in both force fields. Even if they both generate decent densities and enthalpies of vaporization, the other properties (surface tension, dielectric constant, and thermal expansion coefficient) are far from the experimental values.

Molecules containing a nitro group (specially nitromethane, 1-nitropropane, and 2-nitropropane) stick out as a problematic group in GAFF. The charges on nitro groups are high, leading to high density and enthalpy of vaporization.

The standard OPLS/AA parametrization of alcohols has been reported to perform poorly for octan-1-ol. MacCallum and Tieleman[80] therefore derived a specific united atom potential of the molecule where they used modified charges on the headgroup. The OPLS/AA parametrization investigated here gives both too high a density and too high an enthalpy of vaporization, and therefore the other properties investigated for this molecules also deviate from experimental results. Methyl-2-methylprop-2-enaote shows similar problems, and this could probably be corrected in a similar way. It should be noted that, compared to GAFF, the charges on the headgroup in these two molecules are relatively high in OPLS/AA.

Discussion

The development of force fields for molecular simulation is critically dependent on the availability of good reference data, preferably from experimental sources. All force fields, be they empirical, purely derived from quantum-mechanics, or a combination of the two, will eventually have to face the test of comparing predicted to measured values. There is a large amount of literature on force field testing for proteins and peptides,[57,81−87] nucleic acids,[88−91] carbohydrates,[92] specific organic molecules or protein fragments,[20,26,93−97] and ions,[98−101] to list but a few. In addition, there are indirect force field tests, for instance of the binding energy in protein–ligand complexes,[16,102] protein structure prediction,[103] or of force-field-based docking codes.[104−106] It is interesting to mention the industrial fluid properties simulation challenges, which are stimulating modelers to predict properties of liquids by any means, including molecular simulation.[107,108]

Here, we have introduced a benchmark set of 146 liquids in order to assess two popular all atom force fields, OPLS/AA and GAFF, and to set a standard for future force fields. For comparison, we have included an independent density and enthalpy of vaporization data set computed using CGenFF, based on a similar set of molecules.[37] Calculated density, enthalpy of vaporization, heat capacities, surface tension, dielectric constants, volumetric expansion coefficients, and isothermal compressibility from the two force fields are compared to experimental values. Indeed the benchmark is quite revealing, in that systematic deviations can be found and rationalized. The knowledge about such deviations will hopefully be useful for further development of the force fields.

To a first approximation, molecular vibrations can be described as quantum harmonic oscillators.[109] Classical harmonic oscillators do not describe the properties of quantum harmonic oscillators, which makes it necessary to implement corrections in computing for instance heat capacities. The two phase thermodynamics method employed here for estimating cP and cV relies on the force constants of the force field used, and on the effective frequencies in the simulations. The density of states obtained from the velocity autocorrelation is convoluted by a weighting function derived from the partition function for a quantum harmonic oscillator in order to obtain a heat capacity for a corresponding quantum liquid. If a force field would allow one to directly reproduce the “correct” density of states, one could use the much simpler fluctuation formulas, as described by Allen and Tildesley;[38] however, heat capacities computed in this manner overestimate the experimental values by about 100% for OPLS/AA and GAFF (Table S10). Going beyond the harmonic approximation should therefore be considered by force field developers. Despite efforts in the context of the MMF94 force field[110] and the MM3-MM4 family of force fields,[111−113] this has not been widely adopted in the biomolecular simulation community, although the polarizable AMOEBA force field[114] does feature anharmonic bond and angle potentials as well. In principle, it should be advantageous to use for instance Car–Parrinello molecular dynamics,[115] in order to more faithfully represent a liquid than is possible in a classical simulation. This was attempted by Kuo et al. for water.[116] They find a large scatter in cP values due to limited sampling, but also a systematic deviation from the experimental value. Obviously, the computational bottleneck that would be introduced by CPMD or related methods will remain difficult to surmount for the immediate future, and therefore force-field-based methods remain necessary. Nevertheless, it is encouraging that there is a trend to use molecular dynamics simulations based on density functional theory codes to study vibrational properties of biomolecular systems beyond the harmonic approximation.[117−119]

The dielectric constant seems to be the hardest nut to crack. Nonpolarizable force fields (such as GAFF and OPLS/AA) are known to have difficulties in reproducing the dielectric constant and to some extent also the surface tension. In the case of water, for which a large number of force fields have been developed, there are several studies that describe this (for a review, see, for example, Guillot[25]). Improving the dielectric function often turns out to be done at the cost of the enthalpy of vaporization and the free energy of solvation—properties that may be more important to reproduce in biomolecular simulations. In addition to systematic problems, like sampling or the lack of polarization in our simulations,[120] the temperature dependence of the dielectric constant provides both a challenge and an opportunity for future force field development. For most molecules, the temperature dependence is very strong, because molecular motion is the largest factor contributing to ε(0). In his review of water models, Guillot has pointed out that the relation between dielectric constant and other properties is complex, and hence it can be used to test and validate force fields, but not likely as a target for force field optimization.[25]

The benchmark we present here allows one to pinpoint systematic errors in force fields due to the fact that most chemical moieties are represented more than once. The overall performance of GAFF is surprisingly good, seeing that the parameter development was not aimed at liquids. The results from the OPLS/AA force field are slightly better than GAFF, obviously due to the fact that OPLS/AA was parametrized for liquids. The CHARMM generalized force field seems to be even slightly better, at least for density and enthalpy of vaporization.[37] It is reassuring for applications of force field calculations beyond liquids that the parameters in most cases are reasonable; however, the results presented here also show that blind faith in force fields is not warranted in all cases. In Table 2, we list the root-mean-square deviation, as well as the average relative deviation, of the calculated values from the experimental, for each property we have analyzed. Even if our set of molecules is limited to 146, these numbers give a measurement of how well the properties are reproduced in the two force fields, at least for molecules similar to the set presented here.

Wang and Tingjun have recently reported a similar force field test of 71 organic molecules based on the GAFF and OPLS/AA force fields.[32] They report densities and enthalpies of vaporization for these molecules and find small deviations from experimental results that are comparable to our numbers. It is encouraging to note that these authors were able to improve the correspondence to experimental numbers by tuning the Lennard-Jones parameters of some of the atom types. How this affects the other properties that we have studied here, in particular, the dielectric constant and the surface tension, remains to be determined, but it is likely that just tweaking the Lennard-Jones parameters is not sufficient to cure the significant and systematic deviations observed for those properties.

Mobley et al. have performed free energy of solvation (ΔGhyd) benchmarks, reporting a RMS error from experimental numbers of 5.2 kJ/mol for more than 500 molecules.[121,122] This number is comparable to the RMSD of 6.5 kJ/mol we computed for ΔHvap for OPLS/AA (10.6 kJ/mol for GAFF and 4.7 kJ/mol for CGenFF[37]). Since both numbers are to a large extent determined by the intermolecular energies, we can conclude that the RMS error in intermolecular energies for (“small”) organic molecules is 5–6 kJ/mol using state of the art simulations and nonpolarizable force fields. It should be noted that this result may be biased by the choice of test set, as has been shown in the context of the SAMPL contest where hydration energies were to be predicted.[123,124] It was found here that larger molecules with multiple functional groups have similar deviations from the experimental hydration energy—errors up to 10 kJ/mol.[123] It seems plausible that part of this error is due to the simple nonpolarizable water model used, however, since the enthalpy of vaporization is approximately additive (which can be seen by plotting ΔHvap for, e.g., alkanes as a function of the number of carbons), the error per functional group should still be relatively low, less than 5 kJ/mol for most groups. In the present work, we studied pure liquids only, providing a simpler test set than what has been used in previous studies. Further tests on pure liquids and liquid mixtures should provide a more detailed understanding of the predictive power of force field calculations. At the same time, systematic methods for force field development[23,125] could be used for the improvement of classical force fields.