A variational approach to London dispersion interactions without density distortion

TL;DR

This paper presents a variational wavefunction method that accurately models London dispersion interactions between neutral systems without altering their individual densities, enabling precise calculation of dispersion coefficients.

Contribution

It introduces a novel variational approach that captures dispersion forces without density distortion, providing explicit formulas for dispersion coefficients based on ground-state pair densities.

Findings

Accurately computes dispersion coefficients up to C10 for hydrogen atoms.

Achieves errors below 0.2% in simple cases.

Provides a theoretical framework for new approximation methods.

Abstract

We introduce a class of variational wavefunctions that capture the long-range interaction between neutral systems (atoms and molecules) without changing the diagonal of the density matrix of each monomer. The corresponding energy optimization yields explicit expressions for the dispersion coefficients in terms of the ground-state pair densities of the isolated systems, providing a clean theoretical framework to build new approximations in several contexts. As the individual monomer densities are kept fixed, we can also unambiguously assess the effect of the density distortion on van der Waals interactions: for example, we obtain virtually exact dispersion coefficients between two hydrogen atoms up to $C_{10}$, and relative errors below $0.2\%$ in other simple cases.

Tables

Table 1. Table 1: Dispersion coefficients C 6 subscript 𝐶 6 C_{6} , C 8 subscript 𝐶 8 C_{8} and C 10 subscript 𝐶 10 C_{10} between two H atoms computed with 30 terms in the variatonal expression of Eqs. ( 1 )-( 15 ), compared to the reference data of Thakkar Thakkar ( 1988 ) and Masili, et al. Masili and Gentil ( 2008 ) . See also supplementary material sup .

	Thakkar Thakkar (1988)	Masili et al.Masili and Gentil (2008)	This work (30 terms)
$C_6$	$6.4990267054058405\times {10}^0$	$6.4990267054058393\times {10}^0$	$6.4990267054058\mathrm{𝟑𝟗𝟑}\times {10}^0$
$C_8$	$1.2439908358362235\times {10}^2$	n.a.	$1.243990835836223\mathrm{𝟒}\times {10}^2$
$C_{10}$	$3.2858284149674217\times {10}^3$	n.a.	$3.2858284149674217\times {10}^3$

Equations

$${\Psi }_R({\mathbf{r}}_1,{\mathbf{r}}_2)={\Psi }_0^A({\mathbf{r}}_1){\Psi }_0^B({\mathbf{r}}_2)\sqrt{1+J_R({\mathbf{r}}_1,{\mathbf{r}}_2)},$$

$$\int d{\mathbf{r}}_1|{\Psi }_0^A({\mathbf{r}}_1){|}^2{J}_R({\mathbf{r}}_1,{\mathbf{r}}_2)=0\forall {\mathbf{r}}_2$$

$$\int d{\mathbf{r}}_2|{\Psi }_0^B({\mathbf{r}}_2){|}^2{J}_R({\mathbf{r}}_1,{\mathbf{r}}_2)=0\forall {\mathbf{r}}_1,$$

$$J_R({\mathbf{r}}_1,{\mathbf{r}}_2)=\sum _{ij}{c}_{ij,R}{b}_i^A({\mathbf{r}}_1)b_j^B({\mathbf{r}}_2),$$

$$b_i^A(\mathbf{r})=f_i^A(\mathbf{r})-\int \frac{{\rho }_0^A(\mathbf{r})}{N_A}{f}_i^A(\mathbf{r})d\mathbf{r}$$

$$\Delta T(\{c_{ij,R}\})=⟨{\Psi }_R|\hat{T}|{\Psi }_R⟩-T_0^A-T_0^B,$$

$$W_{\mathrm{int},0}(R)=⟨{\Psi }_0^A{\Psi }_0^B|{\hat{H}}_{\mathrm{int}}|{\Psi }_0^A{\Psi }_0^B⟩,$$

$$W_{\mathrm{int},\mathrm{c}}(\{c_{ij,R}\},R)=\int {\hat{H}}_{\mathrm{int}}|{\Psi }_0^A{|}^2|{\Psi }_0^B{|}^2{J}_{R}d{\mathbf{r}}_{1}d{\mathbf{r}}_2,$$

$$E_{\mathrm{int}}(R)=\Delta T(\{c_{ij,R}\})+W_{\mathrm{int},0}(R)+W_{\mathrm{int},\mathrm{c}}(\{c_{ij,R}\},R).$$

