MBD-ML: Many-body dispersion from machine learning for molecules and materials

AI summary

Overview

Research area: Computational physics and chemistry — van der Waals (vdW) dispersion methods, many-body dispersion (MBD), and machine-learned force fields (MLFFs).

Technical level: Advanced. The paper assumes familiarity with density functional theory, dispersion corrections (D3, D4, TS, XDM), the MBD Hamiltonian, equivariant message-passing neural networks, and periodic electronic-structure codes.

Scope (one sentence): The paper introduces MBD-ML, a pretrained message-passing neural network that predicts the atomic polarizability and C6 ratios needed to parametrize the many-body dispersion Hamiltonian directly from atomic structure, and validates it on molecules, dimers, molecular crystals, and inorganic materials.

What This Paper Is About

The many-body dispersion (MBD) method is one of the most accurate and transferable ways to capture vdW interactions, but it needs atomic C6 coefficients and polarizabilities as input, which normally come from a full electronic-structure calculation. That dependence on DFT is the main bottleneck keeping MBD out of large-scale simulations and machine-learning architectures.

This paper removes the bottleneck by training a neural network to predict those atomic inputs from geometry alone, then wiring the model into the libMBD library so that MBD energies, forces, and stress tensors can be produced with no electronic-structure calculation at all.

Key Contributions

1. A pretrained MBD-ML model. An SO3krates equivariant message-passing neural network is re-trained to predict the unit-less ratios α0^r and C6^r — the ratio of the Vydrov-Van Voorhis (VV) atomic coefficients in a molecule or solid to their free-atom values — using DFT+MBD-NL as reference. Unlike the antecedent machine-learned MBD model, which predicted Hirshfeld volume ratios for MBD@rsSCS and was limited to small C, H, N, O molecules, MBD-ML is applicable to both molecules and crystals across over 70 chemical elements.

2. Direct integration with libMBD. Multiplying the predicted ratios by free-atom reference quantities and feeding them to libMBD yields MBD-inclusive total energies, atomic forces, and stress tensors in a single function call through the pymbd Python interface, with no user intervention.

3. Broad validation across five state-of-the-art data sets. QCML (small molecules), OMol25 (large molecules), DES370k (biomolecular and drug-like dimers), OMC25 (molecular crystals), and OMat24 (inorganic materials).

4. A practical test on polymorph ranking and geometry. Crystal structures relaxed at PBE+MBD-ML are compared with PBE+MBD-NL reference structures, and polymorph energy orderings are benchmarked against PBE0+MBD-NL.

Main Findings

• Small molecules (QCML hold-out, N = 83226): RMSE of 0.020 for α0^r and 0.023 for C6^r; MBD energy RMSE 0.228 meV/atom and force RMSE 0.440 meV/Å.
• Molecular dimers (DES370k, N = 309130): errors comparable to or better than the QCML hold-out, with α0^r RMSE 0.025, C6^r RMSE 0.022, MBD energy RMSE 0.165 meV/atom, and force RMSE 0.387 meV/Å — evidence of transferability beyond the training distribution.
• Molecular crystals (OMC25, N = 200 unique-molecule crystals recomputed at PBE0+MBD-NL): ratio RMSE of 0.02–0.03; MBD energies and forces accurate to within 1 meV/atom and 1 meV/Å per atom; stress tensor components, which range only between 0 and 15 meV/ų, predicted with an RMSE of 0.143 meV/ų.
• Force quality versus pairwise methods (825 OMol25 molecules, 3 to 350 atoms): MBD-ML median and mean magnitude deviations of 2–3% with a 95th percentile of 7.7%, and angular deviations averaging below 2° with a 95th percentile of 5.1°. D3 shows 28.2% and 17.4°, and D4 shows 39.5% and 19.8°. Extreme outliers with magnitude deviations above 100% number 3066 atoms (2.4%) for D3 and 4025 (3.2%) for D4, versus 8 atoms (0.006%) for MBD-ML.
• Geometry optimization: PBE+MBD-ML reproduces PBE+MBD-NL crystal structures with RMSDs of 0.001 Å to 0.05 Å and unit-cell volumes differing by less than 1%. Uncorrected PBE overestimates volumes by over 20%, with RMSDs one to two orders of magnitude higher.
• Polymorph energy ranking: 7 of 9 transitions are ranked correctly with errors below 0.3 kJ/mol per molecule, spanning van der Waals, hydrogen-bonded, and mixed interaction types. Two rankings are incorrect: the I→II transition of 2-amino-5-nitropyrimidine (0.5 kJ/mol per molecule) and the I→III transition of maleic hydrazide (0.6 kJ/mol per molecule). Typical polymorph energy differences are on the order of 0–7 kJ/mol.
• Inorganic materials (OMat24 subset, N = 203, HSE06+MBD-NL reference): mean absolute errors of 0.247 for α0^r and 0.269 for C6^r, 10–20 times higher than for all other chemical classes studied. The errors are large enough to introduce false negative MBD eigenvalues, preventing energy and force calculation for most materials in the subset.
• Under-represented elements: alkali (Li, Na, K) and alkaline earth (Be, Mg, Ca) elements each appear in only 10⁻⁵–10⁻⁴ % of the 33.5 million QCML structures (roughly 3–30 structures per element), making them among the rarest elements alongside the lanthanides; these were excluded from the DES370k and OMol25 assessments.
• Anions: the α0^r and C6^r ratios of negatively charged molecules vary substantially with the confinement potential cutoff radius in FHI-aims, and random inspection showed the pathological anions had negative electron affinities. Including anionic molecules deteriorates ratio predictions by more than one order of magnitude, while MBD energies and forces remain nearly unchanged. Negatively charged molecules were therefore excluded from validation.
• Error origin in crystals: a sensitivity analysis found that systematic bias in the predicted ratios influences MBD energies roughly twice as strongly as random noise. Because unit cells often exceed 100 atoms, a per-atom error of 0.8 meV can translate into a total-energy deviation exceeding 80 meV, but such errors cancel almost exactly for energy differences.
• Computational scaling: for water clusters (H2O)n with n = 3–4321 (up to nearly 13,000 atoms) run serially on a single 128-core AMD EPYC Rome node, ratio prediction takes only a few seconds for systems up to 1000 atoms. MBD energy and force evaluation shows the predicted cubic scaling beyond 60 atoms and is faster than the ratio calculation up to a crossover at 250 atoms. Above 1000 atoms, total execution time grows with an approximate power of 1.6, reaching a maximum of 196 seconds for (H2O)4321.
• Damping parameters: β was not re-optimized for MBD-ML; tabulated defaults were used (β_PBE = 0.81, β_PBE0 = 0.83).

