PaperlayerView exact revision on arXiv

Source-linked AI summary

Towards Exact Molecular Dynamics Simulations with Machine-Learned Force Fields

Stefan Chmiela, Huziel E. Sauceda, Klaus-Robert Müller, Alexandre Tkatchenko

arXiv:1802.09238v3physics.chem-ph

TL;DR

Classical interatomic potentials lack accuracy for predictive molecular dynamics, while direct CCSD(T) potential-energy-surface fitting has been practical only for small, rigid molecules. This paper incorporates spatial and temporal physical symmetries into sGDML in a data-driven way, enabling converged MD with fully quantized electrons and nuclei for flexible molecules with up to a few dozen atoms.

  • Problem

    Classical interatomic potentials lack accuracy for predictive MD, while direct CCSD(T) potential-energy-surface fitting has been practical only for small, rigid molecules.

    Source: §I · p. 1

  • Method

    The sGDML approach incorporates relevant molecular spatial and temporal symmetries through data-driven discovery and use of those symmetries.

    Source: §I · p. 1, A. Symmetrized gradient-domain machine… · p. 3

  • Results

    sGDML enables converged MD with fully quantized electrons and nuclei for flexible molecules with up to a few dozen atoms.

    Source: §I · p. 1

  • Takeaways & Limitations

    These simulations pave the way to computing molecular dynamical and thermodynamical properties with an essentially exact underlying potential-energy surface.

    Source: §III · p. 7

    The work focuses on intramolecular forces in small- and medium-sized molecules, and extending sGDML to larger systems remains a challenge.

    Source: §III · p. 9

Abstract

from arXiv · show

Molecular dynamics (MD) simulations employing classical force fields constitute the cornerstone of contemporary atomistic modeling in chemistry, biology, and materials science. However, the predictive power of these simulations is only as good as the underlying interatomic potential. Classical potentials often fail to faithfully capture key quantum effects in molecules and materials. Here we enable the direct construction of flexible molecular force fields from high-level ab initio calculations by incorporating spatial and temporal physical symmetries into a gradient-domain machine learning (sGDML) model in an automatic data-driven way. The developed sGDML approach faithfully reproduces global force fields at quantum-chemical CCSD(T) level of accuracy and allows converged molecular dynamics simulations with fully quantized electrons and nuclei. We present MD simulations, for flexible molecules with up to a few dozen atoms and provide insights into the dynamical behavior of these molecules. Our approach provides the key missing ingredient for achieving spectroscopic accuracy in molecular simulations.

I. INTRODUCION

The introduction frames inaccurate classical potentials and costly or data-hungry alternatives as barriers to predictive, high-accuracy molecular dynamics. It presents sGDML as a symmetry-based approach intended to enable CCSD(T)-level simulations for flexible molecules with up to a few dozen atoms.

I. INTRODUCION

  • A nanosecond CCSD(T) trajectory for one ethanol molecule would take roughly a million CPU years; direct CCSD(T) PES fitting has been practical only for small, rigid molecules.Common AIMD instead uses density-functional approximations, which can yield contrasting predictions.Source: §I · p. 1
  • Direct machine-learning fits typically require thousands or millions of configurations, whereas high-level ab initio reference data are practical for only hundreds of conformations.This data requirement limits the use of those models with high-level methods such as CCSD(T).Source: §I · p. 1
  • sGDML uses data-driven spatial and temporal symmetries to construct high-level force fields and enable converged simulations with quantized electrons and nuclei for molecules of up to a few dozen atoms.The authors report carrying out molecular dynamics at the ab initio coupled-cluster level and present dynamical insights, positioning the approach as a step toward spectroscopic accuracy.Source: §I · pp. 1, 3
  • A multipartite matching algorithm identifies globally consistent atom assignments across conformations, then symmetrizes training through a kernel that combines symmetric variants while ignoring dataset-irrelevant symmetries.The global matching is recovered by synchronizing approximate pairwise assignments and taking their transitive closure.Source: Figure 1 · p. 2

II. RESULTS · A. Symmetrized gradient-domain machine learning

sGDML extends GDML with spatial and molecular dynamic symmetries discovered from MD trajectories, enabling symmetric force-field models that retain the original model’s parameter count. Its matching procedure enforces globally consistent atom assignments so the model can use physically feasible molecular permutations.

