Free energies of hydration from thermodynamic integration: Comparison of molecular mechanics force fields and evaluation of calculation accuracy

Four commonly used molecular mechanics force fields, CHARMM22, OPLS, CVFF, and GROMOS87, are compared for their ability to reproduce experimental free energies of hydration (ΔG_hydr) from molecular dynamics (MD) simulations for a set of small nonpolar and polar organic molecules: propane, cyclopropane, dimethylether, and acetone. ΔG_hydr values were calculated by multiconfiguration thermodynamic integration for each of the different force fields with three different sets of partial atomic charges: full charges from an electrostatic potential fit (ESP), and ESP charges scaled by 0.8 and 0.6. All force fields, except for GROMOS87, give reasonable results for ΔG_hydr · if partial atomic charges of appropriate magnitude are assigned. For GROMOS87, the agreement with experiment for hydrocarbons (propane and ethane) was improved considerably by modifying the repulsive part of the carbon-water oxygen Lennard-Jones potential. The small molecules studied are related to the chemical moieties constituting camphor (C₁₀H_l6O). By invoking force-field transferability, we calculated the ΔG_hydr for camphor. With the modified GROMOS force field, a ΔG_hydr within 4 kJ/mol of the experimental value of −14.8 kJ/mol was obtained. Camphor is one of the largest molecules for which an absolute hydration free energy has been calculated by molecular simulation. The accuracy and reliability of the thermodynamic integration calculations were analyzed in detail and we found that, for ΔG_hydr calculations for the set of small molecules in aqueous solution, molecular dynamics simulations of 0.8–1.0 ns in length give an upper statistical error bound of 1.5 kJ/mol, whereas shorter simulations of 0.25 nm in length given an upper statistical error bound of 3.5 kJ/mol. © 1997 by John Wiley & Sons, Inc.     

An Electrostatic Bond Energy Model for Enthalpies of Formation of Chloroalkanes in Gaseous States

An electrostatic bond energy model is formulated to fit the enthalpies of formation and dipole moments of the alkanes and chloroalkanes. In this model, the charge distributions are calculated by an electrostatic approach similar to the "MSE" method, and the enthalpy of formation of a molecule is the sum of the bond energy terms plus the electrostatic energy of the interactions between the charges on all atoms. All parameters of this model are obtained by parameterization. The calculated dipole moments for 13 chloroalkanes and enthalpies of formation for 19 alkanes and non-geminal chloroalkanes agree with the determined values very well. To calculate the enthalpies of formation of geminal chloroalkanes, a correction mainly attributed to the van der Waals interactions in the geminal substituted group, about 24 kJ/mol per pair of geminal chlorine atoms, is introduced.     

Prediction of SAMPL3 host–guest binding affinities: evaluating the accuracy of generalized force-fields

We used the second-generation mining minima method (M2) to compute the binding affinities of the novel host-guest complexes in the SAMPL3 blind prediction challenge. The predictions were in poor agreement with experiment, and we conjectured that much of the error might derive from the force field, CHARMm with Vcharge charges. Repeating the calculations with other generalized force-fields led to no significant improvement, and we observed that the predicted affinities were highly sensitive to the choice of force-field. We therefore embarked on a systematic evaluation of a set of generalized force fields, based upon comparisons with PM6-DH2, a fast yet accurate semi-empirical quantum mechanics method. In particular, we compared gas-phase interaction energies and entropies for the host-guest complexes themselves, as well as for smaller chemical fragments derived from the same molecules. The mean deviations of the force field interaction energies from the quantum results were greater than 3 kcal/mol and 9 kcal/mol, for the fragments and host-guest systems respectively. We further evaluated the accuracy of force-fields for computing the vibrational entropies and found the mean errors to be greater than 4 kcal/mol. Given these errors in energy and entropy, it is not surprising in retrospect that the predicted binding affinities deviated from the experiment by several kcal/mol. These results emphasize the need for improvements in generalized force-fields and also highlight the importance of systematic evaluation of force-field parameters prior to evaluating different free-energy methods.     

Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models