$${\hat{H}}_{\mathrm{int}}=\frac{{\hat{H}}_{dd}}{R^3}+\frac{{\hat{H}}_{dq}+{\hat{H}}_{qd}}{R^4}+\frac{{\hat{H}}_{qq}}{R^5}+...,$$

$$C_6=\frac{1}{2}{\mathbf{w}}^T{L}^{-1}\mathbf{w},$$

$$w_{ij}=\sum _e{h}_e\int e|{\Psi }_0^A(\mathbf{r}){|}^2{b}_i^A(\mathbf{r})d\mathbf{r}\int e|{\Psi }_0^B(\mathbf{r}){|}^2{b}_j^B(\mathbf{r})d\mathbf{r},$$

$$L_{ij,kl}=\frac{1}{4}\left({\tau }_{ik}^A{S}_{jl}^B+S_{ik}^A{\tau }_{jl}^B\right),$$

$${\tau }_{ij}^A=\int |{\Psi }_0^A(\mathbf{r}){|}^2\nabla b_i^A(\mathbf{r})\cdot \nabla b_j^A(\mathbf{r})d\mathbf{r},$$

$$S_{ij}^A=\int |{\Psi }_0^A(\mathbf{r}){|}^2{b}_i^A(\mathbf{r})b_j^A(\mathbf{r})d\mathbf{r},$$

$${\Psi }_R({\underline{\mathbf{x}}}_A,{\underline{\mathbf{x}}}_B)={\Psi }_0^A({\underline{\mathbf{x}}}_A){\Psi }_0^B({\underline{\mathbf{x}}}_B)\sqrt{1+\sum _{i\in A,j\in B}{J}_R({\mathbf{r}}_i,{\mathbf{r}}_j)},$$

$$\int |{\Psi }_R({\underline{\mathbf{x}}}_A,{\underline{\mathbf{x}}}_B){|}^{2}d{\mathbf{x}}_A=|{\Psi }_0^B({\underline{\mathbf{x}}}_B){|}^2\forall {\mathbf{x}}_B$$

$$\int |{\Psi }_R({\underline{\mathbf{x}}}_A,{\underline{\mathbf{x}}}_B){|}^{2}d{\mathbf{x}}_B=|{\Psi }_0^A({\underline{\mathbf{x}}}_A){|}^2\forall {\mathbf{x}}_A,$$

$$\int {\rho }_0^A({\mathbf{r}}_{i_A})J_R({\mathbf{r}}_{i_A},{\mathbf{r}}_{j_B})\mathrm{d}{\mathbf{r}}_{i_A}=0\forall {\mathbf{r}}_{j_B}$$

$$\int {\rho }_0^B({\mathbf{r}}_{j_B})J_R({\mathbf{r}}_{i_A},{\mathbf{r}}_{j_B})\mathrm{d}{\mathbf{r}}_{j_B}=0\forall {\mathbf{r}}_{i_A},$$

$$w_{ij}=\sum _{e=x,y,z}{h}_e(d_i^A+D_i^A)(d_j^B+D_j^B),$$

$$L_{ij,kl}=\frac{1}{4}\left({\tau }_{ik}^A(S_{jl}^B+P_{jl}^B)+(S_{ik}^A+P_{ik}^A){\tau }_{jl}^B\right),$$

$$d_i^A$$

$$D_i^A$$

$$P_{ij}^A$$

$$C_n[P_0^A,P_0^B]=C_n[{\rho }_0^A,h_{xc}^A,{\rho }_0^B,h_{xc}^B]\to C_n[{\rho }_0^A,{\rho }_0^B].$$

