Chiral force constants: Recommendations for the presentation of internal coordinate force constants

Some force constants associated with the internal coordinates that sense handedness or chirality can have opposite signs in the enantiomers of chiral molecules. Examples of such force constants include interaction force constants between a torsional and stretching or bending internal coordinates. The sign reversal for these force constants in the enantiomers of chiral molecules or in opposite-handed molecular segments is best recognized by labeling them as chiral force constants. Recognition of chiral force constants suggests that certain guidelines are to be followed in the presentation of internal coordinate force constants. © 1993 by John Wiley & Sons, Inc. © John Wiley & Sons, Inc.     

A molecular mechanics force field for lignin

A CHARMM molecular mechanics force field for lignin is derived. Parameterization is based on reproducing quantum mechanical data of model compounds. Partial atomic charges are derived using the RESP electrostatic potential fitting method supplemented by the examination of methoxybenzene:water interactions. Dihedral parameters are optimized by fitting to critical rotational potentials and bonded parameters are obtained by optimizing vibrational frequencies and normal modes. Finally, the force field is validated by performing a molecular dynamics simulation of a crystal of a lignin fragment molecule and comparing simulation-derived structural features with experimental results. Together with the existing force field for polysaccharides, this lignin force field will enable full simulations of lignocellulose.     

End-to-end differentiable construction of molecular mechanics force fields

Molecular mechanics (MM) potentials have long been a workhorse of computational chemistry. Leveraging accuracy and speed, these functional forms find use in a wide variety of applications in biomolecular modeling and drug discovery, from rapid virtual screening to detailed free energy calculations. Traditionally, MM potentials have relied on human-curated, inflexible, and poorly extensible discrete chemical perception rules (atom types) for applying parameters to small molecules or biopolymers, making it difficult to optimize both types and parameters to fit quantum chemical or physical property data. Here, we propose an alternative approach that uses graph neural networks to perceive chemical environments, producing continuous atom embeddings from which valence and nonbonded parameters can be predicted using invariance-preserving layers. Since all stages are built from smooth neural functions, the entire process—spanning chemical perception to parameter assignment—is modular and end-to-end differentiable with respect to model parameters, allowing new force fields to be easily constructed, extended, and applied to arbitrary molecules. We show that this approach is not only sufficiently expressive to reproduce legacy atom types, but that it can learn to accurately reproduce and extend existing molecular mechanics force fields. Trained with arbitrary loss functions, it can construct entirely new force fields self-consistently applicable to both biopolymers and small molecules directly from quantum chemical calculations, with superior fidelity than traditional atom or parameter typing schemes. When adapted to simultaneously fit partial charge models, espaloma delivers high-quality partial atomic charges orders of magnitude faster than current best-practices with low inaccuracy. When trained on the same quantum chemical small molecule dataset used to parameterize the Open Force Field (“Parsley”) openff-1.2.0 small molecule force field augmented with a peptide dataset, the resulting espaloma model shows superior accuracy vis-á-vis experiments in computing relative alchemical free energy calculations for a popular benchmark. This approach is implemented in the free and open source package espaloma, available at https://github.com/choderalab/espaloma.

Graph neural network-based continuous embedding is used to replace a human expert-derived discrete atom typing scheme to parametrize accurate and extensible molecular mechanics force fields.

