libMBD: A general-purpose package for scalable quantum many-body dispersion calculations

,


Introduction
Van der Waals (vdW) dispersion interactions stem from electron correlation induced by the long-range part of the Coulomb force.Since they are attractive under normal circumstances, their importance grows with the characteristic length scale in an atomistic system, making them an often decisive force in large molecular complexes, molecular crystals, nanostructured materials, surface phenomena, and soft matter in general.At the same time, the workhorse electronic structure method for simulating molecules and materials-Kohn-Sham density-functional theory (KS-DFT) in semilocal and hybrid approximations-neglects vdW interactions by construction.As a result, a large number of approaches to mitigate this issue have been developed, and some of them are now used routinely when performing DFT calculations (Grimme et al. 2016;Hermann et al. 2017b).
One such popular approach is the many-body dispersion (MBD) method, which is based on a model Hamiltonian of charged harmonic oscillators that capture the long-range electrodynamic response of an atomic system and yield its vdW energy (Tkatchenko et al. 2012).Its essential offering is that the coarsegraining of the full electronic structure to oscillators makes it efficient, while full many-body treatment of the electronic fluctuations provides consistent description of beyond-pairwise vdW interactions (Ambrosetti et al. 2016).Both these characteristics make the MBD framework especially apt for large-scale model-ing of molecular systems, for which the demand is bound to only increase in the future (Akimov & Prezhdo 2015).Given the existence of MBD integrations with even more efficient substitutes of DFT based on density-functional tight binding (DFTB) (Stöhr et al. 2016) and machine learning (Bereau et al. 2018;Poier et al. 2022a; Unke et al. 2021;Westermayr et al. 2022), as well as the fact that MBD is an approach that is still evolving and being improved and extended (see below), the situation calls for an efficient, flexible, and reusable implementation of MBD that can be easily deployed in any given electronic structure code.Here, we develop and document libMBD, a library package that answers this call.libMBD is written in modern Fortran with additional bindings available for C and Python.It uses MPI/ScaLAPACK parallelization, and implements MBD for both finite and periodic systems, with analytical gradients of the energy with respect to all input parameters, which enables both structure optimization as well as self-consistent vdW energy with respect to the electron density (Ferri et al. 2015).libMBD is designed as a framework for any method based on the MBD Hamiltonian, and has a built-in support for the original MBD formulation (Tkatchenko et al. 2012), the range-separated self-consistently screened variant (MBD@rsSCS) (Ambrosetti et al. 2014), the metal-surface parameterizations (MBD surf ) (Ruiz et al. 2016(Ruiz et al. , 2012)), as well as the recent hybrid of MBD with nonlocal functionals (MBD-NL) (Hermann & Tkatchenko 2020).For convenience, an efficient implementation of the predecessor of MBD, the pairwise TS model (Tkatchenko & Scheffler 2009), is included as well.Finally, libMBD offers not only the vdW energy and its gradi-ents, but can calculate also other quantities derived from the MBD Hamiltonian, namely the eigenstate charge fluctuations and induced charge-density polarization (Hermann et al. 2017a), firstorder Coulomb corrections (Stöhr et al. 2021), the dielectric function (Cui et al. 2020), and the order-by-order expansion of the energy in the dipole interaction (Tkatchenko et al. 2013).
The popularity of MBD as a method is revealed through its independent implementations in numerous electronic structure codes, namely FHI-aims (Blum et al. 2009), VASP (Kresse & Furthmüller 1996), Q-Chem (Epifanovsky et al. 2021), Quantum ESPRESSO (Giannozzi et al. 2020), CASTEP (Clark et al. 2005), and DFT-D4 (Caldeweyher et al. 2019).The original standalone reference implementation of MBD is available for download via Internet Archive (MBD).While some of these implementations introduced methodological innovations-such as the one in VASP, which introduced a reciprocal-space formulation (Bučko et al. 2016), and the one in Quantum ESPRESSO, which introduced analytical gradients (Blood-Forsythe et al. 2016)-all of them lack both in computational efficiency and in capability compared to libMBD.The single exception is the linear-scaling stochastic evaluation of MBD (Poier et al. 2022b), which is more efficient for very large systems, but for the moment lacks analytical gradients.Furthermore, none of these MBD implementations have been designed as a library and their reuse across different electronic structure codes is problematic in theory, and nonexistent in practice.In contrast, libMBD has been already integrated into several such codes, including FHI-aims, VASP, Q-Chem, Quantum ESPRESSO, DFTB+ (Hourahine et al. 2020), and ASE (Larsen et al. 2017).
The flexibility of MBD as a general approach can be seen in the ever-growing list of methods and models that are based on it, and for this reason libMBD has been designed from the start as a framework that is easy to extend and modify.Besides the methods listed above that are already incorporated in libMBD, there are the fractionally ionic variant (MBD/FI) (Gould et al. 2016) and its successor "universal" MBD (uMBD) (Kim et al. 2020), the Wannier-function parameterization (vdW-WF2) (Silvestrelli 2013), beyond-dipole extension (Massa et al. 2021), anisotropic extension (Bandyopadhyay et al. 2022), extension to excited states (Ambrosetti et al. 2022), and the MBD variant of DFT-D4 (Caldeweyher et al. 2019).If implemented in or with libMBD, these methods and models or any future ones would become instantly available across all the electronic structure codes to which libMBD has already been or will be integrated, and could take advantage of parallelization, analytical gradients, and both finite and periodic formulations.Next to computational efficiency, this is the biggest offering of libMBD.
The rest of the paper is structured as follows: Section 2 reviews all the MBD methods implemented in libMBD in terms of equations, as well as the physics of integrating MBD with DFT and its substitutes.Section 3 documents the structure of the code and the public interface served to the electronic structure codes, its scaling properties, numerical convergence, functionality, and existing interfaces.Section 4 describes a number of example MBD calculations and a meta-review of important past applications of MBD.The paper is concluded in Section 5.

