Simultaneous parametrization of torsional and third-neighbor interaction terms in force-field development: The LLS-SC algorithm.

The calibration of torsional interaction terms by fitting relative gas-phase conformational energies against their quantum-mechanical values is a common procedure in force-field development. However, much less attention has been paid to the optimization of third-neighbor nonbonded interaction parameters, despite their strong coupling with the torsions. This article introduces an algorithm termed LLS-SC, aimed at simultaneously parametrizing torsional and third-neighbor interaction terms based on relative conformational energies. It relies on a self-consistent (SC) procedure where each iteration involves a linear least-squares (LLS) regression followed by a geometry optimization of the reference structures. As a proof-of-principle, this method is applied to obtain torsional and third-neighbor interaction parameters for aliphatic chains in the context of the GROMOS 53A6 united-atom force field. The optimized parameter set is compared to the original one, which has been fitted manually against thermodynamic properties for small linear alkanes. The LLS-SC implementation is freely available under http://github.com/mssm-labmmol/profiler.

Evaluation of nine condensed-phase force fields of the GROMOS, CHARMM, OPLS, AMBER, and OpenFF families against experimental cross-solvation free energies.

Experimental solvation free energies are nowadays commonly included as target properties in the validation of condensed-phase force fields, sometimes even in their calibration. In a previous article [Kashefolgheta et al., J. Chem. Theory. Comput., 2020, 16, 7556-7580], we showed how the involved comparison between experimental and simulation results could be made more systematic by considering a full matrix of cross-solvation free energies . For a set of N molecules that are all in the liquid state under ambient conditions, such a matrix encompasses N×N entries for considering each of the N molecules either as solute (A) or as solvent (B). In the quoted study, a cross-solvation matrix of 25 × 25 experimental value was introduced, considering 25 small molecules representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This experimental data was used to compare the relative accuracies of four popular condensed-phase force fields, namely GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In the present work, the comparison is extended to five additional force fields, namely GROMOS-54A7, GROMOS-ATB, OPLS-LBCC, AMBER-GAFF2, and OpenFF. Considering these nine force fields, the correlation coefficients between experimental values and simulation results range from 0.76 to 0.88, the root-mean-square errors (RMSEs) from 2.9 to 4.8 kJ mol-1, and average errors (AVEEs) from -1.5 to +1.0 kJ mol-1. In terms of RMSEs, GROMOS-2016H66 and OPLS-AA present the best accuracy (2.9 kJ mol-1), followed by OPLS-LBCC, AMBER-GAFF2, AMBER-GAFF, and OpenFF (3.3 to 3.6 kJ mol-1), and then by GROMOS-54A7, CHARM-CGenFF, and GROMOS-ATB (4.0 to 4.8 kJ mol-1). These differences are statistically significant but not very pronounced, and are distributed rather heterogeneously over the set of compounds within the different force fields.

Evaluating Classical Force Fields against Experimental Cross-Solvation Free Energies.

Experimental solvation free energies are nowadays commonly included as target properties in the validation and sometimes even in the calibration of condensed-phase force fields. However, this is often done in a nonsystematic fashion, by considering available solvation free energies involving an arbitrary collection of solutes in a limited set of solvents (e.g., water, octanol, chloroform, cyclohexane, or hexane). Here, this approach is made more systematic by introducing the concept of cross-solvation free energies ΔsGA:B⊖ for a set of N molecules that are all in the liquid state under ambient conditions, namely the matrix of N2 entries for ΔsGA:B⊖ considering each of the N molecules either as a solute (A) or as a solvent (B). Relying on available experimental literature followed by careful data curation, a complete ΔsGA:B⊖ matrix of 625 entries is constructed for 25 molecules with one to seven carbon atoms representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This matrix is then used to compare the relative accuracies of four popular condensed-phase force fields: GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In broad terms, and in spite of very different force-field functional-form choices and parametrization strategies, the four force fields are found to perform similarly well. Relative to the experimental values, the root-mean-square errors range between 2.9 and 4.0 kJ·mol-1 (lowest value of 2.9 for GROMOS and OPLS), and the average errors range between -0.8 and +1.0 kJ·mol-1 (lowest magnitude of 0.2 for AMBER and CHARMM). These differences are statistically significant but not very pronounced, especially considering the influence of outliers, some of which possibly caused by inaccurate experimental data.

