Improving the description of interlayer bonding in TiS2 by Density Functional Theory

We investigate energetic and electronic properties of TiS2 , an archetypal van der Waals (vdW) material, from first principles, in the framework of the Density Functional Theory (DFT). In this system a recent experimental study showed a puzzling discrepancy between the distribution of the electron density in the interlayer region obtained by X-ray diffraction data and that computed by DFT, even adopting DFT functionals that should properly include vdW effects. Such a discrepancy could indicate a partial failure of state-of-the-art DFT approaches in describing the weak interlayer interactions of TiS2 and, possibly, of similar systems too. In order to shed light on this issue, we have carried out simulations based on different DFT functionals, basically confirming the mentioned discrepancy with the experimental findings. Subsequently, we have tried to reproduce the experimental interlayer electronic density deformation both by changing the parameters characterizing the rVV10 DFT functional (in such a way to artificially modify the strength of the vdW interactions at short or long range), and also by adopting a modified pseudopotential for Sulfur atoms, involving d orbitals. The latter approach turns out to be particularly promising. In fact, using this novel, more flexible pseudopotential, we obtain not only an electronic density deformation closer to the experimental profile, but also a better estimate of the interlayer binding energy. Interestingly, this improvement in the theoretical DFT description is not limited to TiS2 but also applies to other similar layered systems involving S atoms, such as TaS2 , HfS2 , and MoS2 .


Introduction
First identified in 1873, 1 the van der Waals (vdW) interactions are forces that today attract more interest than ever. They are present everywhere, but their manifestations still pose challenging questions which are relevant for such varied systems as soft matter, surfaces, and DNA, and in phenomena as different as supramolecular binding, surface reactions, and the dynamic properties of water. vdW interactions play a central role even in the characterization of layered materials, such as transition metal dichalcogenides (TMDs). These represent a large class of materials with mild stiffness, which are not as soft as tissues and not as strong as metals. 2 They are two-dimensional systems which have unique properties and are central to current research in solid-state science. In TMDs, a transition metal atom (M) layer is sandwiched between two chalcogen atom (X) layers and it is commonly assumed that the MX 2 slabs are stacked by vdW interactions, whereas the intralayer M-X interactions are instead covalent. The weak, interlayer vdW interactions play a key role in the formation, intercalation, exfoliation and layer-by-layer building of TMD materials, as well as being decisive for their characteristic properties.
Like all non-relativistic electronic effects, the vdW interactions are present in the (unknown) exact Density Functional Theory (DFT) exchange-correlation (XC) functional, but they are not properly described by standard approximate DFT schemes, such as local density approximation (LDA) 3 or semilocal generalized gradient approximation (GGA). 4 Hence the need to introduce truly non local XC functionals, able to accurately account for vdW corrections, arises. Despite a considerable amount of work focused on the development of accurate vdW-corrected DFT methods, a precise characterization of the layer-layer distance and of the interlayer binding energy of vdW layered materials, still represents a challenging issue. [5][6][7][8][9] This is due to the nonlocal and long-range nature of the vdW interaction, as well as to the coexistence of the weak vdW bonding and much stronger intralayer chemical bonding. Fully accounting for the vdW interaction is achievable by high-level methods such as quantum Monte Carlo (QMC), 10 coupled-cluster singles and doubles with perturbative triples (CCSD(T)), 11 and exploiting the adiabatic-connection fluctuation-dissipation theorem within the random-phase approximation (RPA), 12 which are, however, only feasible for limited-size systems because of their high computational cost.
The physical quantity lying at the DFT core, is the electron density (ED), which fully determines the total energy of any many-electron system in its ground state. The ED distribution probably represents the most information-rich observable available since it has become possible to determine EDs from analysis of structure factors obtained from accurate X-ray diffraction data. During the past decade, the accuracy of experimental X-ray diffraction data has increased dramatically owing to the use of high-energy synchrotron sources, which significantly limit systematic errors in the data and, thanks to these modern techniques, nowadays it is possible to evaluate both structural and electronic properties of a wide class of layered materials in an extremely accurately way. In this respect, Medvedev et al. 13 have pointed out that modern DFT functionals are constructed on the basis of empirical fitting on energetic and geometrical benchmarks, neglecting the ED as a parameterization parameter.
If, from one side, this leads to improved predictions of energetic and structural features, from the other, the exact reproduction of the ED has worsened. In a recent study, Kasai et al. 14 studied the archetypal vdW material, TiS 2 , to compare the ED distribution derived from X-ray diffraction data with the theoretical one obtained by several DFT functionals.
They obtained a good agreement between theory and experiment for the description of the intralayer, covalent Ti-S interaction, but, at the same time, significant discrepancies were found for the interlayer vdW interactions, particularly considering the ED distribution in the region between S atoms belonging to adjacent layers. In fact, while the properties at the bond critical point (BCP) were quite similar, noticeable differences were instead observed away from the BCP: considering the theoretical DFT density, the S atoms behave as if they are practically not interacting with the neighboring layer, showing charge concentration and accumulation in the region expected for a lone pair in an sp3 -hybridized S atom with three bonds to nearby Ti atoms and no other bonds; in the experiment, however, an appreciable charge accumulation and concentration was observed in the interlayer region, which was interpreted as a sign of a stronger and more directional interaction between S atoms. Thus, this study suggested the inability of current DFT functionals to accurately describe the interlayer ED for vdW layered materials, thus supporting the conclusions of Medvedev et al.
In order to shed light on this issue, in the first part of this study, we apply the LDA, the semilocal PBE 4 and the vdW-corrected rVV10 15 DFT functionals, with standard Projector-Augmented Wave (PAW) pseudopotentials, 16 to explore the behavior of the interlayer ED in different regimes, obtained by also modifying the parameters of the rVV10 functional, so that to artificially tuning the intensity of the vdW interactions. The basic result is that, in order to produce an ED distribution appreciably closer to the experimental profile, we have to select the rVV10 parameters in such a way to unphysically increase the strength of the vdW interactions, thus leading to an exaggeratedly high interlayer binding energy. In the second part a standard rVV10 approach is used, but a novel pseudopotential is adopted, which explicitly accounts for the 3d orbitals of the S atom. While these orbitals are empty in the isolated, neutral S atom, they can play a role 17 in molecules and condensed-matter systems where S atoms interact with other atoms and form bonds. Interestingly, using this novel, more flexible S pseudopotential, we obtain not only an ED distribution closer to the experimental profile, but also a better estimate of the interlayer binding energy. We finally show that this improvement in the theoretical DFT description is not limited to TiS 2 but also applies to other, similar layered systems involving S atoms, such as TaS 2 HfS 2 , and MoS 2 .