Theory
This section reviews all mathematical equations that define most of the functionality of libMBD.It is structured such that equations in Sections 2.1 to 2.7 are implemented in libMBD, while equations in Section 2.8 must be implemented in the codes with which libMBD is integrated, since their implementation is inherently code-specific.

Many-body dispersion
Any MBD method maps an atomic system to a model Hamiltonian of charged harmonic (Drude) oscillators located at   , with static polarizabilities  0, and frequencies   .In most variants, the oscillators represent atoms, but both finer fragments, such as Wannier functions (Silvestrelli 2013), and coarser fragments, such as whole molecules (Jones et al. 2013), are also possible.As a result of the coarse-graining, only the long-range part of the interaction between oscillators is expected to represent faithfully the original electronic system, justifying a dipole approximation.The MBD Hamiltonian for a finite system of  oscillators then reads Here,   = √   (  −   ) are displacements of the oscillating charges weighted by the oscillator masses   , and  lr is a damped dipole interaction tensor, parametrized through some single-oscillator properties   .The mass is in principle the third parameter (after static polarizability and frequency) required to fully specify a charged quantum oscillator, but it has no effect on the energy under the dipole approximation.Other sets of three independent oscillator parameters can be used as well, for example the  6 dispersion coefficient or charge , but these are always connected through simple relations, such as The damped dipole tensor must reduce to the "bare" dipole tensor  for large distances, Within this general MBD framework, a particular MBD model is then specified through the level of coarse-graining (typically atoms), the parametrization of the oscillators, and the choice of  lr , which necessarily depends on the effective range (degree of non-locality) of the electronic structure method to which MBD is coupled.
A major consequence of the dipole approximation is that the MBD Hamiltonian can be exactly diagonalized, after which it describes a system of uncoupled collective oscillations ξ = ∑   ,   .Here, we use a shorthand notation in which the particle and Cartesian indexes are joined, for example (  )  =   = ()  and (  )  =  , = () , .The MBD energy is then obtained as the change in the zero-point energy of the oscillations (fluctuations) induced by the dipole interaction, with  lr  ≡  lr (  ,   ,   −   ).Note that the diagonalization leads to complex-valued energy if some of the eigenvalues of  are negative.This can result from a failure of the dipole approximation, and in typical scenarios signifies that the particular parametrization of the MBD Hamiltonian does not properly represent the local polarization response of a given molecule or material.

Periodic boundary conditions
In a periodic system, the collective oscillations have an associated wave vector  from the first Brillouin zone (FBZ) as a good quantum number, and the energy per periodic unit cell can be calculated as an integral over the FBZ, evaluated in practice as a sum over a -point mesh, MBD () is obtained from Eqs. ( 5) to ( 7) by making ω,  lr , and  formally dependent on  and considering The infinite sum of the dipole tensor over all periodic copies of the unit cell, indexed by , and skipping the  = ,  =  terms, is evaluated efficiently with Ewald summation (de Leeuw et al. 1980a), Here,  indexes the cells of the reciprocal lattice,   is a reciprocal-lattice vector,  c and  c are the real-space and reciprocal-space cutoffs, respectively,  is the Ewald rangeseparation parameter (see Section 3.3), Ω is the unit-cell volume, and  erfc is the short-range part of the bare dipole tensor, The last term in Eq. ( 10) is the so-called surface term (Ballenegger 2014;de Leeuw et al. 1980b), which does not contribute to the energy (in practice the -point mesh avoids the Γ-point) or most other physical observables, and is stated here only for completeness.

