A fractionally ionic approach to polarizability and van der Waals many-body dispersion calculations

TL;DR

This paper introduces a novel method that explicitly incorporates fractional ionic contributions to improve the accuracy of polarizability and van der Waals many-body dispersion calculations, especially in ionic and nanotechnological systems.

Contribution

It presents a fractionally ionic approach that enhances existing atom-wise many-body van der Waals methods with minimal additional computational cost.

Findings

Significantly improves polarizability predictions in ionic solids by over 65%.

Achieves comparable accuracy to existing methods for non-ionic systems.

Provides better modeling of two-dimensional transition metal dichalcogenides and H₂ interactions with modified coronenes.

Abstract

By explicitly including fractionally ionic contributions to the polarizability of a many-component system we are able to significantly improve on previous atom-wise many-body van der Waals approaches with essentially no extra numerical cost. For non-ionic systems our method is comparable in accuracy to existing approaches. However, it offers substantial improvements in ionic solids, e.g. producing better polarizabilities by over 65% in some cases. It has particular benefits for two-dimensional transition metal dichalcogenides, and interactions of H$_2$ with modified coronenes - ionic systems of nanotechnological interest. It thus offers an efficient improvement on existing approaches, valid for a wide range of systems.

Tables

Table 1. Table 1: Summary of simulation parameters used in this study. See Tab.S5 in the Supporting Information 35 for details on setting used in X23 set calculations.

System	$𝒌$-point mesh	$E_{\text{cut}}$ (eV)
S66$\times$8	$1\times 1\times 1$	1000
X40	$1\times 1\times 1$	1000
H2+coronenes	$1\times 1\times 1$	1000
NaCl	$16\times 16\times 16$	1000
MgO	$16\times 16\times 16$	1000
cryolite	$8\times 8\times 6$	1000
bilayers of dichalcogenides	$16\times 16\times 1$	1000
bulk dichalcogenides	$12\times 12\times 3$	1500
fluorographite	$20\times 20\times 8$	1000

Table 2. Table 2: Mean absolute error (MAE) and mean absolute relative error (MARE) for interaction energies computed for the S66 × \times 8  67 and the X40  68 benchmark sets computed using the MBD@rsSCS and MBD@rsSCS/FI methods. The CCSD(T) results from Refs.  67 , 68 have been used as a reference. Statistics for subsets dimers differing in the nature of dominating interaction is also listed (the numbers in parentheses indicate the number of dimers in the given subset).

set	subset	MBD@rsSCS		MBD@rsSCS/FI
		MAE (meV)	MARE ($\%$)	MAE (meV)	MARE ($\%$)
S66$\times$8	all (66$\times$8)	13.6	9.7	12.3	9.0
	H-bond (23$\times$8)	19.1	8.3	17.8	7.8
	dispersion (23$\times$8)	17.5	6.4	15.4	5.7
	other (20$\times$8)	8.9	6.6	8.1	6.2
X40	all (40)	15.5	10.0	13.9	9.0
	dispersion (4)	3.9	16.4	3.6	15.9
	induction (4)	3.0	8.5	2.5	7.3
	dipole-dipole (2)	8.8	13.6	8.9	13.7
	stacking (2)	38.4	16.8	30.2	13.3
	halogen-bond (14)	5.5	5.0	4.7	4.3
	halogen-$\pi$ (4)	26.4	22.5	23.3	19.3
	H-bond (10)	31.5	8.0	29.3	7.5

Table 3. Table 3: Interaction energy per a H 2 molecule (kJ/mol) for set of complexes of carbonaceous and substituted coronenes. The CCSD(T) results are from Kocman et al.  69

	CCSD(T)	MBD@rsSCS	MBD@rsSCS/FI
Coronene…2H2	$-$4.7	$-$4.5	$-$4.6
CoroB2Li2…2H2	$-$14.3	crash	$-$15.9
CoroB2…2H2	$-$4.9	$-$3.7	$-$3.7
CoroB2Li2…H2ss	$-$10.4	crash	$-$12.9
CoroB2Li2…H2os	$-$5.0	crash	$-$5.6
C-coro…2H2	$-$5.5	$-$4.9	$-$5.0

Table 4. Table 4: Statistics for interaction energy and cell volume for molecular crystals from the benchmark set X23 computed using the MBD@rsSCS and MBD@rsSCS/FI methods. The reference values are from Ref.  71 . Details on energetics and relaxed structrures of individual systems are given in the Supporting Information 35 .

	MBD@rsSCS	MBD@rsSCS/FI
MAE(E) (kJ/mol)	5.7	5.2
MARE(E) ($\%$)	7.4	6.9
MAE(V) (Å3/fu)	1.9	2.0
MARE(V) ($\%$)	1.9	2.0

Table 5. Table 5: Static lattice constants and bulk moduli for crystalline NaCl, MgO and LiF. The zero-temperature experimental values corrected for zero-point phonon effects are from Ref.  75 .

	NaCl		MgO		LiF
	$a$ (Å)	$B_0$ (GPa)	$a$ (Å)	$B_0$ (GPa)	$a$ (Å)	$B_0$ (GPa)
Exp.	5.57	28	4.19	170	3.97	76
PBE	5.55	25	4.25	152	4.07	67
MBD@rsSCS	crash		crash		crash
MBD@rsSCS/FI	5.62	26	4.20	166	4.05	68

Table 6. Table 6: Interlayer distances d 0 subscript 𝑑 0 d_{0} , and interlayer binding energies E b subscript 𝐸 𝑏 E_{b} for AA’ and AB stacking variants of bilayers of transition-metal dichalcogenides from RPA and dispersion corrected DFT calculations. The RPA and pairwise TS results are from He et al. (Ref.  78 ) and Bučko et al. (Ref.  73 ), respectively.

	RPA		TS		MBD@rsSCS		MBD@rsSCS/FI