Molecular mechanics (MM) force fields—physical models that abstract molecular systems as atomic point masses that interact via nonbonded interactions and valence (bond, angle, torsion) terms—have powered in silico modeling to provide key insights and quantitative predictions in all aspects of chemistry, from drug discovery to materials science.^1–9 While recent work in quantum machine learning (QML) potentials has demonstrated how flexibility in functional forms and training strategies can lead to increased accuracy,^10–16 these QML potentials are orders of magnitude slower than popular molecular mechanics potentials even on expensive hardware accelerators, as they involve orders of magnitude more floating point operations per energy or force evaluation.On the other hand, the simpler physical energy functions of MM models are compatible with highly optimized implementations that can exploit a wide variety of hardware,^2,17–21 but rely on complex and inextensible legacy atom typing schemes for parameter assignment:²²• First, a set of rules is used to classify atoms into discrete atom types that must encode all information about an atom''s chemical environment needed in subsequent parameter assignment steps.• Next, a discrete set of bond, angle, and torsion types is determined by composing the atom types involved in the interaction.• Finally, the parameters attached to atoms, bonds, angles, and torsions are assigned according to a look-up table of composed parameter types.Consequently, atoms, bonds, angles, or torsions with distinct chemical environments that happen to fall into the same expert-derived discrete type class are forced to share the same set of MM parameters, potentially leading to low resolution and poor accuracy. Furthermore, the explosion of the number of discrete parameter classes describing equivalent chemical environments required by traditional MM atom typing schemes not only poses significant challenges to extending the space of atom types,²² but optimizing these independently has the potential to compromise generalizability and lead to overfitting. Even with modern MM parameter optimization frameworks^23–25 and sufficient data, parameter optimization is only feasible in the continuous parameter space defined by these fixed atom types, while the mixed discrete-continuous optimization problem—jointly optimizing types and parameters—is intractable.Here, we present a continuous alternative to traditional discrete atom typing schemes that permits full end-to-end differentiable optimization of both typing and parameter assignment stages, allowing an entire force field to be built, extended, and applied using standard automatic differentiation machine learning frameworks such as PyTorch,³² TensorFlow,³³ or JAX³⁴ (Fig. 1). Graph neural networks have recently emerged as a powerful way to learn chemical properties of atoms, bonds, and molecules for biomolecular species (both small organic molecules and biopolymers), which can be considered isomorphic with their graph representations.^35–44 We hypothesize that graph neural networks operating on molecules have expressiveness that is at least equivalent to—and likely much greater than—expert-derived typing rules, with the advantage of being able to smoothly interpolate between representations of chemical environments (such as accounting for fractional bond orders⁴⁵). We provide empirical evidence for this in Section 1.1.Open in a separate windowFig. 1Espaloma is an end-to-end differentiable molecular mechanics parameter assignment scheme for arbitrary organic molecules. Espaloma (extendable surrogate potential optimized by message-passing) is a modular approach for directly computing molecular mechanics force field parameters Φ_FF from a chemical graph such as a small molecule or biopolymer via a process that is fully differentiable in the model parameters Φ_NN. In Stage 1, a graph neural network is used to generate continuous latent atom embeddings describing local chemical environments from the chemical graph (Section 1.1). In Stage 2, these atom embeddings are transformed into feature vectors that preserve appropriate symmetries for atom, bond, angle, and proper/improper torsion inference via Janossy pooling (Section 1.2). In Stage 3, molecular mechanics parameters are directly predicted from these feature vectors using feed-forward neural nets (Section 1.3). This parameter assignment process is performed once per molecular species, allowing the potential energy to be rapidly computed using standard molecular mechanics or molecular dynamics frameworks thereafter. The collection of parameters Φ_NN describing the espaloma model can be considered as the equivalent complete specification of a traditional molecular mechanics force field such as GAFF^26,27/AM1-BCC^28,29 in that it encodes the equivalent of traditional typing rules, parameter assignment tables, and even partial charge models. This final stage is modular, and can be easily extended to incorporate additional molecular mechanics parameter classes, such as parameters for a charge-equilibration model (Section 4), point polarizabilities, or valence-coupling terms for class II molecular mechanics force fields.^30,31Next, we demonstrate the utility of such a model (which we call the extensible surrogate potential optimized by message-passing, or espaloma) to construct end-to-end optimizable force fields with continuous atom types. Espaloma assigns molecular mechanics parameters from a molecular graph in three differentiable stages (Fig. 1):• Stage 1: continuous atom embeddings are constructed using graph neural networks to perceive chemical environments (Section 1.1).• Stage 2: continuous bond, angle, and torsion embeddings are constructed using pooling that preserves appropriate symmetries (Section 1.2).• Stage 3: molecular mechanics force field parameters are computed from atom, bond, angle, and torsion embeddings using feed-forward networks (Section 1.3).Additional molecular mechanics parameter classes (such as point polarizabilities, valence coupling terms, or even parameters for charge-transfer models⁴⁶) can easily be added in a modular manner.Similar to legacy molecular mechanics parameter assignment infrastructures, molecular mechanics parameters are assigned once for each system, and can be subsequently used to compute energies and forces or carry out molecular simulations with standard molecular mechanics packages. Unlike traditional legacy force fields, espaloma model parameters Φ_NN—which define the entire process by which molecular mechanics force field parameters Φ_FF are generated ad hoc for a given molecule—can easily be fit to data from scratch using standard, highly portable, high-performance machine learning frameworks that support automatic differentiation.Here, we demonstrate that espaloma provides a sufficiently flexible representation to both learn to apply existing MM force fields and to generalize them to new molecules (Section 2). Espaloma''s modular loss function enables force fields to be learned directly from quantum chemical energies (Section 3), partial charges (Section 4), or both. The resulting force fields can generate self-consistent parameters for small molecules, biopolymers (Section 5), and covalent adducts (Section 1). Finally, an espaloma model fit to the same quantum chemical dataset used to build the Open Force Field OpenFF (“Parsley”) openff-1.2.0 small molecule force field, augmented with peptide quantum chemical data, can outperform it in an out-of-sample kinase : inhibitor alchemical free energy benchmark (Section A.4 in ESI^†).     