$$\begin{split}&C_{6,ab}=\frac{1}{2}\int_{\Omega_{a}^{A}}\mathrm{d}\mathbf{r}_{1_{A}}\int_{\Omega_{b}^{B}}\mathrm{d}\mathbf{r}_{1_{B}}\hat{H}_{dd}(\mathbf{r}_{1_{A}},\mathbf{r}_{1_{B}})\Big{(}\rho_{0}^{A}(\mathbf{r}_{1_{A}})\rho_{0}^{B}(\mathbf{r}_{1_{B}})\\ &\times J(\mathbf{r}_{1_{A}},\mathbf{r}_{1_{B}})+\int\mathrm{d}\mathbf{r}_{2_{A}}P_{0}^{A}(\mathbf{r}_{1_{A}},\mathbf{r}_{2_{A}})\rho_{0}^{B}(\mathbf{r}_{1_{B}})J(\mathbf{r}_{2_{A}},\mathbf{r}_{1_{B}})\\ &+\int\mathrm{d}\mathbf{r}_{2_{B}}\rho_{0}^{A}(\mathbf{r}_{1_{A}})P_{0}^{B}(\mathbf{r}_{1_{B}},\mathbf{r}_{2_{B}})J(\mathbf{r}_{1_{A}},\mathbf{r}_{2_{B}})\\ &+\int\mathrm{d}\mathbf{r}_{2_{B}}\mathrm{d}\mathbf{r}_{2_{A}}P_{0}^{A}(\mathbf{r}_{1_{A}},\mathbf{r}_{2_{A}})P_{0}^{B}(\mathbf{r}_{1_{B}},\mathbf{r}_{2_{B}})J(\mathbf{r}_{2_{A}},\mathbf{r}_{2_{B}})\Big{)}\end{split}$$

Full text

Within the Born-Oppenheimer approximation, non-relativistic quantum mechanics of Coulomb systems (atoms, molecules) predicts that everything is attracted to everything, since it is always possible to find a relative orientation that lowers the total energy of two neutral molecules separated by a large distance $R$ by at least a term $-C_{6}R^{-6}$, with $C_{6}>0$ Lieb and Thirring (1986); Kaplan (2006). This universal attractive interaction between quantum Coulomb systems, often called London dispersion interaction, is very small compared to the energy scale of a typical chemical (covalent) bond. Yet, it is omnipresent in chemistry, biology, and physics, determining protein folding, the structure of DNA, the physics of layered materials, just to name a few examples. An accurate, efficient, and fully non-empirical treatment of these forces remains an open challenge in many respects, and it is the focus of several on-going research efforts Grimme (2011); Burns et al. (2011); Dobson and Gould (2012); Grimme et al. (2016); Vydrov and Voorhis (2010a); Ferri et al. (2015); Becke (2013).

The physical origin of London dispersion interactions is often discussed in terms of “coupling between instantaneous dipoles” Kaplan (2006); Stone (1996), which are actually correlations in the interfragment part of the pair density (the diagonal of the two-body reduced density matrix), as illustrated, e.g., in Gritsenko and Baerends (2006). The textbook treatment is based on perturbation theory Kaplan (2006); Stone (1996); Margenau et al. (1971), in which excitations on both monomers are coupled, which allows to rewrite the dispersion coefficients in terms of dynamical polarizabilities of the individual systems. In this Letter, inspired by the work of Lieb and Thirring (LT) Lieb and Thirring (1986), we attack the dispersion problem from a variational perspective, using only ground-state properties of the monomers. We construct a class of variational wavefunctions for the long-range interaction between neutral Coulomb systems that are explicitly forbidden to deform the diagonal of the full spatial density matrix (and thus the density) of each individual monomer. This constraint has several advantages: first of all, it is conceptually appealing, as dispersion becomes a monomer-monomer interaction accompanied by a change in kinetic energy only, since all the remaining intra monomer interactions are kept unchanged. Secondly, it highly simplifies the optimization procedure, which can be implemented in a very efficient way and leads to expressions for the dispersion coefficients in terms of the ground-state pair densities of the individual monomers, giving a neat framework to build new approximations. For example, in a density functional theory (DFT) setting, one can use approximations for the exchange-correlation holes of the fragments, providing justification and a route for improving models such as the exchange-hole dipole moment (XDM) Becke and Johnson (2007). It also provides a prescription to obtain atomic dispersion coefficients inside molecules and solids, which are a crucial ingredient in the construction of atomistic force fields and dispersion corrections to DFT calculations (DFT+D).

