Accurate global machine learning force fields for molecules with hundreds of atoms

Stefan Chmiela; Valentin Vassilev-Galindo; Oliver T. Unke; Adil Kabylda; Huziel E. Sauceda; Alexandre Tkatchenko; Klaus-Robert Müller

Source-linked AI summary

Accurate global machine learning force fields for molecules with hundreds of atoms

Stefan Chmiela, Valentin Vassilev-Galindo, Oliver T. Unke, Adil Kabylda, Huziel E. Sauceda, Alexandre Tkatchenko, Klaus-Robert Müller

arXiv:2209.14865v3Sep 29, 2022physics.chem-phphysics.comp-ph

TL;DR

Global ML force fields remain difficult to scale globally because large systems require costly all-atom coupling, while localization can omit long-range interactions. The paper develops an exact iterative sGDML training approach that retains full coupling and evaluates it on molecules of 42–370 atoms, including nanosecond-long simulations.

Problem
Global ML force fields are limited to a few dozen atoms, while localization assumptions for larger molecules can poorly describe non-local interactions.
Method
The paper combines closed-form and iterative training with numerical preconditioning to solve large global sGDML kernel models without approximating the original model.
Results
42–370 atoms are covered in the MD22 benchmark, and the approach supports stable nanosecond-long molecular dynamics simulations for supramolecular complexes.
Takeaways & Limitations
MD22 introduces benchmark challenges in molecular size, flexibility, and non-locality for advancing atomistic MLFF architectures.

Abstract

from arXiv · show

Global machine learning force fields (MLFFs), that have the capacity to capture collective many-atom interactions in molecular systems, currently only scale up to a few dozen atoms due a considerable growth of the model complexity with system size. For larger molecules, locality assumptions are typically introduced, with the consequence that non-local interactions are poorly or not at all described, even if those interactions are contained within the reference ab initio data. Here, we approach this challenge and develop an exact iterative parameter-free approach to train global symmetric gradient domain machine learning (sGDML) force fields for systems with up to several hundred atoms, without resorting to any localization of atomic interactions or other potentially uncontrolled approximations. This means that all atomic degrees of freedom remain fully correlated in the global sGDML FF, allowing the accurate description of complex molecules and materials that present phenomena with far-reaching characteristic correlation lengths. We assess the accuracy and efficiency of our MLFFs on a newly developed MD22 benchmark dataset containing molecules from 42 to 370 atoms. The robustness of our approach is demonstrated in nanosecond long path-integral molecular dynamics simulations for the supramolecular complexes in the MD22 dataset.

I. INTRODUCTION

Global ML force fields can capture interactions across all atoms, but their computational scaling limits global models to molecules of only a few dozen atoms. The paper addresses this by reducing effective degrees of freedom while retaining exact global coupling, enabling systems with several hundred atoms.

I. INTRODUCTION: Localization assumptions reduce degrees of freedom in large structures but can lose long-distance information and truncate long-range interactions.Local models often assume these interactions contribute little to overall dynamics, although long-range effects can matter.
I. INTRODUCTION: Global ML force fields include all interaction scales, but coupling at least a quadratic set of atom-atom interactions limits their computational scalability.This restriction has confined current global models to systems of only a few dozen atoms despite larger ab initio datasets being available.
I. INTRODUCTION: The proposed framework uses spectral analysis to show that large molecules have substantially fewer effective degrees of freedom than N^2.This supports a low-dimensional representation for scaling global MLFF kernel models to large molecules.
I. INTRODUCTION: A two-step procedure solves effective degrees of freedom in closed form, then iteratively converges remaining fluctuations to the exact full-problem solution.The approach lowers memory and computational time requirements simultaneously without localizing atomic interactions.
I. INTRODUCTION: The method targets Gaussian-process models and demonstrates global sGDML force fields for systems containing several hundred fully coupled atoms.The study includes supramolecular complexes, nanostructures, and biomolecular systems, and introduces MD22 with molecules ranging from 42 to 370 atoms.

II. LARGE-SCALE SGDML ALGORITHM

The large-scale sGDML algorithm combines closed-form and iterative GP training to reduce memory and computational demands while retaining globally coupled, exact solutions. It uses low-dimensional kernel structure, conjugate gradients, and Nyström-based preconditioning to scale models beyond previous system-size limits.