A. Symmetrized gradient-domain machine learning

  • sGDML incorporates spatial and molecular symmetries into GDML [47], extending it to force fields that support MD simulations at high-level ab initio accuracy.The symmetries include rotational and translational invariance, reflection, and dynamic non-rigid motions such as methyl-group rotations.Source: A. Symmetrized gradient-domain machine… · pp. 3–4
  • A data-driven algorithm discovers relevant molecular symmetries from MD trajectories, including non-rigid symmetries arising from functional-group rotations and torsional displacements.Such dynamic symmetries are typically not incorporated in traditional force fields and electronic-structure calculations.Source: A. Symmetrized gradient-domain machine… · p. 3
  • Atom correspondence is first estimated by minimizing squared Euclidean distance between adjacency matrices of isomorphic molecular graphs, then refined through eigenvector overlap and Hungarian assignment [55, 56].Pairwise assignments use molecular graph adjacency matrices [9, 13, 47, 53, 54], but must be made globally consistent across the training set [50–52].Source: A. Symmetrized gradient-domain machine… · p. 3
  • A minimum spanning tree of the most confident pairwise assignments provides a transitive closure, yielding a rank-limited, globally consistent permutation matrix across training geometries.The tree minimizes the bipartite-matching cost over the training set.Source: A. Symmetrized gradient-domain machine… · p. 3
  • Symmetric training-set augmentation covers recovered permutational configurations, while coefficient symmetries contract the linear system to its original size without changing the full coefficients.The resulting model has the exact same number of parameters as the original, non-symmetric model.Source: A. Symmetrized gradient-domain machine… · pp. 3–4
  • The symmetric kernel extends GDML to force fields, and the sGDML force estimator combines derivative contributions from training points and symmetry transformations.Training geometries are Boltzmann-subsampled from a 200 picosecond DFT MD trajectory at 500 K; for sGDML@CCSD(T), energy and force labels are recomputed at the higher theory level.Source: A. Symmetrized gradient-domain machine… · p. 4

B. Forces and energies from GDML to sGDML@DFT to sGDML@CCSD(T)

sGDML improves force prediction and data efficiency by incorporating molecular symmetries, and the paper applies it to build high-accuracy models using CCSD(T) forces. The reported CCSD(T) models improve energy prediction over sGDML@DFT while maintaining comparable force accuracy, with benzene as the exception.

B. Forces and energies from GDML to sGDML@DFT to sGDML@CCSD(T)

  • With identical training sets, sGDML improved force prediction by 31.3% to 67.4% for seven molecules, while uracil and salicylic acid showed no improvement without exploitable symmetries.The benefit of symmetry was stronger for force predictions than energy predictions, most clearly for naphthalene, where force predictions improved unilaterally; the authors attribute this partly to omitting an energy penalty from the cost function.Source: B. Forces and energies from GDML to… · pp. 4–5
  • The study trained sGDML@CCSD(T) models for benzene, toluene, ethanol, and malonaldehyde, and used CCSD forces for aspirin.Its stated goal was to recover CCSD(T) force fields for flexible molecules with up to 20 atoms using only a few hundred molecular conformations.Source: B. Forces and energies from GDML to… · pp. 4–5
  • At MAE = 1 kcal mol^-1 Å^-1, sGDML required 200 training examples versus about 800 for GDML, except for aspirin; energy-based ML approaches typically require two to three orders of magnitude more data [47].The paper identifies this force accuracy as the minimum required for reliable molecular dynamics simulations.Source: B. Forces and energies from GDML to… · p. 5
  • For ethanol and toluene, sGDML@CCSD(T) reduced energy prediction error relative to sGDML@DFT by factors of 1.4 and 3.4, respectively; force accuracy was comparable except for benzene.The authors attribute the benzene exception to shifts in potential-energy-surface minima between DFT and CCSD(T), suggesting that CCSD(T) data sampling could be improved through corrected resampling.Source: B. Forces and energies from GDML to… · p. 5

C. Molecular dynamics with ab initio accuracy

Ethanol’s closely spaced conformational minima make global error scores and classical sampling inadequate; accurate reference energies and nuclear quantum effects are needed for predictive MD.