System(stacking)	$d_0$ (Å)	$E_b$ (meV/fu)	$d_0$ (Å)	$E_b$ (meV/fu)	$d_0$ (Å)	$E_b$ (meV/fu)	$d_0$ (Å)	$E_b$ (meV/fu)
MoS2(AA’)	6.27	81.2	6.08	149.2	crash		6.15	90.3
MoS2(AB)	6.17	76.7	6.08	140.4	crash		6.22	80.6
MoSe2(AA’)	6.48	88.4	6.44	148.7	crash		6.42	110.2
MoSe2(AB)	6.47	85.1	6.55	135.6	crash		6.61	88.7
WS2(AA’)	6.24	82.9	6.18	134.2	crash		6.26	76.7
WS2(AB)	6.24	74.3	6.28	124.7	crash		6.37	67.5
WSe2(AA’)	6.50	89.9	6.52	135.8	crash		6.54	93.1
WSe2(AB)	6.54	84.2	6.73	123.5	crash		6.75	76.0
MAE	-	-	0.09	53.7	-		0.10	7.8

Table 7. Table 7: Cell geometry and interlayer binding energy E B subscript 𝐸 𝐵 E_{B} for selected dichalcogenides. Experimental structures and the energies computed at the RPA level  79 are used as a reference. The pairwise TS results are from Ref.  73 .

Compounds	Method	$a$ (Å)	$c$ (Å)	$E_B$ (meV/Å2)
MoS2	Expt./RPA	3.16	12.29	20.5
	TS		12.09	37.9
	MBD@rsSCS	-	crash
	MBD@rsSCS/FI	-	12.28	20.7
MoSe2	Expt./RPA	3.29	12.90	19.6
	TS		12.78	34.0
	MBD@rsSCS	-	crash
	MBD@rsSCS/FI	-	12.78	20.6
WS2	exp./RPA	3.15	12.32	20.2
	TS		12.22	34.2
	MBD@rsSCS	-	crash
	MBD@rsSCS/FI	-	12.46	18.3
WSe2	Expt./RPA	3.28	12.96	20.0
	TS		12.97	31.3
	MBD@rsSCS	-	crash
	MBD@rsSCS/FI	-	13.01	19.8

Equations

$${\alpha }_p^{\text{AIM}}(i\omega ;N_p)={\alpha }_p^{\text{FI}}(i\omega ;N_p)\frac{V_p^{\text{eff}}(N_p)}{V_p^{\text{FI}}(N_p)}.$$

$$P_N=$$

$$\equiv$$

$$C_{6,AB}=\frac{3}{\pi }\int d\omega {\alpha }_{Z_A,N_A}(i\omega ){\alpha }_{Z_B,N_B}(i\omega )$$

$${\alpha }_{Z,N}(i\omega )=$$

$$n^{\prime }(\mathbf{r})=$$

$$n_N(\mathbf{r})=$$

$$n_N^{\prime }(\mathbf{r})=f{n}_{M+1}^{\prime }(\mathbf{r})+(1-f)n_M^{\prime }(\mathbf{r}).$$

$${\chi }_N(i\omega )=$$

$${\alpha }_{Z,N}(i\omega )=$$

$$\alpha (i\omega )=$$

$${\alpha }_{Z,M}(i\omega )=$$

$${\alpha }_{Z,N}(i\omega )=$$

$$+(1-f)\sum _{c=1,2}\frac{a_c(Z,M)}{{\omega }^2+{\Omega }_c(Z,M{)}^2}.$$

$$E_{\text{disp}}=$$

$${\tilde{\alpha }}_p(i\omega )={\alpha }_p(i\omega )\left(1-\sum _q{T}_{SR,pq}(i\omega ){\tilde{\alpha }}_q(i\omega )\right).$$

$$E_{\text{disp}}=$$

$${\tilde{x}}_n=$$

$$E_{\text{disp}}=$$

$$\equiv$$

$${\alpha }_m=\frac{3V}{4\pi }\left(\frac{{ϵ}_{\infty }-1}{{ϵ}_{\infty }+2}\right),$$

Full text

1 Introduction

Over the past decade there has been a resurgence of interest within the electronic structure community in the inclusion of van der Waals (vdW) dispersion forces in ab initio calculations. As a result, a plethora of new approaches has been developed1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 that allow van der Waals forces (more generally implemented as potentials) to be included alongside conventional density functional approximations13, 14, 15. Recent theoretical and experimental progress on “non-additivity” of van der Waals interactions16, 17, 18, 19, 20, 21, 22, 12, 23, 24, 25, 26, 27, 28, 29, 30, 31 has also seen a more fundamental shift in thinking about van der Waals forces and the role they play in nanomaterials.

This reinvigorated interest in vdW forces has both parallelled and been driven by increasing interest in layered “2D” materials which are held together almost exclusively by weak dispersion (vdW) forces. Without accurate models of vdW forces it is very difficult to predict even qualitatively correct properties of such layered systems

• including lattice parameters, binding energies and elastic moduli. Thus the inclusion of vdW effects is vital for the modeling of layered structures; and layered systems provide an important test of their accuracy32.

Some of the more popular approaches for including vdW forces are that of Tkatchenko and Scheffler10 (TS); related many-body dispersion (TS-MBD) schemes11, 20, 33, 34; and self-consistent screening11, 29 (TS-SCS) schemes. These approaches allow vdW forces to be included in calculations via a point-dipole model, with excellent numerical efficiency and generally good qualitative results. The methodology behind these approaches can be summarized as follows:

Pre-calculate $C_{6}$ coefficients and static dipole polarizabilities $\alpha_{0}$ for isolated atoms using high-level electronic structure theory. 2. 2.

Approximate dipole polarizability properties of a molecule or material as a superposition of polarizabilities of free-atoms centered on their nuclei. 3. 3.

Use the free-atom properties, appropriately rescaled by the effective volume of embedded atoms, to calculate a model frequency dependent (dynamic) polarizability (see Sec.S1 in the Supporting Information35). 4. 4.

Use the model polarizability to calculate the dispersion energy. This is done either by summing over pairwise contributions, as in the original TS theory10, or by solving a many-body screening equation that includes non-additive terms, as in many-body approaches11, 20, 33.