We present an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains at 298 K and 1 atm. The hydration free energies DeltaG, entropies DeltaS, enthalpies DeltaH, and heat capacities Deltac(P)() were determined for 10 combinations of force fields and water models. The statistical sampling was extended such that precisions of 0.3, 0.8, 0.8 kJ/mol and 25 J/(mol K) were reached for DeltaG, TDeltaS, DeltaH, and Deltac(P)(), respectively. The three force fields used in this study are AMBER99, GROMOS 53A6, and OPLS-AA; the five water models are SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew. We found that the choice of water model strongly influences the accuracy of the calculated hydration entropies, enthalpies, and heat capacities, while differences in accuracy between the force fields are small. On the basis of an analysis of the hydrophobic analogues of the amino acid side chains, we discuss what properties of the water models are responsible for the observed discrepancies between computed and experimental values. The SPC/E water model performs best with all three biomolecular force fields.     

A molecular mechanics force field for biologically important sterols

A parameterization has been performed of the biologically important sterols cholesterol, ergosterol, and lanosterol for the CHARMM27 all-atom molecular mechanics force field. An automated parameterization method was used that involves fitting the potential to vibrational frequencies and eigenvectors derived from quantum-chemical calculations. The partial charges were derived by fitting point charges to quantum-chemically calculated electrostatic potentials. To model the dynamics of the hydroxyl groups of the sterols correctly, the parameter set was refined to reproduce the energy barrier for the rotation of the hydroxyl group around the carbon connected to the hydroxyl of each sterol. The frequency-matching plots show good agreement between the CHARMM and quantum chemical normal modes. The parameters are tested in a molecular dynamics simulation of the cholesterol crystal structure. The experimental geometry and cell dimensions are well reproduced. The force field derived here is also useful for simulating other sterols such as the phytosterols sigmasterol, and campesterol, and a variety of steroids.     

Atom equivalents for converting DFT energies calculated on molecular mechanics structures to formation enthalpies

New atom equivalents are introduced to convert BP/DN**//MMFF energies into formation enthalpies. As a result of using molecular mechanics structures, poor results are obtained for compounds outside the scope of the force field, such as those bearing  NF₂ groups or some nitrogenous systems. Notwithstanding these limitations, present procedures compare well with the results of previous atom equivalents schemes. Indeed, rms deviations from experiment are below 9 kJ/mol for hydrocarbons, and close to 16 kJ/mol for a variety of compounds reasonably well described by MMFF. The explicit inclusion of thermal and vibrational contributions, using calculated frequencies, does not improve the results. This study demonstrates that cost‐effective approaches to formation enthalpies may be developed on the basis of a combination of DFT with a suitable molecular mechanics force field. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 367–379, 2000     

Is the enthalpy of fusion of tris(acetylacetonato)metal(III) complexes affected by their potential energy in the crystal state?