C. Molecular dynamics with ab initio accuracy

  • Experiment confirms Mt as the ground state and Mg as local minima; their 0.12 kcal mol−1 energy difference and 1.15 kcal mol−1 barrier make a 1 kcal mol−1 target inadequate.The Mg↔Mt rotational barrier is also described as approximately 1.2 kcal mol−1, making ethanol sensitive to nuclear quantum effects.Source: C. Molecular dynamics with ab initio… · p. 5
  • Global MAE and RMSE can overestimate force-field quality because they average out local topological properties relevant to MD trajectory analysis.Source: C. Molecular dynamics with ab initio… · p. 5
  • DFT(PBE-TS) reverses ethanol’s conformer ordering, whereas CCSD(T)/cc-pVTZ places Mt 0.08 kcal mol−1 below Mg, matching experiment.DFT(PBE-TS) instead makes Mg 0.08 kcal mol−1 more stable than Mt.Source: C. Molecular dynamics with ab initio… · p. 5
  • PIMD enables Mg↔Mt transitions and yields statistical weights in excellent agreement with experiment, unlike inadequate classical room-temperature sampling.The ethanol occupation comparison includes classical MD and PIMD for sGDML@CCSD(T) and sGDML@DFT against experiment.Source: C. Molecular dynamics with ab initio… · p. 5

Full MD with sGDML DFT & CCSD(T)

sGDML@CCSD(T) PIMD reveals strong temperature-dependent anharmonicity in ethanol that brings its hydroxyl torsional frequency closer to experiment. Simulations of malonaldehyde and aspirin show that electronic-structure differences can alter molecular-state sampling and vibrational spectra.

Full MD with sGDML DFT & CCSD(T)

  • sGDML@CCSD(T) predicts ethanol hydroxyl-rotation barriers of 1.18 kcal mol−1 for Mt →Mg, 1.19 kcal mol−1 for Mg- →Mg+, and 1.07 kcal mol−1 for Mg →Mt.The figure also compares classical and path-integral state occupations with experimental values.Source: Full MD with sGDML DFT & CCSD(T) · p. 6
  • CCSD(T) yields higher fingerprint-zone frequencies than DFT with similar spectral shapes but slightly different intensities; thermal shifts affect low-frequency, nonsymmetric modes nonuniformly.Each normal mode shifts specifically with temperature rather than by a common scaling factor.Source: Full MD with sGDML DFT & CCSD(T) · pp. 6–7
  • At room temperature, the hydroxyl-rotation mode drops to ∼215 cm−1 from 243.7 cm−1 in the CCSD(T) harmonic approximation, a 12% red-shift toward measured values of ∼202 and ∼207 cm−1.The shift is attributed to strong anharmonicity in the lowest normal mode, which mainly corresponds to hydroxyl-group rotations.Source: Full MD with sGDML DFT & CCSD(T) · p. 7
  • In malonaldehyde, sGDML@CCSD(T) samples a local-minimum structure nearly unvisited by sGDML@DFT despite similar potential-energy surfaces, while their spectra differ most near 500 cm−1.The molecule also fluctuates between symmetric global minima across a crossing barrier of ∼0.2 kcal mol−1.Source: Full MD with sGDML DFT & CCSD(T) · p. 7
  • For aspirin, sGDML@CCSD localizes PIMD sampling at the global minimum, whereas DFT delocalizes it; CCSD raises ester-dihedral barriers by ∼1 kcal mol−1.The DFT model also produces softer interatomic interactions and visibly different vibrational spectra; explicit CCSD(T) calculations corroborate the barrier difference.Source: Full MD with sGDML DFT & CCSD(T) · p. 7

III. DISCUSSION

sGDML enables accurate molecular dynamics for flexible molecules with up to a few dozen atoms by incorporating physical symmetries into a gradient-domain learning model. The discussion highlights its spectroscopic-accuracy goal, limitations of traditional force fields, and challenges in scaling to larger and condensed systems.