One could easily object that this construction cannot give the exact answer. In a widely quoted paper, Feynman Feynman (1939) wrote that van der Waals’ forces arise from charge distributions with higher concentration between the nuclei, with a permanent dipole moment of order $R^{-7}$ being induced on each atom, such that “each nucleus is attracted by the distorted charge distribution of its own electrons” Feynman (1939). Feynman’s conclusion was proven by using the electrostatic Helmann-Feynman theorem for both atoms and molecules Hunt (1990), and has been confirmed with various accurate calculations Allen and Tozer (2002); Strømsheim et al. (2011). Therefore we know that a small density distortion must be there. Nonetheless, it has been observed several times that very accurate dispersion coefficients can be obtained with methods that describe very poorly this density distortion, including cases in which the electrostratic Helmann-Feynman theorem gives zero force at order $R^{-7}$ Hirschfelder and Eliason (1967). The issue has been often debated in the literature, and there is the general feeling that the density distortion, which actually requires expensive calculations to be accurately captured, should have very little influence on dispersion energetics, which is our goal here. Moreover, our construction provides (if combined with accurate pair densities for the monomers), a rigorous upper bound to the best possible variational interaction energy when the density distortion is explicitly forbidden, allowing us to assess its importance in an unambiguous manner, which is per se also conceptually interesting.

The wavefunction for $N=2$ – For illustrative purposes, we discuss first the case of two systems $A$ and $B$, each with $N_{A}=N_{B}=1$ electron, generalizing it right after to the many-electron case. The two atoms are separated by a large distance $R$, and our wavefunction reads

[TABLE]

where ${\bf r}_{1}$ and ${\bf r}_{2}$ are electronic coordinates centered in nucleus $A$ and nucleus $B$, respectively, $\Psi_{0}^{A(B)}({\bf r})$ is the ground-state wavefunction of system $A$ (or $B$) alone, and $J_{R}$, with $|J_{R}({\bf r}_{1},{\bf r}_{2})|\ll 1$, preserves the unperturbed densities at all orders in $R^{-1}$ via the constraints

[TABLE]

which are imposed with an expansion in terms of one-electron functions $b_{i}^{A}({\bf r})$ and $b_{i}^{B}({\bf r})$ and variational parameters $c_{ij,R}$,

[TABLE]

where

[TABLE]

(and similarly for $b_{j}^{B}$), with $\rho_{0}^{A(B)}({\bf r})=|\Psi_{0}^{A(B)}({\bf r})|^{2}$. The $f_{i}^{A(B)}({\bf r})$’s can be chosen in different ways, with the choice influencing the convergence. Our wavefunction does not include antisymmetrization between the electron on $A$ and the one on $B$, which affects the asymptotic (large $R$) interaction energy only with terms that vanish exponentially with $R$. This wavefunction correlates the positions of the two electrons in the two different atoms without changing their one-electron densities.

We now consider the full hamiltonian of the problem, $\hat{H}=\hat{H}_{A}({\bf r}_{1})+\hat{H}_{B}({\bf r}_{2})+\hat{H}_{\rm int}({\bf r}_{1},{\bf r}_{2})$, where $\hat{H}_{A}$ and $\hat{H}_{B}$ are the hamiltonians of each isolated atom and $\hat{H}_{\rm int}$ contains all the interactions between the particles in $A$ (electron and nucleus) with those in $B$, and evaluate its expectation value on our wavefunction, focussing on the interaction energy $E_{\rm int}=E_{0}^{AB}-E_{0}^{A}-E_{0}^{B}$. Thanks to the density constraint of Eqs. (2)-(3), there is no change in the expectation value of the nuclear-electron potential, implying that the interaction energy is entirely determined by the change in kinetic energy $\hat{T}$ with respect to the one of the isolated atoms,

[TABLE]

plus the expectation value of $\hat{H}_{\rm int}$, which splits automatically into the sum of an electrostatic part $W_{\rm int,0}(R)$ and the correlation part $W_{\rm int,c}(\{c_{ij,R}\},R)$,

[TABLE]

As usual Kaplan (2006), we perform a multipolar expansion of $\hat{H}_{\rm int}$,

[TABLE]