The TS, TS-MBD and TS-SCS schemes have shown great success in improving the qualitative and quantitative accuracy of ab initio calculations of many systems36, 37, 38, 39. This reflects their ability to capture Type-A (Dobson-A) and many-body Dobson-B effects, according to the classificiation of Dobson23. However, it was shown in Ref. 39 that the TS and TS-SCS schemes may perform poorly for ionic systems, due to their use of atomic rather than ionic volumes. As a solution, an alternative to stage 3 was proposed39, 40 that accounts for rescaled (via the ratio of the effective embedded volume to the free atom volume) densities of ions, rather than pure atoms. This was found to improve results in many ionic systems. However, while the improvements from the ionic rescaling go some way to improving results for ionic materials, it still fails to capture the full ionic physics of the systems, leading e.g. to poor dipole polarizabilities of ionic systems.

In this work we will go further than Ref. 40 and take into account the effect of both ionic charge and volume on the effective static polarizabilities $\alpha_{0}$ and $C_{6}$ coefficients of the system. More precisely, we will take into account the effect of the ionic charge on the frequency-dependent dipole polarizibilities, drawing from a recently published ionic polarizability dataset produced by two of the authors of the present work41. Our new method thus modifies stages 1, 3 and 4 of TS-derived schemes (vide supra) by including: a) ionic charge in the high-level calculations in stage 1, via a new ionic polarizability dataset41; b) ionic (rather than atomic) volume scaling of the polarizabilities in stage 3, and c) a more accurate model of frequency-dependent polarizabilities in stages 3 and 4.

We will first develop a theory of embedded fractional ions by: i) showing how a dataset of integer ions can be used to generate fractional ions; ii) introducing a “remapping” correction to the MBD equations that corrects for rare, but sizeable, unphysical interactions in some materials. We will finally show that the explicit inclusion of fractional charge in the point-dipole model substantially improves results in various systems, especially in ionic ones. Our new scheme is particularly important for ionic layered transition metal dichalcogenide (TMD) systems, like MoS2 where other schemes perform very poorly. In these systems the ionic scheme gives results on a par with sophisticated electron structure calculations using random-phase approximation (RPA)19, 42, 12, 32.

2 Theory

Our new scheme is closely related to the work of Tkatchenko and co-workers10, 11, 20, as improved by Bučko et al40. Notably, it employs the same iterative Hirshfeld partitioning and the same effective volume scaling as Ref. 40. However, it differs crucially through the polarizability model it employs, because it incorporates fractionally ionic effects and uses a significantly improved screening. The new recipe thus becomes:

Pre-calculate dipole polarizabilities not only for free atoms, but also for non-interacting free ions using high-level electronic structure theory. 2. 2.

Approximate the dipole polarizability properties of a molecule or a material via non-interacting free-ion polarizabilities centered on the nuclei. 3. 3.

Use the free-ion polarizabilities to obtain fractional ionic (FI) polarizabilities $\alpha_{p}^{\text{FI}}$, evaluated at the fractional electron number $N_{p}$, and rescaled111We tested more sophisticated models of rescaling from Ref. 43 but found it did not noticeably improve results. by the effective volume of the embedded charge, to calculate the effective fractional ionic polarizabilities (see Sec.S1 in the Supporting Information35 for further details):

[TABLE] 4. 4.

Use the effective, fractional ionic embedded polarizabilities distributed on the atomic sites to calculate the dispersion energy, either using the additive or many-body method with eigenvalue remapping (discussed later).

These items are treated in detail over the following sections. We first introduce the fractional ionic scheme [Sec. 2.1]; we then discuss the polarizability model used in the fractional ionic scheme [Sec. 2.2]; and we finally explain a “remapping” scheme that ensures that the energy and forces behave reasonably even for some of the systems where the original scheme fails [Sec. 2.3].

2.1 Properties of fractional ions

Before introducing the details of the new method, we should first discuss here what is meant by “fractional ions”. In an embedded system this concept is fairly well defined, at least within partitioning schemes such as Bader 44, Hirshfeld 45, or iterative Hirshfeld 46 atoms-in-molecules (AIM) approaches: the “fractional ion” is defined as a sub-system of a larger molecule or material with an integral nuclear charge $Z$ and a non-integer number of electrons $N$, whose overall charge $f=Z-N$ is thus non-integral in general. This is related to a representative free fractional ion with equal non-integer charge. The use of fractional ions raises interesting fundamental questions, whilst offering significant practical improvements. We thus dedicate the following part of this Section to discuss the nature and practical consequences of using such fractional ions.

Any isolated system with non-integer electron number must be treated via a quantum state ensemble composed of wave functions with different integer electron numbers. This applies at the fully interacting and density functional theory levels47. In an isolated atomic-like system, such as the interaction-free fractional ions that we will use, these states will all share the same form of nuclear potential $-Z/r$. Furthermore, at low temperatures, the ensemble of non-interacting free ions with non-integer $N$ electrons will be formed only from states with $\lfloor{N}\rfloor=M$ and $\lceil{N}\rceil=M+1$ electrons (i.e. the largest integer smaller and the smallest integer larger than $N$). As a consequence, many ground state properties of such systems obey the well-established piecewise linear relationship47, 48

[TABLE]

Examples include the energy $E_{N}$, electron density $n_{N}(\boldsymbol{r})$, and atomic kinetic energy $T_{N}$49.

The $C_{6}$ coefficient is not one of those properties. However, we will argue below that the dipole polarizability at non-integer $N$ is piecewise linear and thus does obey Eq. (2). If we then assume that the usual Casimir-Polder formula

[TABLE]

holds for fractional ions, this gives us a recipe for calculating $C_{6}$ coefficients in terms of the imaginary frequency polarizabilities

[TABLE]

Here $\chi_{N}[v_{Z}](\boldsymbol{r},\boldsymbol{r}^{\prime};i\omega)$ is the density response of an $N$ electron atom (where $N$ can be non-integer) in potential $v_{Z}=-Z/r$.

Our argument is as follows: By definition of the density response $\chi=\frac{\delta n}{\delta v}$, we know that $N$ electrons (kept fixed) in a potential $v_{Z}+\delta ve^{\omega t}$ with $\delta v\to 0$ have a density

[TABLE]

where $n(\boldsymbol{r})$ is the density with $\delta v=0$. In the case of a non-integer electron number ensemble we further know that the density

[TABLE]

is piecewise linear between the electron densities of the adjacent integer systems.