III. DISCUSSION

  • The sGDML model combines spatial, temporal, and local symmetries to simulate flexible molecules with up to a few dozen atoms at high-level ab initio accuracy.Local symmetries are system-size independent and persist when molecules are fragmented, unlike global point-group symmetries, whose availability diminishes as molecules grow.Source: §III · p. 7
  • Spectroscopic accuracy is typically about 0.03 kcal mol−1, and is needed because conformer energy differences are often 0.1–0.2 kcal mol−1.The discussion argues that molecular-PES accuracy should not be set by the 1 kcal mol−1 chemical or thermochemical accuracy convention, which was conceived for thermochemical measurements.Source: §III · p. 7
  • AMBER [61] misses key quantum effects: it lacks ethanol rotor coupling, samples bonded interactions harmonically, and yields a mismatched aspirin PES with spurious minima.For malonaldehyde, AMBER and sGDML distributions can look similar despite arising from different interactions: Coulomb forces versus electron-correlation effects.Source: §III · p. 8, Figure 5 · p. 9
  • Combining physical symmetries improves data efficiency by often more than two orders of magnitude, beyond the factor discussed in studies of invariance constraints.The authors speculate that these gains may reflect a less complex learning problem that is easier to represent.Source: §III · p. 8
  • Extending sGDML to larger and condensed systems requires transferable combinations of trained models, multilevel sampling, and integration with an accurate intermolecular force field [64].Other proposed directions include additional physical priors, PES dimensionality reduction, reaction pathways, and modeling infrared, Raman, and other spectra.Source: §III · p. 9

IV. METHODS · A. Reference data generation

Reference datasets were generated from a 500 K ab initio molecular-dynamics trajectory at the PBE+TS level, then reused geometries were recomputed with coupled-cluster methods for higher-level datasets.

A. Reference data generation

  • DFT training data came from a 200 ps NVT ab initio MD trajectory at 500 K, with forces and energies computed using all-electron PBE plus Tkatchenko-Scheffler van der Waals corrections.The trajectory used a 0.5 fs resolution and a Nosé-Hoover thermostat; the final data were subsampled while preserving the Maxwell-Boltzmann energy distribution.Source: A. Reference data generation · p. 9
  • Coupled-cluster datasets reused the DFT geometries and recomputed energies and forces with CCSD(T), using molecule-specific basis sets and Psi4.Ethanol used cc-pVTZ; toluene and malonaldehyde used cc-pVDZ; aspirin used CCSD/cc-pVDZ.Source: A. Reference data generation · p. 9

B. Molecular dynamics

The molecular dynamics simulations used path-integral molecular dynamics (PIMD) to incorporate quantum nuclear delocalization, with the sGDML model interfaced to i-PI. Ethanol was simulated for 1 ns using 16 beads to obtain converged sampling of the potential-energy surface.

  • PIMD was used to incorporate quantum-mechanical effects associated with nuclear delocalization into the molecular dynamics simulations.PIMD incorporates these effects via Feynman’s path integral formalism.Source: B. Molecular dynamics · pp. 9–10
  • The PIMD simulations used the sGDML model interfaced to the i-PI code.Source: B. Molecular dynamics · p. 10
  • Ethanol was simulated for 1 ns with 16 beads to obtain converged potential-energy-surface sampling, using a 0.2 fs timestep in NVE and NVT ensembles.The timestep was set to ensure energy conservation along the MD.Source: B. Molecular dynamics · p. 10

C. Bipartite matching cost matrix

The method matches eigenvectors of two molecular-graph adjacency matrices by solving an optimal assignment problem with Hungarian-algorithm costs. A penalty term discourages assignments between non-identical nuclei.

C. Bipartite matching cost matrix

D. Training sGDML

sGDML symmetrizes the GDML kernel by summing over relevant atom assignments for each training geometry, approximating a fully symmetrized training set while retaining the original kernel-matrix size.

D. Training sGDML

  • Summing over relevant atom assignments for each training geometry approximates similarities from a fully symmetrized training set while retaining the kernel matrix’s original size.This symmetric kernel formulation is used to symmetrize the GDML model.Source: D. Training sGDML · p. 10
  • The Hessian rows and columns in each summand are permuted so their corresponding partial derivatives align.Source: D. Training sGDML · p. 10
  • At model evaluation, the free kernel variable x is not permuted and the normalization factor is dropped; descriptor inputs are addressed in Supplementary Note 3.Source: D. Training sGDML · p. 10

Related papers

Loading 1802.09238v3…