# How does MBD (MBD-NL) differ from DFT-TS?

Is anyone familiar with Tkatchenko's new MBD-NL method? How does it differ from DFT-TS? I couldn't understand everything from the paper.

• +1, but the D3 dispersion correction is completely different from the TS one, so by understanding how MBD-NL differs from TS you are automatically also getting the answer to how it differs from D3. This question is related: mattermodeling.stackexchange.com/q/63/5. – Nike Dattani Sep 21 '20 at 21:24

The MBD-NL method is an improvement of the 'traditional' MBD (sometimes called MBD@rsSCS) approach, which is quite a bit different from the TS method you've asked. So in order to really appreciate what the MBD-NL method entails, we first need to clarify what even are MBD and TS methods. I will try to point out the main ideas, however to get a full understanding, readong the papers from the group would be beneficial. I understand you've already tried to do that, and it is also my experience that the papers from the groups are relatively unclear. I've this one to be relatively instructive about MBD. It could also be helpful to look at Ph.D. theses from the group as they might be more instructive, they are available here.

## The TS approach

As you are probably well aware, the idea between all of the methods to be discussed is to add van der Waals correlation energies to DFT calculations. Most DFT functionals are simply incapable of modeling van der Waals effect, because the functionals only depend on local quantities, while van der Waals enegies need nonlocal electron correlation description (think: two-point correlation functions). One way to 'fix' this issue is to calculate the dispersion energy correction coarse-grained on atoms. These atomic approaches (such as the DFT-D3 already mentioned in one of the comments) usually all lead to a similar equation

$$E_{\mathrm{disp}} = - \frac12 \sum_{p \neq q} \frac{C_6^{pq}}{R_6^{pq}},$$

that is, you can get part of the missing correlation energy if you calculate pairwise interaction between the atoms using the $$C_6$$ coefficient (you might look at this question for more information here).