But Eq. (6) must hold also for any infinitesimally small change (at least at zero frequency, but most likely at finite frequency too50) to the potential, provided we assume that symmetries are preserved or are correctly accounted for in any approximations employed222Note that our aim here is not to provide a rigorous proof but to provide a justification for the interpolation in practice. Thus we see also

[TABLE]

And if $n^{\prime}$ is piecewise linear and $n$ is piecewise linear, it directly follows that the density response

[TABLE]

must be piecewise linear too, since $n^{\prime}=n+\int\chi\delta v$. Finally, using Eq. (4), we see that

[TABLE]

That is, the polarizability is piecewise linear in the electron number.

This means that we only need to calculate polarizabilities for integer numbers of electrons, and interpolate to find the remaining fractional ionic values. This ability to interpolate from integer cases to non-integer cases is crucial for the success of our method as it keeps the amount of external data sourced from high-level calculations at a manageable size.

2.2 Polarizability model

In the original work of Tkatchenko and Scheffler 11, a simple relationship

[TABLE]

was assumed between the frequency dependent dipole polarizability, the $C_{6}$ coefficient and the static dipole polarizability of the isolated atom. This approximation is based on a one-pole model of the polarizability of an isolated atom, with $\eta$ set consistent with Eq. (3).

In the present work we can no longer rely on Eq. (10) due to the use of linear combinations of polarizabilities, each with a different value of $\eta$. Rather than generalizing Eq. (10), we instead use the two-pole parameterizations of frequency dependent atomic and ionic polarizabilities recently published by two of us41. These models are then combined, as per the piecewise linear formula, into a four-pole model for the fractional ions. This allows us to improve the accuracy of the integer atomic polarizabilities and ensure that our fractional ions are treated correctly.

For each atom and ion with integer electron number $M$ we thus employ a polarizability model

[TABLE]

involving two Lorentzian functions. Then, using the piecewise linearity of $\alpha$, we extend this to arbitrary non-integer electron number $N$ using

[TABLE]

Here we use the properties of systems with neighbouring integer electron numbers $M=\lfloor{N}\rfloor$ and $M+1$, and then linearly interpolate using the remaining fractional charge $f=N-M$. The final model thus employs four Lorentzian functions. The values of parameters $a_{c}(Z,M)$ and $\Omega_{c}(Z,M)$ used in this work are available for neutral atoms and ions of elements from the first six rows of Periodic Table in the Supporting Information of Ref. 41.

Finally, we can perform all integrals analytically in Eq. (3) to derive the $C_{6}$ coefficient (if needed) between two fractionally charged ions $A$ and $B$ in terms of our pretabulated coefficients. Alternatively, in the case of the many-body theory the direct calculation of $C_{6}$ is unnecessary, and instead we can use Eq. (12) directly.

2.2.1 A note on embedded ions

Free-standing ions can be very different to ions embedded in a larger system such as a molecule or a material. For example, the surrounding environment of an anion has a substantial effect on its outermost electron(s). They are only very weakly bound in the free-standing case, with an asymptotic effective potential going to zero as $~{}r^{-3}$. When embedded, the presence of other electrons and nuclei should be represented by an effective confining potential, which can be approximated by a positive power law function of $r$ such as $(r/r_{a})^{\sigma_{a}}$ where $r_{a}$ is an effective embedding radius and $\sigma_{a}$ governs the sharpness of the embedding.

The same electrons that are the most sensitive to the embedding environment, namely electrons in the outermost electronic shell(s), are the ones that contribute the most to the polarizability. Thus the polarizability of an embedded system is highly sensitive to its environment and care must be taken in considering what “ions” we use in our high-level fractional ionic calculations. This is particularly pertinent for open shell systems which are likely to be the most sensitive to the environment.

To this end we use the model polarizabilities from the “minimal chemistry” database of Ref. 41. This provides ab initio dynamic polarizabilities for the first 6 rows of the periodic table that are specifically tailored to provide chemically plausible free anions, calculated as a combination of ensemble DFT51 and time-dependent DFT52, 53. It also allows for the evaluation of double anions, which are problematic in a self-consistent scheme. In short, self-consistent results are included for neutral atoms, cations, and closed shell monoanions; while non self-consistent data (based on the relevant neutral atom) is provided for the remaining anions.

2.3 Eigenvalue remapping

Variants of the MBD method involve (see e.g. discussion in Refs. 20, 54) the solution of the screening equation

[TABLE]

to obtain dispersion energies. Here $\mathsf{A}_{\text{LR}}$ is an atom-wise model of the dipole polarizability in the system which differs only in some details between different schemes (like MBD@rsSCS, and the new MBD@rsSCS/FI model introduced here), and $\mathsf{T}_{\text{LR}}$ is the (long-range) dipole interaction tensor. Calculations thus involve a tractable $3N_{\and}\times 3N_{\and}$ matrix equation, where $N_{\and}$ is the number of atoms/ions in the system. In the case of a periodic system $N_{\and}$ is the number of atoms/ions in the unit cell and we must also include a sum over $\boldsymbol{k}$ points in the irreducible Brillouin zone.

Full details of TS-MBD variants, including the generalization of Eq. (13) to periodic systems, can be found in Ref. 54. In summary, the matrix used in the MBD@rsSCS scheme is defined as $[A_{\text{LR}}(\omega)]_{pq\alpha\beta}=\tilde{\alpha}_{p}(i\omega)\delta_{pq}\delta_{\alpha\beta}$ where $\tilde{\alpha}_{p}(i\omega)$ is the frequency-dependent polarizability of the embedded atom or ion (here $p$ and $q$ refer to atom indices and $\alpha$ and $\beta$ refer to Cartesian indices $x$, $y$ and $z$) obtained by solving the short-range screening equation

[TABLE]

Here $T_{SR}$ is a short-range dipole-dipole interaction tensor. $\tilde{\alpha}(i\omega)$ used to define $A_{LR}$ is one third of trace of the matrix solution of (14).