Theoretical background and Computational details
In DFT the properties of a many-electron system can be determined by introducing suitable functionals of the electron density n(r); the most interesting many-body physical effects are described by the so-called exchange-correlation (XC) functional, related to the XC potential: V xc [n(r)] = V x [n(r)] + V c [n(r)]. In principle V xc describes all the many-body effects, including the vdW interactions, however, in practice, a number of approximation must be introduced, since, unfortunately, the exact form of the XC functional is not known. Nowadays several kinds of approaches to approximately include vdW interaction in DFT calculations are available. Some of these are based on semiempirical, atom-based pair potentials, able to give approximate vdW corrections; 18-20 instead others are built introducing a suitable nonlocal XC density functional, such as the vdW-DF family by Dion et al. 21 and the scheme proposed by Vydrov and Van Voorhis, 22 with the revised, more efficient version rVV10. 15 rVV10 certainly represents one of the best vdW-corrected functionals nowadays available; it takes into account the entire range of vdW interactions at a reasonable computational cost using only the ED n(r) as input. In particular, the rVV10 non-local correlation energy has the following expression: where Φ is the non-local correlation kernel, which is: In the above expression R = |r − r'| and q(r) is a function of the ED and its gradient: where and k(n(r)) = 3πb n(r) 9π 6 The total XC functional is then obtained by adding the non-local correlation energy to the refitted Perdew-Wang exchange functional 23 The rPW86 exchange functional was chosen mainly because it is nearly vdW-free, 23 in such a way to avoid double-counting effects. In eq. (4) the C parameter represents the so-called local band gap and is tuned to give accurate asymptotic vdW C 6 coefficients, thus regulating the behavior of the long-range component of vdW interactions. The fitting procedure, aimed to minimize the average error for a benchmark set of 54 C 6 coefficients, was originally carried out by Vydrov and Van Voorhis, 22 leading to an optimal value of C = 0.0093. Another essential aspect of the VV10 (and its successor rVV10) formalism is the presence of a second adjustable parameter, b (see eq. (5)), which controls the short range behavior of the non-local correlation energy. This means that when E N L c is added to the other energetic terms, the b parameter is adjusted to merge the interaction energy contributions at short and intermediate ranges.
With an empirical fitting procedure on the S22 binding-energy data set, 25 Vydrov and Van Voorhis proposed a value of b = 5.9. 22 After the implementation of the efficient Roman-Perez, Soler 26 interpolation scheme in the reciprocal space, this value was revised by Sabatini and coworkers to b = 6.3. 15 We have chosen rVV10 as the basic DFT functional to investigate the effect of vdW interactions in the ED distribution of TiS 2 (although extensive tests were also performed using other vdW-corrected functionals), also because of the possibility to separately modify the intensity of both short-and long-range vdW interactions in a transparent way by simply tuning the two adjustable parameters b and C mentioned above. First the short-range component is analyzed. This is done by considering two different regimes: the former in which the value of the b parameter is increased to b=10.0, thus damping the intensity of short range vdW interactions, and the latter in which b is reduced to b=1.0, that instead results in a substantial increase of the intensity of the short-range vdW interaction. While tuning the parameters of the rVV10 functional can improve the description of some properties of the TiS 2 system, clearly this could lead to a reduction of the transferability of the modified functional to other systems.
Our ab-initio calculations have been performed using the version 6.5 of the Quantum ESPRESSO (QE) DFT package. 27, 28 We first used for Ti and S atoms standard PAW 16 pseudopotentials taken from the QE database. 29 Subsequently, for further testing, we have generated, for the S atom, a novel pseudopotential, using the atomic QE program. We report DFT results obtained using the LDA, the semilocal GGA PBE functional, and the vdW-corrected functional rVV10, both in its standard form and considering the variants (parameter changes and combination with a novel S pseudopotential) described above. Selfconsistent calculations and relaxation processes, leading to the ground-state equilibrium geometry of bulk TiS 2 , were carried out within the PWscf QE program. The calculations adopted plane-wave and density cutoffs of 100 and 1000 Ry, respectively, to get highlyconverged results, a cold-smearing parameter of 0.01 Ry, and a 12×12×4 k-point mesh for the sampling of the Brillouin Zone (BZ). The ED analysis is carried out using the PostProc QE program, employing the B-Spline method 30 for the n(x, y, z) interpolation; in particular, for the ED profile computed along the S-S, S-center, and Ti-S lines, a real-space grid of 600×600×800 (x, y, z) points, generated inside a parallelepiped containing the axis of interest, is used. To get a reliable estimate of the ∇ 2 n(r) function, in the case of the S-S line, the distance between the parallelepiped axis and the edges of its sides is 0.5 A, for the S-center case 0.25 A, and finally 0.35 A for the Ti-S line. The ED is then processed by using the Critic2 program. 31,32 More specifically, we compute the ED charge deformation, ∆n(r), by considering a spatial mesh of 72×72×108 points inside the unit cell (containing one Ti atom and two S atoms) and taking the difference between the ED n(r) and the sum of the EDs of the isolated Ti and S atoms. ∆n(r) profiles and Laplacian ∇ 2 n(r) maps, in the interlayer and intralayer areas, are then evaluated on the plane of interest using Critic2.
The calculation of the atomic properties of a given S atom is carried out by making an integration of the so-called "atomic basin" of the S atom, represented by the space region whose surface has zero-density gradient flux. This is done again with the Critic2 program, that allows to compute spherical multipolar moments defined as where the integration domain A is the atomic basin volume. In the above equation, n(r) is the charge ED and Y lm are the spherical harmonics, in which the m index runs over the interval [-l ; l]. With this definition, the magnitude of the dipole moment attributed to the S atom is obtained as the square root of the sum of the squared l=1 terms:

Results and discussions
Starting from the experimental lattice parameters, a = 3.398 A and c = 5.665 A, as reported by Kasai et al., 14 we relax the crystal structure in the intralayer atomic plane until the force acting on the system is smaller than 10 −5 Ry/a.u. This relaxation process, carried out with both PBE and rVV10 functionals, leads to very similar results, as expected, 14 since vdW interactions do not substantially influence the geometry of the intralayer plane. By further optimizing the interlayer geometry, the standard rVV10 functional predicts that the distance between two S atoms belonging to adjacent layers is 3.443 A while the intralayer S-Ti distance is 2.421 A, in excellent agreement with the reference DFT values (at the SCAN+rVV10 level 33 ) of 3.446 A and 2.420 A, respectively, reported in ref. 14. As already mentioned, Kasai et al. 14 pointed out that DFT simulations are characterized by a an interlayer S-S ED deformation significantly lower than the experimental one, implying that the X-ray diffrac-  Table S1 of Supporting information. Two non-equivalent (3, -3) CPs are found along the Ti-S line (see top left panel in Figure S1 of Supporting information). The (3, -1) CP is instead located at 1.12 A from the Ti atom. As can be seen, looking at Table S1 of Supporting information, very similar values are predicted by all the DFT functionals applied and most of these are also in good agreement with the experimental data.
A similar analysis can be also applied to characterize the interlayer bonding: the properties of the (3, -3) CP for both the S-S and S-center lines, and for the (3, -1) CP along the the S-S line are summarized in Table 1 and 2, respectively. In the interesting case of the S-S line, the PBE, rVV10, and rVV10(b=10.0) functionals predict a maximum value of the ED at a distance r max ∼ 0.61 A from the reference S atom, with a relative error of ∼ 13 % with respect to the experimental result, 14 while the error is slightly larger with rVV10(b=1.0) and LDA functionals.  The small oscillations of the Laplacian profile are due to the finite-size of the real space mesh considered for the ED computation. In this case ∇ 2 n(rmax) and ∇ 2 n(r min ) are estimated by a parabolic interpolation f (r ∼ c) = c 2 (r−c) 2 +k, where the expansion center is c = rmax or c = r min . c) ED deformation ∆ n(r) computed using the PBE functional and d) the rVV10(b=1.0) functional. Curves obtained by rVV10 and rVV10(b=10.0) functionals are very similar to the PBE ones and therefore are not reported.
The n(r max ) values are all very similar and in reasonable agreement with the experimental data, with the LDA functional that shows the best agreement with the experimental value. The rVV10(b=1.0) result, for which the maximum discrepancy occurs, has a deviation of only ∼ 9 %. The Laplacians ∇ 2 n(r max ) are also in acceptable agreement with the experimental values. Again, the result of the rVV10(b=1.0) functional exhibits the maximum deviation from the experiment, with a discrepancy of ∼ 8 %. This means that, at a distance r max from the reference S atom, along the S-S line, with rVV10(b=1.0), the ED deformation is slightly smaller than for the other functionals.
For the S-center line a direct comparison with experimental data is not available, so we take as reference the theoretical SCAN+rVV10 result reported by Kasai et al. 14 The S-center line is characterized by the same r max value found for the S-S line, again with a small deviation of the value obtained with rVV10(b=1.0) and LDA. Even in this case, the results for the ED and the Laplacians are in line with the theoretical reference data. 14 Interestingly, for each functional employed, a slightly larger n(r max ) value than that relative to the S-S line is observed, with a stronger ED deformation described by the associated Laplacian, indicating that the charge deformation along the S-center line is stronger than for the S-S line. The maximum deviation from the reference values, for both ED and Laplacians, occurs again for the LDA functional with an error of ∼ 8 % and ∼ 28%, respectively. The S-S bond shows a symmetric ED distribution along the line connecting the two S atoms (see Figure 2a) CP accompanied by a smaller Laplacian value ∇ 2 n(r min ) S − S (see Table 2) that indicates a slightly less charge depletion than for the other functionals considered.
The ED deformation map in the interlayer plane, presented in Figure 3, is very similar with the PBE, rVV10, and rVV10(b=10.0) functionals (as well as the LDA case which it is not reported), while the same is not true for rVV10(b=1.0), which predicts a more pronounced ED deformation in the region between the two S atoms belonging to adjacent layers. The Laplacian map of the interlayer plane is presented in Figure S2 of the Supporting information: in this case, no significant differences for the functional schemes adopted up to now are found, as well as for LDA and the other schemes presented below.    Table 3. We optimized the interlayer lattice constant c for each functional adopted (third column of Table 3 can reach chemical accuracy for a broad variety of systems (from small molecules to large nanostructure and periodic materials). 36 In particular, many-body effects were shown to be crucial in two-dimensional vdW materials, 37 as they can qualitatively alter the power law scaling of dispersion interactions. Remarkably (see Table 3) the MBD estimate is very close to the RPA value reported in ref. 34.
As expected, PBE dramatically underestimates the ILBE (the discrepancy is even worse considering the optimal configuration predicted by PBE), a behavior that can be clearly ascribed to the inadequate inclusion of vdW effects by PBE. As can be seen in Table 3 Looking at Table 3 Table 1), confirms that the ED predicted at DFT level along the S-S line is smaller than the one along the S-center line, at variance with the experimental evidence. Small, quantitative differences between our calculations and the theoretical SCAN+rVV10 results reported by Kasai et al. 14 can be easily explained by the different DFT functional adopted and by other technical differences: for example, in ref. 14 a finer mesh of 22×22×13 k-points was used for the sampling of the BZ, this choice being motivated by the need of describing core electron oscillations to carry out a multipole refinement procedure. 14 Coming to a more detailed analysis, we observe (see Figure S2 and Figure S3 of Supporting information) that the maps of ∆n(r) and ∇ 2 n(r) obtained by rVV10(b=10.0) are very similar to those given by the standard rVV10 functional and also by PBE, thus suggesting that the interlayer ED distribution is weakly affected by an augmented damping applied to the shortrange vdW interactions; this is also quantitatively confirmed by the the values of n(r) and ∇ 2 n(r max ) evaluated at the (3, -3) CP (see Table 1 )  Table 1 are in line with the other functionals employed. Similar conclusions hold considering the S-S and S-center lines. Finally, one must point out that, even with this new functional, the ED accumulation along the S-S line is predicted to be smaller than along the S-center line, differently from the experimental findings. Therefore, considering that (see Table 3 ) rVV10opt still largely overestimate the ILBE, although it slightly reduces the dramatic overbinding of rVV10(b=1.0), one can conclude that it is not possible to achieve, at the same time, a description of energetic and ED deformation features, in reasonable agreement with the experimental results, by simply tuning the values of the b and C parameters of the rVV10 functional, which determine the strength of short/mediumand long-range vdW effects.
As an alternative strategy to reduce the discrepancy between the theoretical DFT description of TiS 2 and the experimental evidence, we have also investigated the effect of introducing a modified pseudopotential for S. The S atom is characterized by the [Ne]3s 2 3p4 ground-state configuration, so that in neutral S atom d orbitals are empty, and, in simple molecules, they are nearly unpopulated. Nonetheless, it was shown 17 that they can play a significant role: for instance, even in the small S 2 molecule an appreciable improvement in the binding-energy estimate was observed by adopting for S a pseudopotential built by taking into account d orbitals, which are expected 17 to give S enhanced polarizability and enable a great variety of bonding configurations through hybridization. In fact, although some orbitals (for instance, in this case the d ones) are not bound states of the neutral atom, it is necessary 17 to choose a partially ionized atomic reference configuration, with a charge slightly different from the neutral atom, to get pseudopotentials that are better tailored for applications to solid-state and extended systems, since the effective configuration of these systems differs from that of the neutral atom.
Following this approach, we have generated a novel S ultra-soft pseudopotential, based on the configuration [Ne]3s 2 3p 3.0 3d 1.0 , using the Troullier-Martins pseudization method 39 with non-linear core corrections. 40 Coupling it with the rVV10 functional, we have developed the rVV10d scheme.
First we have applied the rVV10d functional to the simple case of the S 2 molecule, which is in itself an interesting system, being an ubiquitous intermediate in the combustion, atmosphere, and interstellar space. 41 With rVV10d, while the equilibrium distance is unchanged with respect to rVV10, we obtain a substantial improvement of the molecular binding-energy and a slight improvement of the vibrational frequency (see table S3 of Supporting information), thus confirming the importance of including the d orbitals in such a system. In particular, considering the binding energy, the more than 20 % discrepancy with the experimental reference value of the rVV10 original functional, reduces to only about 3 % with rVV10d. We have therefore recomputed the energetic and electronic properties of TiS 2 using the new rVV10d functional.  As can be seen in Table 4, at the (3, -3) and (3, -1) CPs, with rVV10d, the basic ED properties are very similar to those obtained by the standard rVV10 functional. However, the effect of replacing rVV10 with rVV10d is evident looking at Figure 5: one can appreci-ate a significant ED accumulation along the S-S line, which is not present using the other DFT functionals considered, including that (SCAN+rVV10) adopted by Kasai et al. 14 Although the discrepancy with the experimental description is only partially resolved (even with rVV10d the ED accumulation along the S-S line is predicted to be slightly smaller than along the S-center line), a clear improvement is evident, which we attribute to the increased flexibility of the new S pseudopotentials that is able to better reproduce the variety of bonding structures and hybridization processes involving S atoms. with the rVV10d and compared with rVV10. The notation used is the same as in Table 3. Remarkably, with rVV10d, an even more substantial improvement in energetic properties takes place. In fact, while the CE remains very close to that obtained by rVV10, suggesting that d orbitals of S do not play a significant role in the formation of strong, intralayer covalent bonds, the change in the ILBE is instead considerable: assuming the experimental structure, our computed ILBE differs by only ∼ 2 meV/ A 2 from the theoretical references, i.e. high quality, RPA, 34 and MBD estimates, which represents a relative error of about 10 %; such an error is comparable to that obtained by the much more expensive meta-GGA approaches and is much smaller than that (about 47 %) found with the original rVV10 functional which clearly overbinds. As can be seen in Table 5, the good energetic description of rVV10d is preserved even considering a TiS 2 configuration where the experimental lattice constant c is replaced by the c value optimized with rVV10d (whose prediction is only slightly worse than the rVV10 one), thus increasing the confidence in the reliability of this functional.
As far as the dipole moment is concerned, with rVV10d this quantity is not improved with respect to the experimental estimate (see Table 5), which can be rationalized as fol-  lows: as a consequence of applying the modified rVV10d functional, the electronic charge is more widespread around the S atoms, which leads to a better agreement with experimental ED findings and to an improved estimate of the ILBE. However this increased charge delocalization does not occur exclusively in the interlayer region, but also, for instance, in the intralayer region, particularly along lines connecting neighboring S atoms. Therefore, the net delocalization effect (see Figure 5) is approximately isotropic so that the S dipole moment is not appreciably changed from that estimated by the original rVV10 functional.
We have also performed additional calculations on TiS2 by adopting vdW-corrected DFT functionals different from rVV10. Interestingly, while the inclusion of the new pseudopotential involving d orbitals leads to improvements similar to those observed with rVV10 in the description of the ED profiles (with a small charge accumulation in the interlayer region) and in the ILBE estimate, using vdW-corrected DFT functionals based on non-local functionals, such as VDW-DF-cx, 42 the same is not the case when vdW-corrected DFT functionals are adopted, which are more empirical in character (and not directly dependent on the detailed electron density distribution), such as the DFT-D3 approach, 43 which indicates that the new pseudopotential is indeed better suited to describe the subtle charge deformations occurring in this system. Finally, to further assess the validity of this novel approach, we have applied the rVV10d functional to other three different TMDs, whose chalcogen is always the S atom: MoS 2 , TaS 2 , and HfS 2 . The results for the equilibrium lattice parameter c, as well as the ILBE and ILBE* (see Table 6 ), are very promising. The major flexibility of the new pseudopotential, embedded in rVV10d scheme, dramatically improves the energetic description of S-based TMDs, maintaining an appreciable quality of the estimate of the longitudinal lattice parameter c: in fact, although the error relative to c is ∼ 2 % bigger than with rVV10, the mean absolute relative error (MARE) of the ILBE is reduced by about five times; this means that the pronounced overbinding tendency of rVV10 is almost totally eliminated.