Self-consistent screening
When put in an external electric field  ext , a polarizable system partially (or fully) screens the field through induced polarization that is self-consistent with the total field.For charged harmonic oscillators interacting via some dipole interaction tensor  scs , this reads or equivalently in the shorthand notation Here, ᾱ is a nonlocal anisotropic polarizability that describes the dipole on the -th oscillator induced by an external field on the -th oscillator.
In the MBD@rsSCS method, self-consistent screening is used to obtain a refined set of oscillator parameters before they are used in the MBD Hamiltonian.To do this, one can assume a homogeneous external electric field and contract ᾱ to effective local polarizability, and average out the angular dependence to obtain isotropic polarizability, Effective screened oscillator frequencies ω are obtained through by screening the dynamic oscillator polarizability evaluated at imaginary frequency, In practice, the isotropic local polarizability in Eq. ( 17) is evaluated on a quadrature grid of imaginary frequencies, screened through Eq. ( 14) at each grid point, contracted via Eq.( 15), and the integral in Eq. ( 16) is calculated numerically (see Section 3.3) to yield the  6 coefficients and oscillator frequencies.

Analytical gradients
Analytical formulas for the gradients of the MBD energy with respect to all Hamiltonian input parameters can be obtained straightforwardly by application of the derivative chain rule.The two steps related to matrix diagonalization in Eq. ( 6) and matrix inversion and contraction in Eqs. ( 14) and ( 15), which are nontrivial, and for which naive application of the chain rule could result in inefficient evaluation of the gradients are explicitly stated below.When implemented, these formulas can be straightforwardly parallelized and result in a computational cost for the gradients that scales with the same power of system size as the energy evaluation.
The gradient of the MBD energy with respect to any Hamiltonian oscillator parameter   is The gradient of the contracted SCS polarizability with respect to any oscillator parameter   is

MBD properties
The normalized ground-state wave function of the molecular MBD Hamiltonian can be expressed as While only two parameters per oscillator,  0, and   , are needed to get the MBD energy, a third parameter, such as   or   , is needed to evaluate the wave function in real space, in order to convert from   to   .The wave function of the non-interacting oscillators is simply The MBD wave function has been used to calculate two properties of interest-the polarization of the electron density due to vdW interactions (Hermann et al. 2017a), and the first-order correction to the MBD energy from the Coulomb interaction (Stöhr et al. 2021), All these expectation values can be evaluated analytically, and the derivations and final results can be found in Hermann (2018).Instead of diagonalizing the MBD Hamiltonian, the MBD energy can be equivalently obtained as an integral over the imaginary frequency, which also allows for a many-body expansion of the energy in the orders of  lr (Tkatchenko et al. 2013), In practice,  lr can be replaced with the Hermitian  1∕2  lr  1∕2 since both are equivalent under the trace operator.This integral formulation has two applications.First, one can analyze the contribution of the second-order (pairwise) and higher many-body orders to the vdW energy.Second, the eigenvalues   of (i) lr can be heuristically renormalized to obtain a real-valued MBD energy even when the exact diagonalization gives a complex-valued energy due to the dipole collapse (Gould et al. 2016).In particular, one replaces Another property of interest is the macroscopic dielectric constant  M of a periodic system of charged oscillators, which can be calculated via SCS from ᾱ as

Pairwise dispersion
The second (lowest) order of the many-body expansion in Eq. ( 25) corresponds to the familiar pairwise expression for the vdW energy, Here,  is a damping function ( → 1 when   → ∞) and is typically specified directly rather than through  lr .This ansatz has been used in numerous vdW models, among them the predecessor of MBD, the TS method (Tkatchenko & Scheffler 2009).Under periodic boundary conditions, the corresponding lattice sum of Eq. ( 28) converges absolutely, unlike the dipole tensor, but the rate of convergence can be slow enough to present a computational bottleneck.An alternative is to use the Ewald technique (Karasawa & Goddard 1989), with the same notation as in Eq. ( 10).

MBD methods
The general MBD framework discussed up to this point is concretized into an MBD method by specifying at least two aspects.First, how are the MBD oscillators parametrized ( 0 , ).Second, how is  lr chosen to smoothly integrate the long-range correlation from MBD with a given method that accounts for the short-range correlation, typically KS-DFT or its substitute.While independent to a certain degree, these two choices are ultimately related in practice.

TS
In the pairwise TS method (Tkatchenko & Scheffler 2009), each atom is mapped to a single oscillator parametrized through scaling of reference free-atom vdW parameters based on atomic volumes, Here,   is some measure of a volume of an atom in a molecule or material and  ref  is the same measure for a corresponding free atom.Possible choices for the volume measure are discussed in Section 2.8.
The logistic sigmoid function parametrized with volumerescaled vdW radii  vdW is used as the damping function in Eq. ( 28), The damping parameter   adjusts the onset of the vdW interaction and is optimized separately for each individual short-range correlation model, and  = 20.