In conventional MBD@rsSCS the term $\alpha_{p}(i\omega)$ in Eq. (14) comes from Eq. (10). In our new FI scheme, $\alpha_{p}(i\omega)$ is replaced by $\alpha_{Z_{p},N_{p}}(i\omega)$ defined for fractional ions using Eq. (12). The replacement of Eq. (10) by Eq. (12) thus includes embedded ionic charge explicitly [through local electron number $N_{P}$ in Eq. (12)] and implicitly [via the volume scaling in Eq. (1)]. The matrix $[T]_{\text{LR},pq\alpha\beta}$ is used identically in both methods, and is related to the second spatial derivative of a modified Coulomb potential. We refer the reader to Section II of Ref. 54 for furhter details.

Once we have $\mathsf{A}_{\text{LR}}$ and $\mathsf{T}_{\text{LR}}$ we can solve Eq. (13) using the eigenvalues $x_{n}$ of the Hermitian matrix $\mathsf{X}(\omega)=-\mathsf{T}_{\text{LR}}^{\frac{1}{2}}\mathsf{A}_{\text{LR}}(\omega)\mathsf{T}_{\text{LR}}^{\frac{1}{2}}$333Note that identical results are obtained using $\mathsf{A}_{\text{LR}}(\omega)\mathsf{T}_{\text{LR}}$ directly. Indeed our calculations use this direct form. For the manuscript we use the Hermitian form to simplify working.. Here the equalities ${{\text{Tr}}}[f(-\mathsf{A}\mathsf{T})]={{\text{Tr}}}[f(\mathsf{X})]=\sum_{n}f(x_{n})$ give

[TABLE]

The eigenvalues $x_{n}(\omega)$ should obey $x_{n}>-1$ and $\sum_{n}x_{n}=0$ due to constraints and sum rules imposed by the physics of interacting systems.

Unfortunately certain systems, including transition-metal dichalcogenides (discussed in more detail later), highlight a physical and numerical deficiency in the MBD approach that manifests itself most obviously in Eq. (15). In these systems long-range and low-frequency dipole coupling can lead to some values of $x_{n}$ becoming less than $-1$ due to inconsistencies in the model polarizability and screening equations. This problem is related to the well-known phenomenon of “polarization catastrophe” 55, 56 and occurs when atoms are too close together - rendering the Taylor expansion of the Coulomb potential unstable.

To solve the system of self-consistent screening equations in these pathological cases we need to fix these unphysical eigenvalues. Given the success of MBD in many systems, this fix must maintain its good properties in typical cases, but improve on it in the pathological cases. Our solution involves “remapping” the eigenvalues of the system to avoid unphysical cases whilst maintaining continuity in energy and forces with respect to changes in the geometry. Specifically, we replace $\sum_{n}\ln(1+x_{n})$ by $\sum_{n}\{\ln(1+\tilde{x}_{n})-\tilde{x}_{n}\}$ where

[TABLE]

The additional term in the sum ($-\sum_{n}\{\tilde{x}_{n}\}$) is required to correct ${{\text{Tr}}}[\tilde{\mathsf{X}}]$ which, unlike ${{\text{Tr}}}[\mathsf{X}]$, is not guaranteed to be zero.

The remapped eigenvalues $\tilde{x}_{n}$ obey $\tilde{x}_{n}>-1$ for finite $x_{n}$ so that $\ln(1+\tilde{x}_{n})$ is well-defined in all cases. Furthermore, $\tilde{x}_{n}$ is a continuous function of $x_{n}$ with smooth leading derivatives. It thus produces meaningful and smooth energies and forces for all reasonable geometries. The specific form of the ad hoc mapping from $x$ to $\tilde{x}$ is somewhat arbitrary. Equation (16) was chosen as it straightforwardly ensures that: a) $\tilde{x}>-1$ as required, b) $\tilde{x}=x+O(x^{9})$ as $x\to 0$, leading to correct derivatives to all interesting orders and c) $|\tilde{x}-x|<0.01$ for $x>-0.8$, meaning the correction is only significant in extreme cases.

This eigenvalue remapping is mathematically equivalent to transformation of the model polarizability $\mathsf{A}_{\text{LR}}$, albeit in a fairly non-conventional fashion. The final energy expression can thus equivalently be written as

[TABLE]

where $\tilde{\mathsf{A}}_{\text{LR}}(\omega)=-\mathsf{T}_{\text{LR}}^{-\frac{1}{2}}[\sum_{n}\tilde{x}_{n}(\omega)\mathsf{V}_{n}(\omega)]\mathsf{T}_{\text{LR}}^{-\frac{1}{2}}$. Here $[V_{n}]_{ij}=v^{*}_{n,i}v_{n,j}$ for $\boldsymbol{v}_{n}$ the eigenvector of $\mathsf{X}$ with eigenvalue $x_{n}$.

Eigenvalue remapping is an integral part of our new FI scheme, as it ensures that the effective polarizability tensor obeys the necessary constraints. This leads to improved dispersion physics and consequentially better energies when remapping is needed, and almost no change when it is not.

We note that eigenvalue remapping can be optionally linked with the original MBD@rsSCS scheme. This variant, denoted by MBD@rsSCS+ER, is able to fix the pathological behavior of the MBD@rsSCS method in certain, albeit in not all, cases. However, when the solution exists, the results obtained by MBD@rsSCS+ER are generally poorer than those obtained by the new MBD@rsSCS/FI method. Some numerical MBD@rsSCS+ER results are presented in the Supporting Information35. Performance for “well-behaved” systems is affected by eigenvalue remapping only marginally (e.g. Sec. 3.2.1) in both the MBD@rsSCS+ER and MBD@rsSCS/FI approaches.

3 Numerical results

We examine the usefulness of the proposed method in a number of practical applications. First, we will demonstrate that the dipole polarizabilities computed using fractional charges (Eq. (12)) outperforms the approach based on neutral atoms. Second, we combine the frequency dependent dipole polarizibilities with the many-body dispersion correction method of Ambrosetti et al. 33 (PBE+MBD@rsSCS) and show that this new version of the method performs significantly better, especially for for systems with strong ionic behavior.

All calculations were carried out with the periodic DFT code VASP 5.4.157, 58, 59, modified to include the FI approach. The FI approach will be available in the next official release of VASP. A patch is presently available by request to the authors. Setting keyword IVDW=263 activates the FI treatment. For the sake of completeness we note that the combination of keywords IVDW=202 and ITIM=1 activates conventional MBD54 with the eigenvalue remapping, albeit we do not recommend this approach for routine applications.