Determination of the molecular force constants of colloid particles

The possibility of analysis of the shape of colloid particles in solutions of latex and French and Moldovan wines with a flow-type ultramicroscope and digital photo camera was studied. The micrographs taken can be used to determine distances between particles by comparison with particle radii found by an independent method. The concepts of Deryagin-Landau-Vervey-Overbeek’s physical theory of stability of colloid systems were used to calculate the molecular force constants for wine particles.     

New out-of-plane angle and bond angle internal coordinates and related potential energy functions for molecular mechanics and dynamics simulations

With currently used definitions of out-of-plane angle and bond angle internal coordinates, Cartesian derivatives have singularities, at ±π/2 in the former case and π in the latter. If either of these occur during molecular mechanics or dynamics simulations, the forces are not well defined. To avoid such difficulties, we provide new out-of-plane and bond angle coordinates and associated potential energy functions that inherently avoid these singularities. The application of these coordinates is illustrated by ab initio calculations on ammonia, water, and carbon dioxide. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1067–1084, 1999     

The use of internal cartesian coordinates for describing molecular vibrations

An analysis of the influence of isotope substitution on the system of electronic-nuclear equations for an arbitrary molecular system was used as a basis for formulating invariance conditions with respect to isotope substitution of the potential energy surface written in the Cartesian coordinates rigidly bound with the center of mass of the molecule (internal Cartesian coordinates). This property of the potential function obviates the necessity of using curvilinear natural coordinates, which can be replaced by Cartesian coordinates, in theoretical studies of the vibrational spectra of molecules and their isotopomers and in solving the direct and inverse anharmonic problems. An equation for the quantum-mechanical Hamiltonian of a normal molecule in internal Cartesian coordinates was obtained.     

Use of the virial theorem to calculate molecular force constants. Harmonic force constants of the water molecule

The quantum mechanical virial theorem is used to calculate the harmonic part of the potential function for the water molecule on the basis of wave functions obtained in the self-consistent field approximation.     

Strained ring carbonyl force constants from molecular orbital calculations

Molecular orbital theory is used to evaluate the carbonyl force constants in several compounds. The carbonyl force constants k_CO are found to be essentially unchanged in the series: acetone, cyclopentanone, and cyclobutanone; however, the value of k_CO is calculated to be 5% higher in the case of cyclopropanone.     

Theoretical study of dynamical properties on reaction path in molecular internal coordinates

The dynamical properties on reaction path (IRC) in internal coordinates have been obtained, which in. clude ω_K (frequencies orthogonal to IRC), L_K (vibrational modes), B_(KF) (coupling constants between the IRC and vibra tions orthogonal to it), B_(KL) (coupling constants between every two vibrations orthogonal to IRC). A set of theory of teac. tion path in molecular internal coordinates has been also constructed. The dynamical properties, including ω_K, B_(KF) B_(KL) of the reaction H~1O~2H~3 H~4→H~1O~2 H~3H~4 have been calculated, which explicitly explain the interaction, chang ing trend and contribution of each chemical bond (including bond angle) in the reaction.     

Derivatives of molecular Vibrational force constants with a singular Jacobian