Correction: Developing force fields when experimental data is sparse: AMBER/GAFF-compatible parameters for inorganic and alkyl oxoanions.

Correction for 'Developing force fields when experimental data is sparse: AMBER/GAFF-compatible parameters for inorganic and alkyl oxoanions' by Sadra Kashefolgheta et al., Phys. Chem. Chem. Phys., 2017, 19, 20593-20607.

Developing force fields when experimental data is sparse: AMBER/GAFF-compatible parameters for inorganic and alkyl oxoanions.

We present a set of Lennard-Jones parameters for classical, all-atom models of acetate and various alkylated and non-alkylated forms of sulfate, sulfonate and phosphate ions, optimized to reproduce their interactions with water and with the physiologically relevant sodium, ammonium and methylammonium cations. The parameters are internally consistent and are fully compatible with the Generalized Amber Force Field (GAFF), the AMBER force field for proteins, the accompanying TIP3P water model and the sodium model of Joung and Cheatham. The parameters were developed primarily relying on experimental information - hydration free energies and solution activity derivatives at 0.5 m concentration - with ab initio, gas phase calculations being used for the cases where experimental information is missing. The ab initio parameterization scheme presented here is distinct from other approaches because it explicitly connects gas phase binding energies to intermolecular interactions in solution. We demonstrate that the original GAFF/AMBER parameters often overestimate anion-cation interactions, leading to an excessive number of contact ion pairs in solutions of carboxylate ions, and to aggregation in solutions of divalent ions. GAFF/AMBER parameters lead to excessive numbers of salt bridges in proteins and of contact ion pairs between sodium and acidic protein groups, issues that are resolved by using the optimized parameters presented here.

Is the Solution Activity Derivative Sufficient to Parametrize Ion-Ion Interactions? Ions for TIP5P Water.

Biomolecular processes involve hydrated ions, and thus molecular simulations of such processes require accurate force-field parameters for these ions. In the best force-fields, both ion-water and anion-cation interactions are explicitly parametrized. First, the ion Lennard-Jones parameters are optimized to reproduce, for example, single ion solvation free energies; then ion-pair interactions are adjusted to match experimental activity or activity derivatives. Here, we apply this approach to derive optimized parameters for concentrated NaCl, KCl, MgCl2, and CaCl2 salt solutions, to be used with the TIP5P water model. These parameters are of interest because of a number of desirable properties of the TIP5P water model, especially for the simulation of carbohydrates. The results show, that this state of the art approach is insufficient, because the activity derivative often reaches a plateau near the target experimental value, for a wide range of parameter values. The plateau emerges from the interconversion between different types of ion pairs, so parameters leading to equally good agreement with the target solution activity or activity derivative yield very different solution structures. To resolve this indetermination, a second target property, such as the experimentally determined ion-ion coordination number, is required to uniquely determine anion-cation interactions. Simulations show that combining activity derivatives and coordination number as experimental target properties to parametrize ion-ion interactions, is a powerful method for reliable ion-water force field parametrization, and gives insight into the concentration of contact or solvent shared ion pairs in a wide range of salt concentrations. For the alkali and halide ions Li+, Rb+, Cs+, F-, Br-, and I-, we present ion-water parameters appropriate at infinite dilution only.

Computational Replication of the Primary Isotope Dependence of Secondary Kinetic Isotope Effects in Solution Hydride-Transfer Reactions: Supporting the Isotopically Different Tunneling Ready State Conformations.

We recently reported a study of the steric effect on the 1° isotope dependence of 2° KIEs for several hydride-transfer reactions in solution (J. Am. Chem. Soc. 2015, 137, 6653). The unusual 2° KIEs decrease as the 1° isotope changes from H to D, and more in the sterically hindered systems. These were explained in terms of a more crowded tunneling ready state (TRS) conformation in D-tunneling, which has a shorter donor-acceptor distance (DAD) than in H-tunneling. To examine the isotopic DAD difference explanation, in this paper, following an activated motion-assisted H-tunneling model that requires a shorter DAD in a heavier isotope transfer process, we computed the 2° KIEs at various H/D positions at different DADs (2.9 Å to 3.5 Å) for the hydride-transfer reactions from 2-propanol to the xanthylium and thioxanthylium ions (Xn(+) and TXn(+)) and their 9-phenyl substituted derivatives (Ph(T)Xn(+)). The calculated 2° KIEs match the experiments and the calculated DAD effect on the 2° KIEs fits the observed 1° isotope effect on the 2° KIEs. These support the motion-assisted H-tunneling model and the isotopically different TRS conformations. Furthermore, it was found that the TRS of the sterically hindered Ph(T)Xn(+) system does not possess a longer DAD than that of the (T)Xn(+) system. This predicts a no larger 1° KIE in the former system than in the latter. The observed 1° KIE order is, however, contrary to the prediction. This implicates the stronger DAD-compression vibrations coupled to the bulky Ph(T)Xn(+) reaction coordinate.