Conclusions
We have presented the results of a first principles, DFT investigation, of the energetic and electronic properties of TiS 2 , that is a system where a puzzling discrepancy, between the distribution of the electron density obtained by X-ray diffraction data and that computed by state-of-the DFT schemes, has been recently reported, mainly concerning the interlayer region between S atoms belonging to adjacent layers. We basically confirm this observation, using the LDA, the (semilocal GGA) PBE, the (vdW-corrected) rVV10 DFT functionals, and also a modified rVV10 functional (rVV10(b=10.0)) where the b parameter has been increased so that short/medium-range vdW interactions are artificially damped. If instead the b parameter is considerably reduced to b=1.0 (rVV10(b=1.0)) significant changes occur: a stronger ED deformation in the interlayer region is accompanied by a dramatic (and physically unrealistic) overestimate of the CE and ILBE terms. Substantially decreasing (by a 4.5 factor) the other (C) rVV10 parameter and further adjusting the b parameter (rVV10(b=1.1,C=0.00207) functional) only slightly reduces the dramatic overbinding of rVV10(b=1.0) scheme. Interestingly, a more substantial improvement, in the direction of making the theoretical description in closer agreement with the experiment, is obtained by adopting for S a modified, more flexible, pseudopotential, involving d orbitals (rVV10d functional). Although the rVV10d estimate of the atomic S dipole moment still remains substantially larger (see Table 5) than the experimental one, a consistent improvement in the description of the ED deformation in the interlayer region and of the ILBE is achieved.
Probably, the residual discrepancy with the experiment could be only eliminated by an higher-level theoretical, first-principle approach, such as a Quantum Monte Carlo scheme, in principle able to include the whole set of correlation effects. Such an approach would be however extremely expensive computationally because very small ED deformations could only be accurately reproduced by very long simulations to get a sufficiently small statistical error associated to differential densities. In any case, although all fine details in the ED deformation map cannot be captured by rVV10d, nonetheless we suggest that this DFT functional indeed represents a reasonable compromise between accuracy and efficiency for improving the theoretical description of TiS 2 . Additional investigations indicate that similar improvements, particularly in the evaluation of energetic terms, can be obtained by DFT calculations based on the rVV10d functional, applied also to other TMD materials containing S atoms and characterized by the same structure of TiS 2 , such as TaS 2 , HfS 2 , and MoS 2 .