For the semi-local approximation we use the Perdew-Burke-Ernzerhof (PBE) functional 60. The electronic energies were converged to self-consistency with an accuracy of $10^{-7}$ eV. The $\boldsymbol{k}$-point grids and plane-wave cutoffs used in calculations are summarized in Tab. 1 and Tab.S5 in the Supporting Information35. Drawings of structures presented in this work have been created using the program VESTA 61.

As with conventional TS and its variants, the choice of semi-local approximation is somewhat flexible. For this work we focus on PBE as it allows ready comparison with past calculations. However, the method presented here can, in principle, be coupled with alternative approximations like accurate meta-GGAs62, 63, 64. In such cases the adjustable parameter $\beta$ would need to be recalculated to properly account for the improved short-range interaction in meta-GGA functionals.

3.1 Polarizabilities of atoms in cubic ionic crystals

On the basis of experimental data, polarizabilities of cubic crystals can be obtained from the Clausius-Mossotti equation

[TABLE]

where $V$ is the volume per formula unit (fu) and $\epsilon_{\infty}$ is the measured high-frequency dielectric constant. In this work we use experimental results for from Refs. 65, 66 to compute $\alpha_{m}$ for a set of 18 cubic crystals that ranges from 6.2 au (LiF) to 61.9 au (CsI).

Using the static AIM polarizabilities, polarizability of atoms in crystals are computed by solving the self-consistent screening (SCS) equation described in the Supporting Information35. The following three models to determine AIM polarizabilities (differing in the type of reference object (i.e. atom or fractional ion) and partitioning scheme used) have been considered: rescaling of polarizabilities of neutral atoms using effective volumes obtained from Hirshfeld partitioning (TS), rescaling polarizabilities of neutral atoms using volumes from iterative Hirshfeld partitioning (HI), and rescaling polarizabilities of fractional ions using volumes from iterative Hirshfeld partitioning (FI).

The TS model is used in the MBD@rsSCS scheme of Ambrosetti 33, while the FI model is the basis of our new scheme. The polarizabilities computed using all three models are compared with experimental data in Figure 1, the numerical values for all systems are compiled in Tab. S1 in the Supporting Information35. The TS model strongly overestimates polarizabilities for all systems: mean absolute error (MAE) and mean absolute relative error (MARE) being as large as $41.5$ au/fu (where again fu stands for formula unit) and $157.4\%$, respectively. The use of iterative Hirshfeld partitioning leads to significant improvement but the error still remains very large (MAE=13.7 au/fu, MARE$=43.5\%$). Finally, replacing neutral atoms as reference by fractional ions (FI) leads to a further significant reduction of error with respect to experiment (MAE=6.5 au/fu, MARE$=22.9\%$).

These results also highlight the importance of a full fractional ionic treatment for e.g. self assembly or surface physics problems. The interaction of atoms, layers and molecules with bulks is highly dependent on the crystal polarizability. Here the FI treatment gives a twofold improvement over the HI method and a sixfold improvement over standard TS.

3.2 Molecular systems

3.2.1 Benchmark sets S66$\times$8 and X40

The benchmark set S66$\times$8 of Řezáč et al. 67 consists of configurations from dissociation curves for 66 molecular dimers. For each dimer, eight points differing in the distance between monomers ranging between 0.90$\times$R0 to 2.00$\times$R0 (R0 is the ground-state separation) are considered. The interaction types covered by this benchmark set include hydrogen bonding, dispersion interactions and other interaction types such as X-H$\cdots\pi$ (X=C,N,O) interactions, as well as nonspecific interactions of polar molecules. The list of all 66 dimers is presented in Tab.S2 in the Supporting Information35. Minimization of error in the interaction energy of the S66$\times$8 set with respect to high-level reference data has been used by Ambrosetti et al. 33 to optimize the value of free parameter of the MBD@rsSCS method ($\beta$) and we have followed the same strategy to determine $\beta$ for the MBD@rsSCS/FI method.

The optimal value $\beta$=0.83 found for MBD@rsSCS/FI is identical to that used in the model based on neutral reference atoms. As shown in Tab. 2, the method based on fractional ions yields results that are consistently better for all interaction types covered by this benchmark set than those obtained by the MBD@rsSCS. The difference is particularly notable for the dispersion interaction dominated subset where MAE decreased from 17.5 meV to 15.4 meV. It is also encouraging that the MBD@rsSCS/FI method improves over the MBD@rsSCS for all separations of monomers (see Tab. S3 in Supporting Information35). The only exception is the largest distance (R=2.00$\times$R0) where the interaction energies are small in most cases and where the performance of MBD@rsSCS and MBD@rsSCS/FI methods is very similar.

The set of noncovalent molecular dimers X40 of Řezáč et al. 68 was designed to benchmark the quality of theoretical methods for description of halogen atoms and covers the interaction types such as London dispersion, induction, dipole-dipole, stacking, halogen bonds, halogen-$\pi$ bonds, and hydrogen bonds.

As a reference, the interaction energies computed at the CCSD(T) level are used. The computed interaction energies for all systems are compiled in Tab. S4 in the Supporting Information35, the statistics obtained by the two variants of the MBD@rsSCS method is compared in Tab. 2. The use of the model based on fractional ions leads to improvement for all interaction types with exception of the dipole-dipole interaction where the results are practically identical to those of the original method. The most significant improvement is found for stacking interactions (MAE reduced from $16.8\%$ to $13.3\%$) and halogen-$\pi$ bonds ($22.5\%$ (MBD@rsSCS) vs. $19.3\%$ (MBD@rsSCS/FI)). The decrease of error for the whole set is also quite significant (MAE: 13.9 meV vs. 15.5 meV; MARE: $9.0\%$ vs. $10.0\%$).

3.2.2 Interaction of H2 with substituted coronenes