The main question is then: what is the $$C_6$$ coefficient of an atom in a molecule? Up to this point, most pairwise atomic methods agree, and this is the point where they start to diverge.The Tkatchenko-Scheffler (TS) method starts from the idea that each atom can be considered a Drude oscillator (that is, a harmonic oscillator with a charge, representing the fact that the potential energy is supposed to have a minimum at some distance from the nucleus after all, so a Taylor expansion there would yield a harmonic oscillator). This idea is relatively useful, and ideed the equation for the dispersion coefficient can also be derived from interacting harmonic oscillators (see the paper I've linked in the beginning). In order to get the $$C_6$$ coefficient of a given interaction, the so-called Casimir-Polder formula is then used

$$C_6^{pq} = \frac{3}{\pi} \int_0^{\infty} \mathrm{d} \omega \alpha_p(i \omega) \alpha_q(i \omega)$$

So then if we have the (dipole) polarizabilities of the atoms in your system, we can calculate everything. The main insight of the TS method is the idea to calculate the polarizability of atoms in the molecules (or solid state) by re-scaling the free atomic polarizability by calculating the relative volume of the atom in the molecule to the free atomic volume

$$\alpha = \frac{V}{V_{\mathrm{free}}} \alpha_{\mathrm{free}}$$

So then all you have to do is calculate the atomic volume in your system in question, which you can done however you want to do. Usually it's done with the so-called Hirshfeld partitioning.

So, in conclusion, the TS method does the following things:

1. Runs a 'normal' DFT calculation to get the electron density $$n(r)$$.
2. This electron density is then partitioned to atoms using for example a Hirshfeld scheme.
3. From the partitioned electron density, you can define a 'bound' atomic volume, which, together with the known free atomic parameters, gives you the polarizability of your atoms.
4. Modeling the atoms with harmonic oscillators to account for the frequency dependence, the Casimir-Polder integral is carried out.
5. Having obtained the $$C_6^{pq}$$ coefficients for all interacting pairs, you get the TS correction to the van der Waals energy.

## The MBD approach

The TS approach has been a huge success in describing van der Waals energies in DFT, however it has some shortcomings: only pairwise interactions are accounted for; and when calculating the polarizability, only the local electron density is considered. Many-body dispersion was intended to fix these issues by considering non-pairwise interactions as well as long-range screening effects.

The main idea is the following: you first do a TS calculation, which then basically gives you parameters for the Drude oscillator model of your atoms. Now that you've done that, you don't have to worry about the exact chemical nature of your system, you only have a bunch of interacting oscillating dipoles. These dipoles will then interact with each other, and this interaction is described by the dipole-dipole coupling tensor $$T_{ab}(r)$$. For more mathematical details, I suggest you look at page 8 of this paper. In practice, the self-consistent screening is done by solving the following Dyson-like equation

$$\alpha_{a,\mathrm{MBD}}(i \omega) = \alpha_{a,\mathrm{TS}}(i \omega) - \alpha_{a,\mathrm{TS}}(i \omega) T_{ab} \alpha_{b,\mathrm{MBD}}(i \omega)$$

If we solve this equation self-consistently, we will end up with dipole polarizabilities that contain the screening from all other oscillators in your system. To get an energy, we then give a Hamiltonian, in quantum mechanical sense, for the whole system $$H = -\frac12 \sum_n \nabla^2_n + \frac12 \sum_n \omega_n \mu_n^2 + \sum_{p \neq q} \omega_p \omega_q \mu_p T \mu_q,$$

where we take $$\mu$$ to be the massless normal coordinates we usually define for an ensemble of oscillators. The last term then contains a $$T$$ interaction tensor, which describes the interaction between Drude oscillators - the exact details are found in the paper linked above.

For this Hamiltonian, if we find the energies $$\lambda_m$$, then the correlation energy correction will be given as

$$E_{\mathrm{corr}} = \frac12 \sum_m^{3N} \sqrt{\lambda_m}-\frac32 \sum_m^N \omega_m$$

This is the main idea of MBD, and this algorithm outlined here is sometimes called MBD@rsSCS, if the self-consistent screening is done in real space.

## MBD-NL

So now we understand that TS puts non-interacting dipoles on the atomic position, and MBD couples these dipoles. What is MBD-NL then? In fact, it is just a small modification of the MBD scheme outlined above. As I've said, in MBD (or rather, already in TS), you need to partition the electron density to atoms in order to define atomic properties. Originally, Hirshfeld partitioning was used, which is convenient to use, and it is a mathematically nice operation, however it is not build on any real physical picture. Because of this, if you want to reproduce reference data with the TS method, you don't really get good result with the Hirshfeld scheme (look at Fig. 1. here).

Realizing these shortcomings, the MBD-NL approach set up some physical checks in the algorithm, namely

1. Instead of partitioning the electron density with a Hirshfeld approach, the polarizability itself is being partitioned, based on the VV10 (Vydrov and van Voorhis) formula
2. Some technical parameters used in the MBD implementation were re-designed, for example the cutoff functions were fit to more physical systems; and the van der Waals radii were re-evaluated based on new theoretical progress.

However, these modifications do not affect the main MBD idea, it just makes parametrizing the Hamiltonian slightly better. So to answer your question, the methodological difference between TS and MBD-NL is approximately is the same as between TS and MBD, namely the fact that MBD contains long-range screening as well as many-body interaction terms in the correction.

## Technical notes

If you are interested in trying MBD, it is already implemented in many different DFT codes (VASP,AIMS,Q-Chem,Quantum Espresso just to name a few), but a direct Python and Fortran interface is also available in a library called LibMBD. Sadly the documentation of this code is also a bit lackluster, however based on some insider information I've heard there's a technical paper being currently written on this library which will (hopefully) explain why and how it works, as well I think there's also a more practical paper being worked on that would show how to actually use it in real DFT calculations. However, I am not too sure about these information, so take it with a grain of salt.