Counter examples are presented that negate the generality of Gans' claim that a singular jacobian demands the existence of one off-diagonal force constant F_st possessing the extremal property ?F_rr/?F_st = ∞ for all of the n diagonal force constants F_rr (r = 1, > 62; n) in a symmetry species containing n vibrations.     

Geometry optimization of molecular clusters and complexes using scaled internal coordinates

Scaled internal coordinates are introduced for use in the geometry optimization of systems composed of multiple fragments, such as solvated molecules, clusters, and biomolecular complexes. The new coordinates are related to bond lengths, bond angles and torsion angles by geometry-dependent scaling factors. The scaling factors serve to expedite the optimization of complexes containing outlying fragments, without hindering the optimization of the intramolecular degrees of freedom. Trial calculations indicate that, at asymptotic separations, the scaling factors improve the rate of convergence by a factor of 4 to 5.     

A molecular mechanics force field for the cobalt corrinoids

A force field for the cobalt (III) corrinoids (derivatives of vitamin B₁₂) for use with a modified version of the molecular mechanics program 2(87) has been developed empirically around 19 cobalt corrinoid crystal structures. Bond lengths, bond angles and torsional angles are reproduced with r.m.s. differences of 0.01 Å, 2.4 °, and 4.2 °, respectively, within the standard deviation of the mean of these parameters found in the solid state. The axial ligand occupying the lower coordination site in the cobalamins, 5,6-dimethylbenzimidazole, is shown to have very limited rotational freedom and is constrained by the downward-pointing b and d propionamide side chains of the corrin ring. Strain-energy profiles for rotation of the side chains of the corrin ring show the existence of several local energy minima and this explains the observed variability in the orientations of these side chains in the solid state. The known change in conformation which occurs in the C ring when the e side chain is epimerized from the lower to the upper face of the corrin ring in cyano-13-epicobalamin is correctly predicted, provided the starting conformation of the C ring is unbiased. A study of cyano-8-epicobalamin indicates that an analogous conformational change does not occur in the B ring and the epimerized d side chain assumes an equatorial orientation relative to the corrin ring. Parameters for the Co---C bond in alkylcobalamins were developed and the structure of methyl- and adenosylcobalamin are accurately reproduced. An examination of the strain energy consequences of rotation of the adenosyl ligand about the Co---C bond identifies a number of low-energy conformations at least two of which, in which adenosyl lies over the “southern” and “eastern” portions of the corrin ring, respectively, have been previously deduced from NMR observations. Coordinated neopentyl in neopentylcobalamin is much more hindered to rotation about the Co---C bond and the lowest conformation finds two γ(C) atoms straddling the upwardly projecting C46 methyl group of the corrin.     

Thermodynamic molecular mechanics force field: Modified QCFF program

The QCFF program originated by Warshel and Karplus^4a was modified to compute accurate thermodynamic properties S^o, C, (H – H)/T, and ΔH for various acyclic and cyclic alkenes and alkadienes. Modifications consisted of adjusted bond angle, dihedral angle, bond stretch, and bond energy parameters that improved calculated vibrational frequencies, zero point energies, and thermodynamic functions. Supplemental torsional potential energy functions that were added to existing torsional functions led to greatly improved relative conformer energies and ΔH values. It was shown that inclusion of hindered internal rotation leads to significantly better agreement of calculated thermodynamic functions with observed values for acyclic alkenes at high temperatures. The calculated thermodynamic properties of the alkenes and alkadienes were deemed sufficiently accurate for calculation of standard enthalpies and Gibbs free energies of gas phase chemical reactions at various temperatures. © 1994 by John Wiley & Sons, Inc.     

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.     

Partial hessian fitting for determining force constant parameters in molecular mechanics

We present a new protocol for deriving force constant parameters that are used in molecular mechanics (MM) force fields to describe the bond‐stretching, angle‐bending, and dihedral terms. A 3 × 3 partial matrix is chosen from the MM Hessian matrix in Cartesian coordinates according to a simple rule and made as close as possible to the corresponding partial Hessian matrix computed using quantum mechanics (QM). This partial Hessian fitting (PHF) is done analytically and thus rapidly in a least‐squares sense, yielding force constant parameters as the output. We herein apply this approach to derive force constant parameters for the AMBER‐type energy expression. Test calculations on several different molecules show good performance of the PHF parameter sets in terms of how well they can reproduce QM‐calculated frequencies. When soft bonds are involved in the target molecule as in the case of secondary building units of metal‐organic frameworks, the MM‐optimized geometry sometimes deviates significantly from the QM‐optimized one. We show that this problem is rectified effectively by use of a simple procedure called Katachi that modifies the equilibrium bond distances and angles in bond‐stretching and angle‐bending terms. © 2016 Wiley Periodicals, Inc.     