Supporting information 6 Intralayer bonding
The analysis of the Ti-S covalent bond, which characterize the intralayer interactions, is presented here. The electron density (ED) (panel a) ) and Laplacian (panel b) ) profiles computed along the Ti-S line, as well as the ED deformation (panel c) ) and the laplacian map (panel d) ) of the Ti-S intralayer plane, are reported. Since the results achieved with all the functional schemes employed are very close to each other, we show only the PBE results as an example ( Figure S1). In Table S1, the characterization of the two (3, -3) CPs and the (3, -1) CP is presented. As can be seen, very similar values are predicted by all the DFT functionals applied and these are also in good agreement with experimental data, with the exception of the LDA estimate of the Laplacian of (3, -3) CP #1, that exhibits a larger error and the Laplacian values at the (3, -1) CP, which are also substantially overestimated.

ED interlayer Laplacians
Here we present the Laplacian maps of the interlayer plane of interest, computed for all the functional employed (the LDA map is very similar to the ones shown).

S 2 molecule characterization
The calculations for the S 2 molecule, were carried out using the same approach adopted for the simulations of the TiS 2 system (see section 2 in the main text); the only differences were the chosen supercell (a cube with sides of 13.2 A) and the sampling of the BZ limited to the Γ point. Calculations were spin-polarized and the S 2 vibrational frequency was computed using the PHonon program of QE.

Acknowledgements
We acknowledge funding from Fondazione Cariparo, Progetti di Eccellenza 2017, relative to the project: "Engineering van der Waals Interactions: Innovative paradigm for the control of Nanoscale Phenomena".