where $d$ stands for dipole, $q$ for quadrupole, etc. The variational problem then takes a simplified form, as detailed in the supplementary material sup , with each order $R^{-2n}$ in the energy having its optimal $c_{ij,R}=c_{ij}\,R^{-n}$. The physics described by this wavefunction can be understood by looking at the leading order: our wavefunction gives a $\Delta\,T(\{c_{ij}\})$ which is positive (as it costs kinetic energy to correlate the electrons) and quadratic in the variational parameters $c_{ij}$, while $W_{\rm int,c}(\{c_{ij}\},R)$ is linear in the $c_{ij}$, negative (as correlation occurs exactly to lower the interaction term), and has a prefactor $R^{-3}$ from (10). This guarantees that there are always optimal variational parameters $c_{ij}$, all proportional to $R^{-3}$, which minimize $E_{\rm int}$, whose minimum becomes the familiar $-C_{6}/R^{6}$. This structure repeats at each order in $1/R$. In LT Lieb and Thirring (1986) and subsequent works in a similar spirit,Anapolitanos and Sigal (2017); Anapolitanos and Lewin (2018) the energy cost to correlate the electrons, quadratic in the variational parameters, comes from the full hamiltonian of the monomers, while with our construction we force it to be of kinetic energy origin only.

Solving the variational equations leads to explicit expressions for all the even dispersion coefficients. For example, $C_{6}$ is given by

[TABLE]

where the vector $\mathbf{w}$ has contracted indices $I=ij$,

[TABLE]

with $e=x,y,z$, and $h_{e}=(1,1,-2)$ when the bond is along the $z$-axis. The matrix $L$ has two contracted indices $I=ij$ and $K=kl$,

[TABLE]

with

[TABLE]

and similarly for $B$. Analogous expressions hold for all the even dispersion coefficients (see sup ).

We have started with a very simple choice for the $f_{i}^{A(B)}({\bf r})$’s (just powers of the distance $r$ from the nucleus times spherical harmonics sup ), finding that the convergence is already fast, as shown in Fig. 1 for the case of $C_{6}$ for two H atoms. The values of $C_{6}$, $C_{8}$ and $C_{10}$ for the same case are reported in Table 1, where they are compared with reference data from Ref. 22. In the supplementary material sup , we also report even second-order coefficients up to $C_{30}$, which also agree within numerical accuracy with those of Ref. 22. From Table 1 we see that our wavefunction yields essentially the exact lowering in energy due to dispersion without any distortion of the density of each individual atom.

The many-electron case – We now turn to the many-electron case, by considering two systems with $N_{A}$ and $N_{B}$ electrons, for which our wavefunction reads

[TABLE]

where ${\bf x}_{A(B)}$ indicates all the electronic coordinates (including spins) of the electrons in system $A$ (or $B$), and, again, $\Psi_{0}^{A(B)}$ is the ground state wavefunction of system $A$ (or $B$) alone. We choose $J_{R}$ such that the diagonal of the full many-body density matrix of each system is not allowed to deviate from its ground-state expression:

[TABLE]

which is enforced by choosing for $J_{R}$ the same form as Eq. (4), i.e., by imposing that

[TABLE]

with $\rho_{0}^{A(B)}$ the ground-state one-electron densities of the two systems. Again, it is only the interfragment pair density and the off-diagonal elements of the first-order reduced density matrix that are allowed to change to lower the energy, implying that also in the many-electron case the interaction energy is given by Eq. (9), since the expectation value of the electron-electron interaction inside each monomer is kept unchanged by virtue of Eqs. (17)-(18). The variational minimization of the interaction energy is similar to the $N=2$ case. For example, $C_{6}$ takes the same form as Eq. (11), where now $\mathbf{w}$ and $L$ read

[TABLE]

where (with the same expressions for $B$)

[TABLE]

and where $\tau_{ij}^{A(B)}$ and $S_{ij}^{A(B)}$ are defined in Eqs. (14)-(15). We see that the dispersion coefficients $C_{n}$ for the many-electron case become explicit functionals of the spin-summed ground-state pair densities $P_{0}^{A}(\mathbf{r}_{1_{A}},\mathbf{r}_{2_{A}})$ and $P_{0}^{B}(\mathbf{r}_{1_{B}},\mathbf{r}_{2_{B}})$ of the two monomers, $C_{n}=C_{n}[P_{0}^{A},P_{0}^{B}]$. The variational solution for the optimal $c_{ij}$ can also be turned into a Sylvester equation, which can be solved in an efficient way Bartels and Stewart (1972); sup .