Ionic Liquids as Electrolytes for Electrochemical Double-Layer Capacitors: Structures that Optimize Specific Energy.

Key parameters that influence the specific energy of electrochemical double-layer capacitors (EDLCs) are the double-layer capacitance and the operating potential of the cell. The operating potential of the cell is generally limited by the electrochemical window of the electrolyte solution, that is, the range of applied voltages within which the electrolyte or solvent is not reduced or oxidized. Ionic liquids are of interest as electrolytes for EDLCs because they offer relatively wide potential windows. Here, we provide a systematic study of the influence of the physical properties of ionic liquid electrolytes on the electrochemical stability and electrochemical performance (double-layer capacitance, specific energy) of EDLCs that employ a mesoporous carbon model electrode with uniform, highly interconnected mesopores (3DOm carbon). Several ionic liquids with structurally diverse anions (tetrafluoroborate, trifluoromethanesulfonate, trifluoromethanesulfonimide) and cations (imidazolium, ammonium, pyridinium, piperidinium, and pyrrolidinium) were investigated. We show that the cation size has a significant effect on the electrolyte viscosity and conductivity, as well as the capacitance of EDLCs. Imidazolium- and pyridinium-based ionic liquids provide the highest cell capacitance, and ammonium-based ionic liquids offer potential windows much larger than imidazolium and pyridinium ionic liquids. Increasing the chain length of the alkyl substituents in 1-alkyl-3-methylimidazolium trifluoromethanesulfonimide does not widen the potential window of the ionic liquid. We identified the ionic liquids that maximize the specific energies of EDLCs through the combined effects of their potential windows and the double-layer capacitance. The highest specific energies are obtained with ionic liquid electrolytes that possess moderate electrochemical stability, small ionic volumes, low viscosity, and hence high conductivity, the best performing ionic liquid tested being 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide.

Computational replication of the abnormal secondary kinetic isotope effects in a hydride transfer reaction in solution with a motion assisted H-tunneling model.

We recently reported abnormal secondary deuterium kinetic isotope effects (2° KIEs) for hydride transfer reactions from alcohols to carbocations in acetonitrile (Chem. Comm. 2012, 48, 11337). Experimental 2° KIE values were found to be inflated on the 9-C position in the xanthylium cation but deflated on the ß-C position in 2-propanol with respect to the values predicted by the semi-classical transition-state theory. No primary (1°) isotope effect on 2° KIEs was observed. Herein, the KIEs were replicated by the Marcus-like H-tunneling model that requires a longer donor-acceptor distance (DAD) in a lighter isotope transfer process. The 2° KIEs for a range of potential tunneling-ready-states (TRSs) of different DADs were calculated and fitted to the experiments to find the TRS structure. The observed no effect of 1° isotope on 2° KIEs is explained in terms of the less sterically hindered TRS structure so that the change in DAD due to the change in 1° isotope does not significantly affect the reorganization of the 2° isotope and hence the 2° KIE. The effect of 1° isotope on 2° KIEs may be expected to be more pronounced and thus observable in reactions occurring in restrictive environments such as the crowded and relatively rigid active site of enzymes.

Imbalanced tunneling ready states in alcohol dehydrogenase model reactions: rehybridization lags behind H-tunneling.

The secondary kinetic isotope effects for the hydride transfer reactions from aliphatic alcohols to two carbocations (NAD(+) models) in acetonitrile were determined. The results suggest that the hydride transfer takes place by tunneling and that the rehybridizations of both donor and acceptor carbons lag behind the H-tunneling. This is quite contrary to the observations in alcohol dehydrogenases where the importance of enzyme motions in catalysis is manifested.