Model construction: sGDML incorporates conservation-law-derived constraints and invariances through a Gaussian-process kernel trained on force examples.The approach uses force data available from electronic-structure calculations and exploits GP closure under linear transformations.
Iterative solver: Conjugate-gradient descent selects conjugate search directions and optimal step sizes, making the iterative training scheme parameter-free.Fast-decaying eigenvalue spectra can allow early truncation of the CG expansion and super-linear convergence in practice.
Preconditioning: Increasing the number of inducing points strengthens preconditioning and can accelerate convergence, but increases construction time, memory requirements, and per-iteration cost.The spectrum of the preconditioned kernel is increasingly attenuated as inducing points increase, reducing its condition number.
Preconditioning: Preconditioning changes convergence control from the condition number of Kλ to that of P^-1Kλ while preserving the solver's symmetric positive-definite requirements.This addresses the slow or unstable convergence caused by ill-conditioned kernel matrices.
Preconditioning: The Nyström preconditioner approximates the relevant eigenspace using k columns selected with approximate statistical leverage scores instead of cubic-cost optimal rank-k decompositions.The approximation costs O(k^2m), while applying and inverting the resulting preconditioner have runtime and memory complexities O(mk^2) and O(mk).

III. ASSESSMENT OF LARGE-SCALE MOLECULAR FORCE FIELDS

The assessment examines preconditioned iterative training and large-molecule scalability, showing reliable convergence, reduced cost with suitable preconditioners, and force predictions across MD22 systems up to 370 atoms.

Convergence properties: Above 10% preconditioning, iterative training has lower memory and runtime complexity than the closed-form solver.Increasing preconditioner strength eventually reverses the runtime benefit because constructing P^-1 scales cubically.
Scalability: MD22 spans four biomolecular and supramolecular classes from 42-atom peptides to a 370-atom double-walled nanotube.The benchmark was designed to test scaling to systems containing several hundred atoms.
Scalability: Force RMSE was targeted at around 1 kcal mol^-1 Å^-1, requiring only a few hundred training points for some large systems but several thousand for others.The buckyball catcher and double-walled nanotube reach the target with a few hundred points, whereas DHA, stachyose, and Ac-Ala3-NHMe require several thousand.
Representation of non-local interactions: All atoms contribute to sGDML predictions, and ring-rotation energy changes are delocalized across the donor-bridge-acceptor molecule.This behavior is consistent with a global model rather than one partitioning energy into localized atomic contributions.

IV. MOLECULAR DYNAMICS

The study evaluates large-scale sGDML force fields in nanosecond molecular-dynamics and path-integral simulations. The simulations reach equilibrium, reveal rotation-driven fluctuations in the nanotubes, and show that nuclear quantum effects improve selected vibrational frequencies.

Simulation setup: Nanosecond classical MD and PIMD simulations of a hydrogen-saturated double-walled carbon nanotube were run at 300 K with 0.2 fs timesteps.The PIMD simulations used 16 beads and a Langevin thermostat.
Simulation stability: After thermalization, the simulations reach equilibrium at roughly 500 ps for classical MD and 100 ps for PIMD.Equilibrium was assessed using cumulative potential-energy behavior.
Geometry fluctuations: PIMD produces RMSD peaks exceeding 3.0 Å, driven by relative rotation of the inner nanotube with respect to the outer nanotube.The RMSD variations strongly correlate with the relative rotation angle Φ.
Geometry fluctuations: PIMD samples relative rotation angles from 40° to 80° more evenly, whereas classical MD shows pronounced peaks near 60° and 80°.The differences support smoothing of the nanotube rotational profile by nuclear quantum effects.
Vibrational spectra: Nuclear quantum effects shift the =C-H stretching peak closer to 3000 cm^-1, correcting the classical MD value of approximately 3100 cm^-1.Low-frequency long-range vibrational spectra are generally consistent between MD and PIMD.

V. CONCLUSION

The work enables global sGDML force fields to scale to larger systems and training sets without approximating the original model. It also creates opportunities for broader MLFF development and benchmarks.

V. CONCLUSION: An iterative scheme applies sGDML to significantly larger systems and training sets without introducing approximations to the original model.Numerical preconditioning reduces the learning problem’s conditioning number and enables rapid convergence of conjugate-gradient iterations.
V. CONCLUSION: Global interactions are represented on equal footing with local interactions, supporting long-timescale molecular dynamics for systems with far-reaching correlation lengths.
V. CONCLUSION: The MD22 benchmark introduces challenges involving molecular size and flexibility for atomistic models and future MLFF architectures.
V. CONCLUSION: Kernel-based MLFFs can use GPU parallelism and deep-neural-network software infrastructure, while kernel principles can inspire architectures such as transformers.

Appendix A: MD22 dataset properties

The MD22 datasets were generated with specified electronic-structure, simulation, thermostat, and sampling settings. The listed systems use PBE+MBD calculations and span several molecular complexes.

Appendix A: MD22 dataset properties: MD22 reference energies and atomic forces were calculated with FHI-aims at the PBE+MBD level of theory, using i-PI for molecular dynamics.The trajectories were sampled every 1 fs.
Appendix A: MD22 dataset properties: All listed MD22 systems were simulated at 500 K with either global Langevin or Nosé–Hoover thermostats.
Appendix A: MD22 dataset properties: Buckyball catcher and double-walled nanotube simulations used the light basis setting and Nosé–Hoover thermostats with coefficient 1700.
Appendix A: MD22 dataset properties: The dataset properties table organizes each system by dataset, formula, size, and energy and force ranges and variances.