Liquid-phase force field and dipole moment derivatives with respect to internal coordinates of benzene

This paper presents an analysis of the infrared vibrational intensities found for C(6)H(6), C(6)D(6) and C(6)H(5)D in the liquid phase, motivated in part by the quite marked intensity differences between the fundamentals of C(6)H(6) and C(6)D(6) in the liquid, and between corresponding vibrations in the liquid and gas phases. The analysis is carried out under the harmonic approximation and results from a determination of the force field for liquid C(6)H(6), C(6)D(6) and C(6)H(5)D. The force constants for the liquid-phase are presented and compared to those in the literature for the gas-phase. Previously reported experimental intensities are used along with the eigenvectors of the force field analysis to determine the dipole moment derivatives with respect to symmetry and internal coordinates. The dipole moment derivatives with respect to internal coordinates obtained are partial differentialmicro/ partial differentials=0.38+/-0.02DebyeA(-1), partial differentialmicro/ partial differentialt=0.24+/-0.01, partial differentialmicro/ partial differentialbeta=0.26+/-0.01, and partial differentialmicro/ partial differentialgamma=0.64+/-0.03DebyeA(-1). There is very little difference between the dipole moment derivatives with respect to internal coordinates obtained from non-linear least squares fitting of the two D(6h) isotopomers and those obtained from non-linear least squares fitting of the three isotopomers. The results show that there is significant intensity sharing in the CH stretch region of C(6)H(5)D between the fundamental and combination bands.     

MATCH: an atom-typing toolset for molecular mechanics force fields

We introduce a toolset of program libraries collectively titled multipurpose atom-typer for CHARMM (MATCH) for the automated assignment of atom types and force field parameters for molecular mechanics simulation of organic molecules. The toolset includes utilities for the conversion of multiple chemical structure file formats into a molecular graph. A general chemical pattern-matching engine using this graph has been implemented whereby assignment of molecular mechanics atom types, charges, and force field parameters are achieved by comparison against a customizable list of chemical fragments. While initially designed to complement the CHARMM simulation package and force fields by generating the necessary input topology and atom-type data files, MATCH can be expanded to any force field and program, and has core functionality that makes it extendable to other applications such as fragment-based property prediction. In this work, we demonstrate the accurate construction of atomic parameters of molecules within each force field included in CHARMM36 through exhaustive cross validation studies illustrating that bond charge increment rules derived from one force field can be transferred to another. In addition, using leave-one-out substitution it is shown that it is also possible to substitute missing intra and intermolecular parameters with ones included in a force field to complete the parameterization of novel molecules. Finally, to demonstrate the robustness of MATCH and the coverage of chemical space offered by the recent CHARMM general force field (Vanommeslaeghe, et al., J Comput Chem 2010, 31, 671), one million molecules from the PubChem database of small molecules are typed, parameterized, and minimized.     

Molecular mechanics calculations for conformational energies and torional force constants in fluoro propenes and butenes

Experimental data on conformational energies of the molecules FH₂CHCCH₂, FH₂CFCCH₂, FH₂C(CH₃)C&.dbnd;CH₂ trans-FH₂CHCC(CH₃)H have been used to establish parameter values for the nonbonding atom ⋯ atom interaction F ⋯ C(sp²) within the Morse potential formulation. Torsional potentials have been calculated for the four molecules mentioned above and in addition for cis- and trans-FH₂CHCCHF, (FH₂C)₂CCH₂, cis-FH₂CHCCHCH₂F, CH₃FCHHCCH₂ and FH₂CCH₂HCCH₂. Calculated results have been compared with experimental values. Torsional force constants for the molecules have been obtained. A comparison between fluoro, chloro and bromo compounds is presented.