The investigation of hydrogen storage devices represents one of the important applications of dispersion corrected DFT methods. Kocman et al. 69 proposed a set of six adsorption complexes consisting of H2 molecules and coronene-derived structures as a benchmark for methods used in this type of applications. The coronene-derived substrate models contain two purely carbonaceous molecules (Coronene…2H2 and C-coro…2H2) and four boron substituted coronene molecules (CoroB2Li2…2H2, CoroB2…2H2, CoroB2Li2…H2os, CoroB2Li2…H2ss) out of which three models contain two additional Li+ cations exhibiting relatively strong affinity to molecular hydrogen. The reference interaction energies determined at the CCSD(T) level range between 4.7 kJ/mol (Coronene…2H2) to 14.3 kJ/mol (CoroB2Li2…2H2) and the two extreme cases are shown in Fig. 2 (all structures are presented in Figs. S1-S3 in the Supporting Information35).

The MBD@rsSCS method based on the neutral atomic reference works quite well for the models that do not contain Li atoms. As shown in Tab. 3, the error with respect to the CCSD(T) reference is only 1.2 kJ/mol or less in these cases. The presence of Li atoms, however, causes serious numerical issues. In particular, solution of the screening equation Eq. (14) yields negative polarizabilities for some atoms and these can not be handled in the present version of the MBD@rsSCS scheme.

This problem is caused by a very large AIM polarizabilities of the Li atoms that is a consequence of fundamental property of the Hirshfeld partitioning that the atoms-in-molecules are as similar to the reference objects (neutral atoms) as possible in the information-theoretical sense 70. Replacing the neutral atomic reference by fractional ions fixes this problem: in accord with expectation, the AIM polarizabilities of the Li in substituted coronenes ($\sim$5 au) are now closer to polarizability of non-interacting Li+ (0.2 au 41) than to that for free Li0 (164 au 41). Importantly, the quality of the MBD@rsSCS/FI results for the Li-containing models is consistent with that for the purely carbonaceous models (see Tab. 3).

3.3 Crystalline materials

3.3.1 Benchmark set X23 consisting of molecular crystals

In this section we will discuss the performance of different variants of the MBD@rsSCS method in optimization of systems from the X23 set of Reilly et al. 71 (based on previous work of Otero-de-la-Roza and Johnson 72) consisting of molecular crystals with cohesion dominated by London dispersion and H-bonding interactions. Full optimization of atomic and lattice degrees of freedom has been performed and cohesive energies were computed for the final geometries. Computational details and additional results are presented in the Supporting Information35.

In our previous work 54 we have shown that the neutral atom based model provides already a very good agreement with the reference energetic and structural data. The replacement of neutral atoms by fractional ions in the polarizability model leads to a modest improvement in energetics (MARE$=7.4\%$ (MBD@rsSCS) vs. $6.9\%$ MBD@rsSCS/FI) and almost negligible variation in predicted lattice structure, see Tab. 4 and the Supporting Information35 (Tabs.S6 and S7).

3.3.2 Ionic crystals

Crystalline sodium chloride (space-group $Fm\bar{3}m$) represents highly ionic material in which the Born effective charges of Na and Cl ions (Z∗=$\pm$1.1 e) are close to the formal charges $\pm$1 e. 73 In our previous work 73 we have shown that the Hirshfeld partitioning scheme based on neutral atoms strongly underestimates the ionic charges (qAIM=$\pm$0.2 e). Consequently, and in contrast to expectation, the predicted properties of Na and Cl embedded in the interacting systems are more similar to those of neutral atoms than to those of ions.

In this case, the pairwise Tkatchenko-Scheffler method based on neutral atoms 36 leads to a strong overestimation of interaction between atoms in NaCl and a strong underestimation of the lattice constant [5.34 Å (TS) vs. 5.57 Å (exp.)]73. The consequences for the MBD@rsSCS scheme based on the same polarizability model are even more serious: the method suffers from numerical issues (already discussed in Sec. 3.2.2) that prevent its application to systems of this type. The use of fractional ions in the MBD@rsSCS/FI scheme fixes this problem and the computed lattice constant and bulk modulus are in very good agreement with experiment, see Tab. 5.

Magnesium oxide (MgO) is another example of ionic material which we will consider here. The crystalline MgO is isostructural to NaCl and this material is frequently used as a catalyst or adsorbent in many applications in chemistry74. The standard MBD@rsSCS fails when used to optimze MgO, while the MBD@rsSCS/FI method yields a lattice constant and bulk modulus that are in good agreement with experimental data75. Similar quality results are found for LiF.

In all three tested ionic solids PBE is accurate, yielding lattice constants and bulk moduli only slightly worse than MBD@rsSCS/FI. However, when considering interactions with surfaces, such as adsorption or surface reactions74, it is important to treat the whole system self-consistently using a method with vdW corrections. Unlike MBD@rsSCS, MBD@rsSCS/FI is thus accurate enough to use for surface studies.

Next we shall consider mineral cryolite (Na3AlF6), which is an example of ionic system with a complex structure. The low temperature phase $\alpha$ that we will discuss here has a symmetry $P2_{1}/n$ (monoclinic) and consists of distorted octahedra AlF6 and NaF6 linked via common corners occupied by the fluorine atoms, see Fig. 3. The vacancies between these octahedra are occupied by cations Na+.

A recent theoretical study 76 showed that the semilocal DFT (PBE) tends to overestimate the cell volume by $\sim~{}3\%$ and the error even increases in the high-temperature phase $\beta$ and in the liquid phase. The missing dispersion interactions have been suggested as a possible reason for this volume overestimation. The lattice geometry determined using the MBD@rsSCS/FI ($a=5.39$ Å, $b=5.61$ Å, $c=7.78$ Å, $\beta=90.3$∘) is close to the experimental geometry for T=0K 77 ($a=5.38$ Å, $b=5.58$ Å, $c=7.69$ Å, $\beta=90.3$∘) and the error in the computed cell volume is $<1\%$. By contrast, the attempted MBD@rsSCS calculation failed for the same reason as in the case of Li-substituted coronenes and NaCl.

3.3.3 Dichalcogenides

In our previous work, we have shown that the conventional MBD@rsSCS scheme works very well in description of the structure and energetics of graphene bilayers and graphite 54. Here we focus on the properties of bilayers and bulk systems of dichalcogenides (MoS2, MoSe2, WS2, and WSe2), which can be considered as van der Waals systems with a partially ionic character.