Appendix B: Donor-bridge-acceptor dataset

The donor-bridge-acceptor reference data were generated by normal-mode sampling with randomized individual-bond rotations. Energies and forces used the semi-empirical GFN2-xTB method.

Appendix B: Donor-bridge-acceptor dataset: Reference configurations were generated by normal-mode sampling at 300 K with random rotations applied to individual bonds.
Appendix B: Donor-bridge-acceptor dataset: Energies and forces were calculated using the semi-empirical GFN2-xTB method.
Appendix B: Donor-bridge-acceptor dataset: The energy profile varies a single-bond rotation while all other bond distances and angles remain fixed at equilibrium values.

Appendix C: Scaling to larger training dataset sizes

The iterative solver extends kernel-based sGDML training to larger datasets while reducing memory demands. The authors use large training sets to demonstrate solver stability, not practical force-field reconstruction.

Appendix C: Scaling to larger training dataset sizes: 0.15 per mil of the closed-form solver’s memory was required for the largest molecule at 50k training points.This allowed the learning problem size to increase by two orders of magnitude on the same hardware.
Appendix C: Scaling to larger training dataset sizes: Models were trained for all MD17 molecules with dataset sizes of 1k and 50k using a preconditioner based on 75 inducing points.
Appendix C: Scaling to larger training dataset sizes: All models reached a training error of 10^-4 kcal mol^-1 A^-1 after approximately 800 to 2700 iteration steps.The range was approximately 800 steps for benzene to approximately 2700 for ethanol.
Appendix C: Scaling to larger training dataset sizes: The exercise demonstrates iterative-solver stability for large-scale problems rather than practical force-field reconstruction with massive training sets.The authors state that generating such reference datasets is infeasible for practical force-field reconstruction.

Appendix D: Parametric complexity of the models

sGDML force fields use non-parametric Gaussian-process models, yet their parameter counts are substantially lower than those of comparable neural-network force fields.

sGDML force fields are based on non-parametric Gaussian-process models whose parameter sets adapt to training-data complexity.Unlike parametric models, their parameter set is not fixed in advance.
In practice, sGDML force fields have around one order of magnitude lower parametric complexity than neural-network force fields at comparable accuracy.
The largest sGDML model for MD17-aspirin uses 63k parameters with 1k training points and 630k with 10k training points.The aspirin system contains 21 atoms.
Comparable neural-network force fields use 500k–3M parameters with 1k training points and 120k–12M with 10k training points on MD17.These ranges cover the cited NewtonNet, SpookyNet, NequIP, SchNet, and ForceNet models.

Appendix E: Cholesky decomposition

The appendix describes kernel-system solution and preconditioning strategies for Gaussian-process learning, combining Cholesky-based factorization with a Nyström approximation selected through leverage scores.

The Gaussian-process system α = (K + λI)−1y is solved through Cholesky decomposition into lower-triangular factors, followed by forward and backward substitution.
The Nyström method projects the symmetric positive-semidefinite kernel matrix onto selected columns to form a memory-efficient, easily invertible preconditioner.The approximation is used for Kλ in the learning problem.
Approximate λ-ridge leverage scores estimate column importance without directly inverting K + λI, whose inversion is as expensive as the original learning problem.
Inducing columns sampled from the leverage-score distribution approximate Kλ with constant probability and small relative error versus the best rank-k eigendecomposition approximation, at lower computational cost.

Appendix G: Implementation details

The implementation avoids storing the full kernel matrix, constructs the preconditioner through numerically robust factorizations, and uses initialization and restart choices to support stable iterative convergence.

Memory requirement: Kernel matrix–vector products are computed on the fly with NumPy LinearOperator, reducing stored Kαt memory complexity to O(m).The factorized P−1 is retained in memory.
Preconditioner construction: The preconditioner uses an indirect Cholesky decomposition through thin QR factorization because direct Woodbury evaluation can be numerically unstable.The indirect decomposition avoids squaring the condition number of the relevant term.
Preconditioner construction: Preconditioner construction maintains O(mk) memory complexity, and P−1 is applied through LinearOperator without expanding the matrix product.
Initial guess α0: The solver uses α0 = 0 because non-zero initialization can adversely affect convergence of Krylov subspace methods.
CG solver restarts: CG restarts are triggered when accumulated rounding errors make optimization steps non-orthogonal and slow or stall convergence.The implementation monitors average residual-reduction progress within a rolling time window.
Simulation diagnostics: The supplementary simulations track RMSD, cumulative potential energy, and molecular spectra for classical MD and PIMD trajectories.The cited figures cover RMSD distributions, energy convergence, and velocity-autocorrelation spectra for the buckyball catcher and double-walled nanotube.