In this Article we present enthalpies of fusion and melting points obtained from new thermochemical measurements of tris(acetylacetonato)metal(III), M(acac)(3), complexes (M = Fe, Al, Cr, Mn, Co) using differential scanning calorimetry (DSC) and evaluate them in relation to their different values found in the literature. An enthalpy of fusion of 27.67 kJ mol(-)(1) was derived for Mn(acac)(3) from a symmetrical DSC thermogram captured for the first time. The enthalpy value was indirectly confirmed with the solubility measurements of Mn(acac)(3) in acetylacetone. A hypothesis has been stated that the enthalpy of fusion and the potential energy of M(acac)(3) in the crystal state may be related. To calculate molecular in-crystal potential energy, in this Article we proposed a molecular mechanics model for the M(acac)(3) class of compounds. Nine X-ray crystal structures of M(acac)(3) complexes (M = Fe, Al, V, Mn, Co, Cr, Sc) were included in the modeling. The conformational potential energy was minimized for a molecule surrounded by other molecules in the crystal lattice. The partial charges from two schemes, the electrostatic potential (ESP) fit and the natural population analysis (NPA), were used to construct two types of force fields to examine which force field type would yield a better fit with the experimental thermal properties. The final force fields were named FF-ESP and FF-NPA. Both force field sets reproduced well the experimental crystal data of nine M(acac)(3) complexes as well as of tris(3-methyl-2,4-pentanedionato-O,O')cobalt(III). Only in-crystal potential energies derived by FF-NPA yielded a significant correlation (correlation coefficient R = -0.71) with the measured enthalpies of fusion. The enthalpy of fusion for Co(acac)(3) could not be determined experimentally because of simultaneous decomposition and fusion, and it is predicted to be 33.2 kJ mol(-)(1) from the correlation regression line.     

Ab-initio molecular geometry and normal coordinate analysis of tetrahydrothiophene molecule.

The molecular geometry of tetrahydrothiophene (THT) was quantum mechanically calculated using the split valence 6-31G** basis set. Electron correlation energy has been computed employing MP2 method. The molecule showed a twist form puckered structure with a twist torsion angle of 13 degrees and has a total energy of -347,877.514 kcal/mol of which a 436.715 kcal/mol electron correlation energy. The envelope form of the molecule showed an inter-plane angle of 22 degrees and has a total energy of -347,874.430 kcal/mol involving -436.558 kcal/mol electron correlation energy. The normal coordinates of the molecule were theoretically analyzed and the fundamental vibrational frequencies were calculated. The IR and laser Raman spectra of THT molecule was measured. All the observed vibrational bands including combination bands and overtones were assigned to normal modes with the aid of the potential energy distribution values obtained from normal coordinate calculations. The molecular force field was determined by refining the initial set of force constants using the least square fit method instead of using the less accurate scaling factor methods. The determined molecular force field has produced simulated frequencies which best match the observed values. The lowest-energy modes of vibration were two molecular out-of-plane deformations, observed at 114 and 166 cm(-1). The barrier of ring twisting estimated from the observed ring out-of-plane vibrational mode at 114 cm(-1) was estimated.     

Ab-initio molecular geometry and normal coordinate analysis of pyrrolidine molecule.

The molecular geometry of pyrrolidine was quantum mechanically calculated using the split valence 6-31G** basis set. Electron correlation energy has been computed employing MP2 method. The molecule showed an envelope form puckered structure with inter-plane angle of 36.4 degrees and has a total energy of -132976.80 kcal mol(-1) of which a -464.86 kcal mol(-1) electron correlation energy. The twist form of the molecule showed a twist angle of 10.2 degrees from planarity and has a total energy of -132976.05 kcal mol(-1) involving -464.097 kcal mol(-1) electron correlation energy. The normal coordinates of the molecule were theoretically analyzed on the basis of the Cs point symmetry of the envelope form. Using initial set of force constants obtained from the ab-initio calculations the fundamental vibrational frequencies were computed. The IR and laser Raman spectra of Pyrrolidine molecules were measured. All the observed vibrational bands including combination bands and overtones were assigned to normal modes with the aid of the potential energy distribution values obtained from normal coordinate calculations. The molecular force field was obtained by refining the initial set of force constants using the least square fit method. The molecular force field was determined by refining the initial set of force constants using the least square fit method instead of using the less accurate scaling factor methods. The determined molecular force field has produced simulated frequencies best match to the observed values. The low frequency molecular out-of-plane deformation modes were observed in both infrared and Raman spectra at 298 and 163 cm(-1). The barrier of ring twisting estimated from the observed ring out-of-plane vibrational mode at 163 cm(-1) was found 3.1 kcal mol(-1).     

Computational study of hydrogen binding by metal-organic framework-5

We report the results of quantum chemistry calculations on H(2) binding by the metal-organic framework-5 (MOF)-5. Density functional theory calculations were used to calculate the atomic positions, lattice constant, and effective atomic charges from the electrostatic potential for the MOF-5 crystal structure. Second-order M?ller-Plesset perturbation theory was used to calculate the binding energy of H(2) to benzene and H(2)-1,4-benzenedicarboxylate-H(2). To achieve the necessary accuracy, the large Dunning basis sets aug-cc-pVTZ, and aug-cc-pVQZ were used, and the results were extrapolated to the basis set limit. The binding energy results were 4.77 kJ/mol for benzene, 5.27 kJ/mol for H(2)-1,4-benzenedicarboxylate-H(2). We also estimate binding of 5.38 kJ/mol for Li-1,4-benzenedicarboxylate-Li and 6.86 kJ/mol at the zinc oxide corners using second-order M?ller-Plesset perturbation theory. In order to compare our theoretical calculations to the experimental hydrogen storage results, grand canonical Monte Carlo calculations were performed. The Monte Carlo simulations identify a high energy binding site at the corners that quickly saturated with 1.27 H(2) molecules at 78 K. At 300 K, a broad range of binding sites are observed.     

Molecular mechanics (MM3) calculations on lithium amide compounds

The MM3 force field has been extended to deal with the lithium amide molecules that are widely used as efficient catalysts for stereoselective asymmetric synthesis. The MM3 force field parameters have been determined on the basis of the ab initio MP2/6-31G* and/or DFT (B3LYP/6-31G*, B3-PW91/6-31G*) geometry optimization calculations. To evaluate the electronic interactions specific to the lithium amides derived from the diamine molecules properly, the Lewis bonding potential term for the interaction between the lithium atom and the nonbonded adjacent electronegative atom such as nitrogen was introduced into the MM3 force field. The bond dipoles were evaluated correctly from the electronic charges on the atoms calculated by fitting to the electrostatic potential at points selected. The MM3 results on the molecular structures, conformational energies, and vibrational spectra show good agreement with those from the quantum mechanical calculations.     

Dissolved organic carbon--contaminant interaction descriptors found by 3D force field calculations

Enthalpies of transfer at 300 K of various partitioning processes were calculated in order to study the suitability of 3D force fields for the calculation of partitioning constants. A 3D fulvic acid (FA) model of dissolved organic carbon (DOC) was built in a MM+ force field using AMI atomic charges and geometrical optimization (GO). 3,5-Dichlorobiphenyl (PCB14), 4,4'-dichlorobiphenyl (PCB15), 1,1,1-trichloro-2,2-bis-(4-chlorophenyl)-ethane (PPDDT) and 2-chloro-4-ethylamino-6-isopropylamino-s-triazine (Atrazine) were inserted into different sites and their interaction energies with FA were calculated. Energies of hydration were calculated and subtracted from FA-contaminant interactions of selected sites. The resulting values for the enthalpies of transfer from water to DOC were 2.8, -1.4, -6.4 and 0.0 kcal/mol for PCB 14, PCB15, PPDDT and Atrazine, respectively. The value of PPDDT compared favorably with the experimental value of -5.0 kcal/mol. Prior to this, the method was studied by the calculation of the enthalpies of vaporization and aqueous solution using various force fields. In the MM + force field GO predicted enthalpies of vaporization deviated by +0.7 (PCB14), +3.6 (PCB15) and -0.7 (PPDDT)kcal/mol from experimental data, whereas enthalpies of aqueous solution deviated by -3.6 (PCB14), +5.8 (PCB15) and +3.7 (PPDDT) kcal/mol. Only for PCB14 the wrong sign of this enthalpy value was predicted. Potential advantages and limitations of the approach were discussed.     

Vibrational spectroscopic determination of local solvent electric field, solute-solvent electrostatic interaction energy, and their fluctuation amplitudes

IR probes have been extensively used to monitor local electrostatic and solvation dynamics. Particularly, their vibrational frequencies are highly sensitive to local solvent electric field around an IR probe. Here, we show that the experimentally measured vibrational frequency shifts can be inversely used to determine local electric potential distribution and solute-solvent electrostatic interaction energy. In addition, the upper limits of their fluctuation amplitudes are estimated by using the vibrational bandwidths. Applying this method to fully deuterated N-methylacetamide (NMA) in D(2)O and examining the solvatochromic effects on the amide I' and II' mode frequencies, we found that the solvent electric potential difference between O(═C) and D(-N) atoms of the peptide bond is about 5.4 V, and thus, the approximate solvent electric field produced by surrounding water molecules on the NMA is 172 MV/cm on average if the molecular geometry is taken into account. The solute-solvent electrostatic interaction energy is estimated to be -137 kJ/mol, by considering electric dipole-electric field interaction. Furthermore, their root-mean-square fluctuation amplitudes are as large as 1.6 V, 52 MV/cm, and 41 kJ/mol, respectively. We found that the water electric potential on a peptide bond is spatially nonhomogeneous and that the fluctuation in the electrostatic peptide-water interaction energy is about 10 times larger than the thermal energy at room temperature. This indicates that the peptide-solvent interactions are indeed important for the activation of chemical reactions in aqueous solution.     

聚(对乙烯基苄基己内酰胺)树脂对非水体系中苯酚的吸附热力学

测定了聚（对乙烯基苄基己内酰胺）吸附树脂对正庚烷溶液中苯酚的吸附等温线,发现吸附等温线符合Freundlich等温吸附方程,相关系数大于0．99．根据热力学函数关系计算了等量吸附焓、Gibbs吸附自由能和吸附熵,等量吸附焓在30～40kJ／mol之间,推测吸附过程为氢键吸附．通过比较树脂对正庚烷中苯酚和邻硝基苯酚的吸附性能的差别,进一步论证了树脂对非水体系中酚类物质的吸附是基于氢键作用的机理．同时,聚（对乙烯基苄基己内酰胺）树脂对正庚烷溶液中苯酚的吸附相对水溶液中苯酚的吸附有更低的吸附焓变、吸附自由能变和吸附熵变．     

Quantum mechanical investigation of the inner-sphere reorganization energy of cyclooctatetraene/cyclooctatetraene radical anion. Part I

The inner-sphere reorganization energy of the electron self-exchange of the couple cyclooctatetraene/cyclooctatetraene radical anion has been investigated by quantum mechanical calculations. The more stable Jahn Teller distorted B2g conformation of the radical anion has been used in this study. Two different theories have been applied in this first part. The harmonic approximation in the classical Marcus scheme has been modified by using projected force constants, which are obtained from the complete force constant matrix and the geometry changes of the molecule during the ET (introduced by Mikkelsen). A different approach (introduced by Nelsen) combines the different energies of the neutral and radical anion with and without relaxation corresponding to the vertical ionization potential and the vertical electron affinity. The electronic energies of the neutral molecule and the radical anion differ dramatically applying three different levels of quantum mechanical calculations (UAM1, UB3LYP, PMP2 with three different basis sets with and without diffuse functions). Nevertheless the Nelsen method gives almost consistent results for the inner-sphere reorganization energies: 120.1 kJ/mol for semiempirical UAM1 method, 159.3 kJ/mol, 156.4 kJ/mol and 158.3 kJ/mol for density functional UB3LYP/6-31G*, UB3LYP/6-31++G* and UB3LYP/AUG-cc-pVDZ calculations and 192.5 kJ/mol for ab-initio PMP2/6-31G* investigations, respectively. These values are in agreement with earlier experimental work supposing the total reorganization energy to be larger than 38 kcal/mol assuming an electron self-exchange rate of 10(4) M(-1) s(-1). The simple harmonic approximation of Marcus relation has not yet been applied for a molecule like cyclooctatetraene with large torsional geometry changes. Using the projected force constants after scaling, considerably different results for the inner-sphere reorganization energy have been calculated: 738.1 kJ/mol for the UB3LYP/6-31G*, 743.3 kJ/mol for UB3LYP/6-31++G* and 759.1 kJ/mol for UB3LYP/AUG-cc-pVDZ level of theory. Comparison with our concentration dependent EPR experiments are controversial to the earlier experimental results, but the latter supports the assumption that the electron self-exchange occurs in a time scale so that the molecules cannot complete their vibrational motions. Therefore the projected Marcus relation is not valid for cyclooctatetraene/cyclooctatetraene radical anion including a large torsional change during the electron transfer.     

Dissolved organic carbon--contaminant interaction descriptors found by 3D force field calculations

Enthalpies of transfer at 300 K of various partitioning processes were calculated in order to study the suitability of 3D force fields for the calculation of partitioning constants. A 3D fulvic acid (FA) model of dissolved organic carbon (DOC) was built in a MM+ force field using AM1 atomic charges and geometrical optimization (GO). 3,5-Dichlorobiphenyl (PCB14), 4,4'-dichlorobiphenyl (PCB15), 1,1,1-trichloro-2,2-bis-(4-chlorophenyl)-ethane (PPDDT) and 2-chloro-4-ethylamino-6-isopropylamino- s -triazine (Atrazine) were inserted into different sites and their interaction energies with FA were calculated. Energies of hydration were calculated and subtracted from FA-contaminant interactions of selected sites. The resulting values for the enthalpies of transfer from water to DOC were 2.8, m 1.4, m 6.4 and 0.0 kcal/mol for PCB14, PCB15, PPDDT and Atrazine, respectively. The value of PPDDT compared favorably with the experimental value of m 5.0 kcal/mol. Prior to this, the method was studied by the calculation of the enthalpies of vaporization and aqueous solution using various force fields. In the MM+ force field GO predicted enthalpies of vaporization deviated by +0.7 (PCB14), +3.6 (PCB15) and m 0.7 (PPDDT) kcal/mol from experimental data, whereas enthalpies of aqueous solution deviated by m 3.6 (PCB14), +5.8 (PCB15) and +3.7 (PPDDT) kcal/mol. Only for PCB14 the wrong sign of this enthalpy value was predicted. Potential advantages and limitations of the approach were discussed.     

Accurate prediction of thermodynamic properties of alkyl peroxides by combining density functional theory calculation with least-square calibration

Owing to the significance in kinetic modeling of the oxidation and combustion mechanisms of hydrocarbons, a fast and relatively accurate method was developed for the prediction of Delta(f)H(298)(o) of alkyl peroxides. By this method, a raw Delta(f)H(298)(o) value was calculated from the optimized geometry and vibration frequencies at B3LYP/6-31G(d,p) level and then an accurate Delta(f)H(298)(o) value was obtained by a least-square procedure. The least-square procedure is a six-parameter linear equation and is validated by a leave-one out technique, giving a cross-validation squared correlation coefficient q(2) of 0.97 and a squared correlation coefficient of 0.98 for the final model. Calculated results demonstrated that the least-square calibration leads to a remarkable reduction of error and to the accurate Delta(f)H(298)(o) values within the chemical accuracy of 8 kJ mol(-1) except (CH(3))(2)CHCH(2)CH(2)CH(2)OOH which has an error of 8.69 kJ mol(-1). Comparison of the results by CBS-Q, CBS-QB3, G2, and G3 revealed that B3LYP/6-31G(d,p) in combination with a least-square calibration is reliable in the accurate prediction of the standard enthalpies of formation for alkyl peroxides. Standard entropies at 298 K and heat capacities in the temperature range of 300-1500 K for alkyl peroxides were also calculated using the rigid rotor-harmonic oscillator approximation.