Methodology in Plain English

The MBD method treats each atom's valence electrons as a charged particle attached to a positively charged pseudonucleus by a spring — a coupled quantum Drude oscillator. The whole system of oscillators is described by a Hamiltonian whose inputs are, per atom, a static polarizability α0 and an effective excitation frequency ω derived from C6. Diagonalizing that Hamiltonian gives the dispersion energy, and the computational cost scales as the cube of system size, though with a much smaller prefactor than Kohn-Sham DFT.

The catch is that α0 and C6 have traditionally been computed from an electronic density. The authors instead train a neural network to predict two unit-less ratios per atom: the VV functional's α0 and C6 for the atom in its environment, divided by the corresponding free-atom values. These ratios typically fall in the range of roughly 0 to 2 and are relatively short-ranged, which suits a semi-local message-passing architecture; they are also comparatively insensitive to the choice of exchange-correlation functional.

The model is an SO3krates network with two message-passing layers and a 4 Å cutoff radius, giving an effective receptive field of 8 Å, trained on the QCML data set of over 30 million molecules with up to 8 non-hydrogen atoms across 79 chemical elements. It takes atomic positions, atomic numbers, total charge, and spin as input. At run time, libMBD multiplies the predicted ratios by the free-atom reference values and proceeds with the standard MBD machinery. The authors then benchmark against ab initio MBD-NL reference values on held-out molecules, dimers, crystals, and materials, and separately compare geometry optimizations, polymorph rankings, and force deviations against D3 and D4.

Why This Matters

Impact on research. The accuracy of MBD-NL has until now come at the cost of a full electronic-structure calculation. MBD-ML decouples many-body dispersion from any underlying electronic-structure code, which makes state-of-the-art vdW treatment available to empirical force fields, machine-learned force fields, and high-throughput data-generation workflows. It also sidesteps the documented limitations of the pairwise D3 correction that currently dominates large-scale applications.

Real-world applications:

• Pharmaceutical and organic-crystal polymorph screening, where relative stabilities of molecular crystal forms often fall below 1 kJ/mol per molecule.
• Drug design and biomolecular modeling, where dispersion governs protein folding and molecular recognition.
• Catalysis and surface chemistry, where adsorption energies depend on dispersion.
• Battery and energy-materials applications, where layered, polymeric, and porous material cohesion is dispersion-driven.

Industry relevance. Any simulation stack built on FHI-aims, ASE, or a DFT code can call the pymbd interface with a single function call, and any MLFF pipeline can in principle adopt the model for a more accurate dispersion term. The OMC25 polymorph results — structural RMSDs of 10⁻³ to 10⁻² Å and ranking errors of at most 0.6 kJ/mol — are directly relevant to industries that depend on crystal-structure prediction.

Future Directions

• Extend training data to inorganic solids. The OMat24 results (MAE of 0.247 and 0.269, false negative MBD eigenvalues) point to a mismatch between the extended periodic environments of inorganic materials and the organic molecules in QCML; representative inorganic subsets would need to be added, which the authors state will be pursued in future work.
• Cover alkali and alkaline earth elements. These are among the rarest elements in QCML, and the authors state the limitation can be overcome by fine-tuning the model on a specific data set of interest.
• Establish reliable training data for molecular anions. The pathological behavior observed is a known difficulty of describing negatively charged molecules that is exacerbated within the MBD formalism; it will require identifying suitable basis set settings and density functional approximations.
• Re-optimize damping parameters and reuse the model with other functionals. Since the ratios are only weakly sensitive to the density functional approximation, applying MBD-ML with other functionals currently requires only changing the damping parameter β, which was held at tabulated default values here.

Target Audience

Researchers and practitioners in computational chemistry, condensed matter physics, and materials science who work with dispersion corrections, machine-learned force fields, or molecular crystal structure prediction. It is also relevant to developers of electronic-structure codes and simulation libraries who want to add many-body dispersion without a DFT dependency, and to method developers interested in how equivariant neural networks can target physically motivated intermediate quantities rather than energies and forces directly.

Authors’ abstract

Van der Waals (vdW) interactions are essential for describing molecules and materials, from drug design and catalysis to battery applications. These omnipresent interactions must also be accurately included in machine-learned force fields. The many-body dispersion (MBD) method stands out as one of the most accurate and transferable approaches to capture vdW interactions, requiring only atomic $C_6$ coefficients and polarizabilities as input. We present MBD-ML, a pretrained message passing neural network that predicts these atomic properties directly from atomic structures. Through seamless integration with libMBD, our method enables the immediate calculation of MBD-inclusive total energies, forces, and stress tensors. By eliminating the need for intermediate electronic structure calculations, MBD-ML offers a practical and streamlined tool that simplifies the incorporation of state-of-the-art vdW interactions into any electronic structure code, as well as empirical and machine-learned force fields.