MBD@rsSCS
The MBD@rsSCS method (Ambrosetti et al. 2014) is the many-body extension of the TS method, in which the oscillators interact through the dipole tensor derived from the full Coulomb interaction of two oscillator (Gaussian) charge densities with widths  2  , The oscillator width is related to its polarizability through certain self-consistence requirements (Mayer 2007), To integrate this full-range interaction with KS-DFT while avoiding double-counting of the short-range correlation,  GG is rangeseparated, and the short-range part is used in SCS to refine the TS oscillator parameters from Eq. ( 30 Here,  is a damping function of the same form as  in Eq. ( 31), but with  = 6 and   usually denoted as .

MBD-NL
The MBD-NL method (Hermann & Tkatchenko 2020) removes the need for SCS, and the MBD Hamiltonian is parametrized directly with Here, a polarizability functional of the density is used to obtain vdW parameters ( df 0 ,  df ) for an atom in a molecule or material and the corresponding free atom (see Section 2.8).MBD-NL uses the same parametrization of  lr in Eq. ( 35) as MBD@rsSCS, but with vdW radii in Eq. ( 31) replaced with

Integration with short-range models
MBD models the long-range part of electron correlation and must be always coupled with a model that captures the rest of the electronic energy, including the short-range correlation, Typically, this short-range model is some semilocal or hybrid version of KS-DFT, but it can also be its computationally more efficient substitute such as DFTB or a machine-learned interatomic potential.In any case, the short-range model and MBD must be integrated in two ways.First, the damping parameter in an MBD method (Section 2.7) must be adjusted to avoid double counting of correlation in the intermediate range.This is typically done by fitting on the S66x8 dataset of vdW dimers of small organic molecules (Brauer et al. 2016).Second, the short-range model is used to parametrize the MBD Hamiltonian, that is, to obtain the vdW parameters of each oscillator.How this is done depends on a particular model and the various options are described below.In the context of libMBD, the following equations must be necessarily evaluated in the code that implements the short-range model and couples to libMBD.
Hirshfeld volumes This parametrization uses Eq. ( 30), in which the atomic volumes are calculated from the Hirshfeldpartitioned electron density obtained from a DFT calculation, Here  ref  are electron densities of free atoms located at  = 0.A straightforward extension is to use not only neutral free atoms, but also free ions (Bučko et al. 2014) for the density partitioning, which improves description of ionic compounds.This was further extended to use not only ion densities, but also ionic reference vdW parameters by using a piecewise linear dependence of atomic polarizabilities on the charge (Gould et al. 2016).This results in atomic   (i) that cannot be expressed through a single oscillator with Eq. ( 17), requiring either direct use of Eq. ( 24) or calculation of effective oscillator frequencies via Eq.( 16).

MBD-NL
Electron density from a KS-DFT calculation is also used in MBD-NL to parametrize the MBD Hamiltonian.Here, a renormalized VV10 polarizability functional of the density (Vydrov & Van Voorhis 2010) is coarse-grained via Hirshfeld partitioning, and the resulting atomic dynamic polarizabilities are reduced to effective oscillator frequencies ( df  ) via Eq.( 16).The freeatom reference values ( df,ref 0, ,  df,ref

𝑖
) are obtained by evaluating the functional on free-atom densities.The renormalized VV10 functional uses a semilocal cutoff function  (see Hermann & Tkatchenko 2020, Eq. 7), which is a function of the local ionization potential  and the iso-orbital indicator , which in turn are functions of the electron density, its gradient, and the kinetic energy density, quantities readily available in any code that implements meta-GGA functionals.

Charge population analysis
In the absence of a real-space representation of the electron density, the above-mentioned methods are not applicable.Many electronic structure methods exist where the electron density is not at all or only rarely evaluated explicitly.This includes approximate methods such as DFTB, as well as other density-matrix-based approaches.Stöhr et al. (2016) found that a correlation between effective atomic polarizability of an atom in a molecule can be established based not only on the effective volume of an atom, but also on measures derived from the density matrix represented in a local atomic orbital basis set, |  ⟩ = ∑     |  ⟩, where single-electron states |  ⟩ are expanded in atomic orbital basis functions |  ⟩.When scaling freeatom vdW parameters (Eq.30) with the degree of hybridization between atoms as measured by the sum over the on-site component of the Mulliken charges (or equivalently the atom-projected static atomic polarizabilities can be predicted that are similar to the ones obtained by rescaling with Hirshfeld volumes (Eq.40).

Machine-learned interatomic potentials Several works have
shown that the effective atoms-in-molecules volume scaling ratios as a function of the atomic positions can be accurately represented by machine-learning regression (Bereau et al. 2018;Poier et al. 2022a;Westermayr et al. 2022).Once a continuous representation of the scaling ratio is constructed, an MBD description of the long-range correlation energy can be straightforwardly coupled to short-range force fields or atomistic machine-learned potentials.
In cases where these potentials are directly trained on semilocal approximations to DFT, the choice of the range-separation parameter remains identical to the optimal choice for the underlying density-functional approximation.

Code structure and API
libMBD implements all equations in Sections 2.1 to 2.7 in Fortran 2008.To enable straightforward and rapid development of new variants and extensions of MBD, libMBD has been designed from the start to maintain a close correspondence between code and physics (Fig. 1).Equations are documented in source code with L A T E X right next to the individual procedures that implement them and rendered in auto-generated documentation (https://libmbd.github.io).Gradients of the energy are computed in a forward mode, where each individual function that Eqs. ( 24) to ( 26) mbd_methods.f90mbd.f90 Eq. ( 30) Eqs. ( 14), ( 15), ( 19) Eq. ( 36) Eqs. ( 2), ( 16), ( 17), ( 33) for each q-point in Brillouin zone (Eq.8): Eqs. ( 5) to ( 7), ( 18 Figure 1 | Mapping code to physics.libMBD was designed for ease of development, which is achieved through tight correspondence between code and physics.An overview of an MBD@rsSCS calculation is shown on the background of a map of MBD source files with the individual implemented equations. computes a value from its arguments also optionally returns the partial derivatives of the output value with respect to some (or all) of its input arguments.Parallelization is available either by trivial distribution of the -point mesh (Eq.8) or by distribution of atom-indexed matrices (, , ᾱ) via BLACS.Linear-algebra operations are handled via LAPACK, ScaLAPACK, or ELPA (Marek et al. 2014)  API (Listing 1) is based on two derived types, mbd_input_t mbd_calc_t, where the former is a plain data structure of all input parameters used to initialize the latter, which is an opaque handle that provides access to computation and results.The Python API (Listing 2) is based on the MBDGeom class that represents a given geometry and provides access to computation via its methods.80% of the code base of libMBD is covered by tests, which include both unit tests of analytical gradients of individual procedures as well as regression tests of the total energies and gradients.libMBD is built with CMake.

Scalability
For libMBD to be practical, it needs to keep up with the codes that implement the corresponding short-range models (DFT and its substitutes) in terms of scalability to large systems, and scalability to a large number of processors.libMBD uses only process parallelization through the MPI and BLACS, and no thread parallelization.
Tested on various supercells of crystalline urethane (unit cell with 26 atoms), with or without periodic boundary conditions, the MBD method and its implementation in libMBD scale with the number of atoms, , as follows (Fig. 2).For small , the computational cost is dominated by the construction of the vari- ous matrices and as such grows quadratically with .At around a few hundred atoms (finite system) or a few thousand atoms (periodic system), the matrix operations that scale cubically with  (multiplication, diagonalization, inversion) start to dominate.The MBD analytical gradients as implemented by Blood-Forsythe et al. (2016) add one extra factor of  to the scaling (making it quartic) due to explicit evaluation of the gradient of ᾱ (see Eqs. 44 and 45 therein).libMBD avoids this by calculating only the gradient of the fused evaluation of ᾱ and its contraction to atoms (Eq.19), and as a result calculation of the gradients changes only the prefactor in the scaling behavior.Tested on a 5 × 5 × 5 supercell of urethane (3,250 atoms), libMBD exhibits good strong scaling with the number of processor cores  (Fig. 3).The construction of the various matrices (, , ᾱ) is trivially parallelizable and has perfect strong scaling up to  = 4096.The evaluation of the MBD gradients requires only very little inter-process communication and shows similarly good behavior.The distributed matrix inversion (SCS) and diagonalization (MBD) as provided by ScaLAPACK and ELPA exhibits strong scaling up to  = 256 on our particular setup, at which point the scaling plateaus.The evaluation of SCS gradients requires only matrix multiplications and contractions, but its current simple implementation in libMBD requires an unnecessary amount of inter-process communication, making the strong scaling plateau already at  = 64 for finite systems.Note that this slight inefficiency still allows the calculation of the MBD energy and gradients in ∼100 s at  = 256, and furthermore it affects only the MBD@rsSCS method, not the newer MBD-NL method which does not use SCS.Nevertheless if better strong scaling is required in the future, the inter-process communication in the evaluation of the SCS gradients can be massively optimized.

Convergence and default parameters
There are a few critical parameters that determine the numerical accuracy of the calculated MBD energy with respect to the exact energy as defined in Sections 2.1 to 2.7.The SCS is evaluated and integrated on an imaginary-frequency grid (Eq.16), for which libMBD uses a Gauss-Legendre quadrature transformed from The MBD energy converges exponentially with the grid size, and using the default 15 grid points, it is converged to relative error of 10 −8 (Fig. 4a).The frequency grid size is the only parameter for calculations on finite systems.For periodic systems, there are additional parameters related to converging the infinite lattice sums involved.
The MBD energy per unit cell of a periodic system is calculated by integration over the FBZ (Eq.8).libMBD uses a uniform quadrature grid in the reciprocal space, offset from the origin since the MBD energy has a removable singularity at the Γ-point (see Section 2.2).The convergence of the MBD energy with the -point grid size is polynomial and system-dependent (see Fig. 4b Table 1 | Available features in libMBD.

Python bindings periodic parallel gradients
✓  Includes forces, stress tensor, and energy gradients with respect to the vdW parameters (for self-consistence).
for example convergence on urethane), and as such the grid needs to be specified by the user.
The Ewald summation of the dipole interactions between the origin cell and the infinite crystal (Eq.10) is controlled by three interconnected parameters-the range-separation parameter , and the real-and reciprocal-space cutoffs  c ,  c .libMBD uses the following default prescription for the parameters, where Ω is the unit-cell volume.The MBD energy converges exponentially with the cutoffs and the default values ensure convergence to relative error of 10 −10 (Fig. 4c).The range-separation parameter does not affect accuracy, but through its determination of the cutoffs affects computational cost.The default value in libMBD ensures an optimal balance between the costs of the realand reciprocal-space parts (Fig. 4d).

Functionality
libMBD implements all equations in Sections 2.1 to 2.7, with all physical quantities accessible via the Python bindings (Table 1).Namely, this includes the calculation of MBD and TS energies, evaluation of SCS, evaluation of the MBD energy via the imaginary-frequency integration (Eq.24), access to MBD eigenstates, and evaluation of the MBD polarization density and Coulomb correction.All energy evaluations are implemented for both finite and periodic systems, parallelized, and with available analytical gradients (except for the imaginary-frequency integration).Evaluation of MBD properties such as the Coulomb correction or density polarization is currently implemented only for finite systems and is not parallelized.

Existing interfaces
libMBD was originally developed alongside FHI-aims (Blum et al. 2009) and was tightly integrated into it.Due to specific build requirements of FHI-aims, libMBD is currently distributed with FHI-aims in a lightly pre-processed form.The current API of libMBD was then based on the preexisting framework for interfacing external libraries in DFTB+ (Hourahine et al. 2020), which was the first third-party program in which libMBD was integrated.As a result, the DFTB+ integration is very lightweight and can serve as a model for interfaces with other programs.libMBD has been since integrated into Quantum ESPRESSO (Giannozzi et al. 2020), VASP (Kresse & Furthmüller 1996), and Q-Chem (Epifanovsky et al. 2021) without having to change anything in the existing source code, demonstrating its universality.While libMBD is typically integrated into a larger electronic structure program as part of a DFT or similar calculation, this is not the only use case.Especially when developing new methods or when investigating unexpected behavior of existing methods, it is practical to be able to evaluate MBD separately from DFT, and to be able to do so as flexibly as possible.For this reason, libMBD was from the start developed together with a Python interface called pyMBD.Beyond using pyMBD directly (Listing 2), it can also serve as a bridge between libMBD and other Python programs.One such existing example is the interface to ASE (Westermayr et al. 2022).

Applications
To demonstrate the functionality of libMBD, here we report on novel calculations done on supramolecular complexes and review selected previous applications of MBD.

Interaction energies
We apply the MBD scheme via the libMBD interface within the FHI-aims (Blum et al. 2009) software package to benchmark the interaction energies of structures in the L7 dataset (Sedlak et al. 2013) against existing data in the literature.We compare the interaction energies of the L7 complexes as calculated with MBD against higher-level theories, in particular diffusion Monte Carlo (DMC) (Foulkes et al. 2001;Reynolds et al. 1982) data from Al-Hamdani et al. (2021) and Benali et al. (2020) and local natural orbital coupled cluster with single, double, and perturbative triple excitations [LNO-CCSD(T)] (Rolik & Kállay 2011;Rolik et al. 2013) data from Ballesteros et al. (2021).Both DMC and LNO-CCSD(T) have been shown to accurately predict intermolecular interactions of small organic molecules.We also compare MBD-calculated interaction energies against DFT-D3 ( Grimme et al. 2010) data from Xu et al. (2021) and DFT-D4 data from Caldeweyher et al. (2017Caldeweyher et al. ( , 2019)).Interaction energies were calculated as the difference between the total energy of the structure and the sum of the monomer energies.
DFT calculations were performed with the FHI-aims software package, using the PBE (Perdew et al. 1996) and HSE06 (Heyd et al. 2003;Krukau et al. 2006) density-functional approximations.The standard screening parameter of 0.11  0 −1 was used for HSE06, and the numeric atomic orbitals were represented using a 'tight' basis set (Blum et al. 2009;Havu et al. 2009).The total energy, sum of eigenvalues, and charge density criteria were set to 1 × 10 −6 eV, 1 × 10 −2 eV, and 1 × 10 −5 eV  0 −3 , respectively.Input and output files for all DFT+MBD calculations are available as a dataset in the NOMAD electronic structure data repository (Chaudhuri et al. 2023).
DFT+MBD performs fairly well for the supramolecular complexes with respect to the higher-level LNO-CCSD(T) and DMC methods (Table 2).For the majority of L7 complexes (C3A, C3GC, GCGC, GGG, and PHE), DFT+MBD is within chemical accuracy (±1 kcal mol −1 ) of LNO-CCSD(T) and/or DMC.For complexes where the difference between dispersioncorrected DFT and higher-level interaction energies is greater than 1 kcal mol −1 , both PBE+MBD and HSE06+MBD generally outperform PBE+D3 with respect to the higher-level methods, as can be seen for both C2C2PD and CBH.
LNO-CCSD(T) and FN-DMC are in close agreement for five of the L7 structures (Table 2).However, their interaction en-ergies differ by 1.1 kcal mol −1 for C2C2PD, 2.2 kcal mol −1 for C3GC, and 7.6 kcal mol −1 for the C 60 @[6]CPPA complex.Interestingly, PBE0+D4 is in close agreement with LNO-CCSD(T) across all complexes, including the three aforementioned structures, with a mean absolute deviation of 1.1 kcal mol −1 .In contrast, all the DFT+MBD approaches (PBE+MBD, PBE0+MBD, and HSE06+MBD) are much closer to FN-DMC, especially for the C 60 @[6]CPPA complex.The differences between LNO-CCSD(T) and DFT+MBD, and between FN-DMC and DFT+D for the C2C2PD, C3GC and C 60 @[6]CPPA complexes makes it difficult to evaluate a reliable reference interaction energy for these three complexes.This disparity should motivate the continued development of dispersion-correction methods in order to better characterize the interactions within dispersion-dominated complexes.

MBD properties
The MBD method is based on an effective Hamiltonian and as such provides access to more than just energies (Section 2.5).The coupled zero-point oscillations of the electron density and the density polarization have been used to reinterpret - stacking (Carter-Fenk & Herbert 2020;Hermann et al. 2017a).Both can be readily calculated with libMBD using the pyMBD interface.
The C3A dimer from the L7 dataset exhibits typical behavior of these quantities in a - complex.The low-energy zero-point MBD oscillations display long-range order spanning the whole complex, and corresponding to simple electrostatic patterns such as a dipole moment across a whole molecule interacting with an opposite dipole moment on the other molecule (Fig. 6a).The vdW interaction in general polarizes the electron density out-ofplane in - complexes.The C3A dimer demonstrates that while  the coupled oscillations are completely delocalized even in the case of differently sized interaction partners, the density polarization is mostly localized in the overlap region between the two molecules (Fig. 6b).

Previous applications of MBD
The MBD scheme has been used to identify chemical and thermodynamic trends across a broad range of materials and complexes, and has often succeeded even when other vdW-inclusive methods have failed (Stöhr et al. 2019;Xu et al. 2020).
On many occasions, MBD has been able to accurately predict the thermodynamic stabilities of polymorphs often exhibited by molecular crystals, which are typically governed by non-covalent interactions (Stöhr et al. 2019;Xu et al. 2020).The ability to correctly identify the energetic rankings of such vdW-bound systems has applications in numerous fields such as pharmacy and organic electronics (Bučar et al. 2015;Stöhr et al. 2019).One well-studied case is oxalic acid (Stöhr et al. 2019), for which many vdW-inclusive methods do not correctly predict the relative stabilities of its two polymorphs, whereas PBE+MBD agrees with experimental results (Reilly & Tkatchenko 2015).Another example is coumarin, where the inclusion of beyond-pairwise interactions significantly improves the predicted energy ranking of crystalline polymorphs compared to the pairwise TS method (Shtukenberg et al. 2017).Furthermore, the explicit accounting for many-body interactions elucidates the relative prevalence of a particular aspirin polymorph ('form 1') over the other in nature, whereas most other electronic structure methods (both with and without vdW corrections) predict both polymorphs to be energetically degenerate (Reilly & Tkatchenko 2014).
MBD has also been shown to accurately calculate interaction energies for supramolecular host-guest complexes (Stöhr et al. 2019).In particular, DMC calculations show that the binding energies of the C 70 -fullerene to [10]-and [11]-cycloparaphenylene are degenerate (within DMC uncertainty) (Hermann et al. 2017a).Only the explicit inclusion of many-body interactions alongside DFT correctly predicts this degeneracy, whereas pairwise or twoand three-body vdW models show a clear preference for the 10membered ring (Antony et al. 2015;Hermann et al. 2017a).
Other applications that the MBD approach has been used for include layered materials (Bučko et al. 2016) and molecular switches (Chen et al. 2019).In the former case, the heat of adsorption of toluene on graphene, as calculated using PBE+MBD, was found to be in good agreement with experiment (Bučko et al. 2016).Furthermore, the balance between Pauli repulsion, chemical binding, electrostatic interactions, and MBD forces was observed to be critical for achieving molecular switches between the two bistable configurations of gold hexamers on single-walled carbon nanotubes, while the many-body effects were controlled by the former two interactions (Chen et al. 2019).
Many-body dispersion effects have also been shown to be important for the description of hybrid organic-inorganic in-terfaces where MBD clearly outperforms additive pairwise dispersion schemes.For example, MBD was shown to improve the adsorption energetics and geometries of adsorbates such as xenon, graphene, 3,4,9,10-perylene-tetracarboxylic dianhydride (PTCDA), and 7,7,8,8-tetracyanoquinodimethane (TCNQ) on silver surfaces, as they reduce the overbinding typically found in pairwise additive dispersion-correction schemes (Blowey et al. 2020;Maurer et al. 2015).Recently, PBE+MBD-NL provided an accurate prediction of the energetic ranking and the adsorption height of different adsorption phases of TCNQ on Ag(100) when compared to X-ray standing wave measurements (Sohail et al. 2023).MBD also permits changes in the anisotropic polarizability tensor, the description of adsorbate vibrations and adsorbate-surface interaction screening to be captured (Maurer et al. 2015).Furthermore, MBD interactions were shown to be a key factor behind the stabilization of a PTCDA molecule being brought to a standing configuration on Ag(111) surfaces using a scanning probe microscopy tip (Knol et al. 2021).At a temperature of 5 K, the MBD-calculated lifetime of the standing configuration was calculated to be 1.7 × 10 5 s, which was consistent with experimental data; however, the vdW surf -calculated lifetime was evaluated to be only 0.1 s (Knol et al. 2021).In addition, the MBD-calculated energy barrier heights were around 170% larger than the vdW surf -calculated barrier heights (Knol et al. 2021).This is because the screening of long-range vdW interactions via non-additive many-body interactions reduces the moleculesurface attraction, which ultimately stabilizes the molecule in the upright geometry due to higher-order contributions to the vdW energy contained in MBD that counteract the pairwise contributions within vdW surf .

Conclusions
libMBD is a mature, efficient, yet flexible library that offers all the functionality one might need to model vdW interactions in molecules and materials with MBD.It can be used to enhance any electronic structure code that can interface with Fortran, C, or Python external libraries with a powerful implementation of the MBD and TS methods, provided the said code can calculate the Hirshfeld volumes or evaluate the VV polarizability functional (for MBD-NL).With the recent advent of machine learning models for Hirshfeld volumes and atomic polarizabilities in general, it is even possible to calculate MBD energies with libMBD without running a single electronic structure calculation.libMBD is hosted and developed on GitHub under the Mozilla Public License 2.0, is distributed with the ESL Bundle (Oliveira et al. 2020), and can be easily installed from Conda-Forge and PyPI.We hope that libMBD either standalone or embedded in a growing number of third-party programs will equip researchers with a new computational tool to study molecules and materials.

Figure 3 |
Figure 3 | Strong scaling.(a) 5×5×5 urethane supercell used for benchmarking.(b) libMBD parallelizes well to hundreds of CPU cores, and is mostly limited only by the parallelization of the matrix operations in external libraries.The current simple implementation of the SCS forces uses an unnecessary amount of inter-process communication, resulting in a relatively early plateau for finite systems.This will be addressed in future work.

Figure 4 |
Figure 4 | Numerical convergence.Where possible, libMBD offers sane defaults for parameters that ensure numerical convergence.All errors are relative.Dotted vertical line denotes the default value of a parameter.(a) The MBD energy converges exponentially with respect to the frequency grid (Eq.16).(b) Polynomial convergence with respect to the size of the -point grid (Eq.8).(c) Exponential convergence with respect to the real-and reciprocal-space cutoffs,  c ,  c , in the Ewald summation (Eq.10).(d) Optimal range separation in Ewald summation to minimize computational cost.

Figure 6 |
Figure 6 | MBD wavefunction.libMBD gives access not only to energies, but also to other expectation values of the MBD wavefunction.Here calculations on the π-π-stacked circumcoronene⋯adenine dimer (C3A) of the L7 dataset.(a) One of the low-energy zero-point oscillations of the MBD Hamiltonian (Eq.6), displayed as instantaneous dipole displacements of the electronic clouds around each atom.(b) Polarization of the electron density due to the vdW interaction (Eq.22), with density accumulation and depletion colored as orange and blue, respectively.
), while the long-range part is used in the MBD Hamiltonian,