Being derived from a wavefunction, the expression for the interaction energy is variational (i.e., it provides a lower bound to the exact $C_{6}$) as long as we use exact (or very accurate) pair densities for the monomers. As an example, with the same simple choice for the $f_{i}^{A(B)}({\bf r})$ used for two H atoms, we have computed $C_{6}$ for He-He using the accurate variational wavefunction for the He atom of Freund, Huxtable and Morgan Freund et al. (1984) (see also Ref. 26). The convergence is again fast, very similar to the one of Fig. 1, yielding $C_{6}=1.458440$ with an error of 0.17% with respect to the accurate value 1.460978 of Ref. 27. For the He-H case we have a similar error of 0.15% with respect to Ref. 27. Whether these residual errors are due to the reduced variational freedom of our wavefunction or to inaccuracies in the multiple moment integrals of the correlated He pair density of Ref. 25 is for now an open question, but we see that the results can be rather accurate, at a computational cost which is determined by the individual monomer calculations. This opens many possibilites, as our expressions can be coupled with any wavefunction method for the monomers, provided it gives access to the second-order reduced density matrices and their integrals with the $f_{i}^{A(B)}({\bf r})$, whose optimal choice needs to be assessed in a systematic way.

Lower-level approximations – We can of course use pair densities from lower level theories, which opens another whole realm of possibilities, bearing in mind that the expressions will be no longer variational for the interaction energy. For example, in Fig. 2 we show the result for $C_{6}$ for Ne-Ne using the Hartree-Fock (HF) pair density, including the basis set dependence. We see that the error is around 5%, as HF is a good approximation for the Ne atom, and that we are below the exact interaction energy. It is clear that a detailed study of the basis set dependence is also needed, which should be coupled to an optimal choice for the $f_{i}^{A(B)}({\bf r})$. The use of Kohn-Sham orbitals from different approximate functionals to determine the exchange pair densities is also worth to be investigated.

Within DFT, one could make even simpler approximations for the exchange-correlation holes $h_{xc}^{A(B)}$ of the monomers, transforming the dispersion coefficients into density functionals,

[TABLE]

For example, the Becke-Roussel Becke and Roussel (1989) exchange-hole model could be used to estimate higher moments of the exchange pair density Gori-Giorgi et al. (2009), in a spirit similar to the XDM model Becke and Johnson (2007). The main difference is that now the polarizability of the monomers is no longer needed, as the expression for the dispersion coefficients only requires the exchange-correlation holes.

An important ingredient in atomistic force fields models and DFT+D are the dispersion coefficients $C_{n,ab}$ between pairs of atoms $a\,b$, with $a$ in system $A$ and $b$ in system $B$, which are usually computed from atomic polarizabilities. Our variational wavefunction can provide a framework to define and compute them. One should obtain atomic volumes $\Omega_{a}^{A}$ and $\Omega_{b}^{B}$ using one of the various definitions (e.g. Hirshfeld partitioning Hirshfeld (1977)), and then obtain the pair dispersion coefficients using, order by order, the energy density coming from our wavefunction with the $c_{ij}$ optimized for the $AB$ interaction, restricting the integration on the atomic volumes. For example for $C_{6}$ (with $J$ below being the optimized factor at leading order),

[TABLE]

and similarly for the higher order coefficients.

Conclusions – In conclusion, the class of variational wavefunctions we have introduced here provides a neat theoretical framework to build new approximations to tackle dispersion interactions, which are reduced to a competition between kinetic energy and monomer-monomer interaction only. Although this class of variational wavefunctions violates the virial theorem (similarly to second-order perturbation theory), and although we know that for the exact wavefunction the density must be distorted, the results can be very accurate (or even exact), with a computational efficiency that is promising for the many-electron case. Moreover, we know what is being approximated, so that one could also devise a scheme in which the monomer (pair) densities are relaxed in a second step, or perturbatively, if this is needed. The realm of possibilities open is vast, ranging from the use of various wavefunction methods to compute the monomer pair densities or approximations thereof, to the definition and calculation of atomic-pair dispersion coefficients. The influence of the level of theory on the calculation of the monomer pair densities, the basis set, and different choices for the $f_{i}^{A(B)}({\bf r})$ need all to be investigated. We have also shown that exact dispersion coefficients up to $C_{10}$ between two hydrogen atoms can be obtained without introducing any density distortion, which is interesting by itself. The framework can be also generalised to many monomers.

Acknowledgements – Financial support from European Research Council under H2020/ERC Consolidator Grant corr-DFT (Grant Number 648932) and the Netherlands Organisation for Scientific Research under Vici grant 724.017.001 is acknowledged.