For the bilayers of dichalcogenides discussed here, the two most favorable stacking patterns AA’ and AB 78 are considered and the RPA results of He et al. 78 are used as a reference. In the calculations, the interlayer distance was varied while the geometry of individual layers was fixed at the value determined experimentally for the bulk systems  79.

Because of the problems with incorrect eigenvalue spectrum of the $\mathsf{A}_{\text{LR}}(\omega)\mathsf{T}_{\text{LR}}$ matrix (see Sec. 2.3), the original variant of the MBD@rsSCS method can not be applied directly to the systems discussed here, hence we discuss only results obtained using MBD@rsSCS/FI method in which this problem is overcome, as described in Sec. 2.3.

The computed interaction energy profiles for the AA’- and AB-stacked MoS2 are compared with the RPA reference of Hu et al. 78 in Fig. 4; the profiles for the other systems considered here (MoSe2, WS2, and WSe2) are presented in Figs. S4-S7 in the Supporting Information35. As evident from the comparison of the energy profiles and from the numerical values of binding energies and equilibrium interlayer separations (see Tab. 6), the MBD@rsSCS/FI method reproduces the RPA results very well and the improvement over the pairwise version of the method (PBE+TS 36) is remarkable.

The calculations of ground-state structure and binding energies of bulk dichalcogenides with AA’ stacking have been performed via series of single-point energy calculations in which the lattice parameter parallel with the stacking direction ($c$) was varied in the interval between $\sim$10.7 Å to 44.2 Å. In these calculations, all intralayer structural parameters and lattice vectors parallel to the layers were fixed at their experimental values.

The computed results presented in Tab. 7 are qualitatively consistent with our results for bilayers and the MBD@rsSCS/FI results are reasonably close to the RPA reference. The standard MBD@rsSCS method has the same numerical issues (from the polarization catastrophe) as in the case of bilayers. We thus cannot calculate or report meaningful results from it.

3.3.4 Graphite fluoride

In fluoride derivative of graphite, the fluorine atoms are covalently bound to carbon atoms, see Fig. 5. In this work we use the MBD@rsSCS approach to compute binding energy for the structure with the AA stacking of layers which has recently been identified as the most stable structural arrangement for this material 80. The simulation cell used in calculations contained four atoms belonging to a single layer and the interlayer separation was controlled via the length of the lattice vector normal to the layer ($c$). Following Lazar et al. 80, atomic positions and length of lattice vectors parallel with the layer ($a$=2.582 Å) were fixed during the calculation.

The dependence of the binding energy on the interlayer distance computed using the MBD@rsSCS and MBD@rsSCS/FI methods is compared with the RPA reference 80 in Fig. 6. The MBD@rsSCS and MBD@rsSCS/FI methods predict a similar ground state value for the interlayer separation (5.91 Å (MBD@rsSCS), 5.89 Å (MBD@rsSCS/FI)) that are close to the RPA reference (5.88 Å). The binding energy is somewhat underestimated in both cases but the MBD@rsSCS/FI value (28.7 meV/fu) is closer to the RPA reference (34.7 meV/fu) than the MBD@rsSCS result (25.3 meV/fu).

4 Limitations and perspectives

Before concluding, we discuss a few limitations of our FI approach which warrant study in future works. Many of these limitations are shared by other van der Waals functionals.

Firstly, the FI approach neglects quadropolar and higher order van der Waals terms. This can possibly be remedied by approximating these terms81 or pre-calculating them directly in higher level theory as in Ref. 41. However, we note that in layered geometries the higher order contributions appear to be heavily screened in the near-contact regime31, perhaps explaining the success of methods which neglect them.

Secondly, the lack of long-range fluctuations in the scheme means it cannot properly capture Dobson-C23 effects, a property it inherits from the MBD scheme. This is likely to be important in long-range forces from metallic systems16, 22 and systems with unusual topological properties26. However, there is some evidence (e.g. 29) that some of these effects might be captured via the coupled fluctuating dipole model’s inclusion of Dobson-B terms, at least in the medium range.

Finally, the eigenvalue remapping scheme is somewhat ad hoc. Further testing of how the polarization catastrophe behaves in exemplar (non-remapped) systems should shed light on how to include more physically intuitive corrections. Work in this direction is being pursued.

5 Conclusions

In this work we introduced a scheme (MBD@rsSCS/FI) for calculating polarizabilities and dispersion forces in electronic systems. Our approach uses the properties of fractional ions (FIs) in a point-dipole dispersion approximation based on the MBD approach. Key equations describing our approach are given in Eqs. 1, 4 and 16, and related discussion.

Our fractional ionic approach vastly improves results for ionic systems in bulk and layered materials, compared to related approaches10, 11, 20, 40. It is particularly successful in transition-metal dichalcogenides like MoS2, and in interactions of H2 with modified coronenes, systems of nanotechnological interest. It marginally outperforms similar schemes in most of the non-ionic system sets tested, and is of only marginally lower accuracy in the remainder.

Our approach is thus a suitable alternative to existing schemes in systems where they work, while offering substantial improvements in ionic systems where they do not. Therefore we strongly advocate in favour of the PBE+MBD@SCS/FI method in routine applications (using keyword IVDW=263 in VASP).

{acknowledgement}

T.G. received computing support from the Griffith University Gowonda HPC Cluster and a Griffith University International Travel Fellowship. Part of this work has been supported by the VASP project. T.B. is grateful to University of Lorraine for invited professorships during the academic year 2015-16 and he acknowledges support from Project 643 No. APVV-15-0105. Calculations were performed using computational resources of University of Vienna, and supercomputing infrastructure of Computing Center of the Slovak Academy of Sciences acquired in projects ITMS 26230120002 and 26210120002 supported by the Research and Development Operational Program funded by the ERDF. J.G.A. is thankful for the partial support of this research by the European Union and the State of Hungary, co-financed by the European Social Fund in the framework of TÁMOP 4.2.4. A/2-11-1-2012-0001 National Excellence Program. S. L. acknowledges HPC resources from GENCI-CCRT/CINES (Grant x2017-085106).

{suppinfo}

We include with this manuscript supporting information including:

• A PDF document that includes: further theoretical details of methods; values of computed atoms-in-molecule polarizabilities; and numerical results for S22 and X40 benchmark sets.