Van der Waals effects on grazing incidence fast atom diffraction for H/LiF(001)

We theoretically address grazing incidence fast atom diffraction (GIFAD) for H atoms impinging on a LiF(001) surface. Our model combines a description of the H-LiF(001) interaction obtained from Density Functional Theory calculations with a semi-quantum treatment of the dynamics. We analyze simulated diffraction patterns in terms of the incidence channel, the impact energy associated with the motion normal to the surface, and the relevance of Van der Waals (VdW) interactions. We then contrast our simulations with experimental patterns for different incidence conditions. Our most important finding is that, for normal energies lower than 0.5 eV and incidence along the<100>channel, the inclusion of Van der Waals interactions in our potential energy surface yields a greatly improved accord between simulations and experiments. This agreement strongly suggests a non-negligible role of Van der Waals interactions in H/LiF(001) GIFAD in the low-to-intermediate normal energy regime.


I. INTRODUCTION
The extraordinary sensitivity of grazing incidence fast atom diffraction (GIFAD or FAD) has turned this phenomenon into a powerful surface analysis technique, positioning it also as a most useful tool for testing potential energy surfaces (PESs) [1][2][3][4].
When atomic projectiles in the keV energy range grazingly impinge on a crystal surface along a low-index crystallographic direction, scattering proceeds under axial surface channeling conditions [5]. The fast motion along the channel is, on a first approach, sensitive only to the periodic-PES average in this direction. The associated energy E is thus essentially conserved, and motions parallel and perpendicular to the channel get decoupled from each other. The scattering process can then be projected into the plane normal to the channeling direction, where motion proceeds with an energy E ⊥ in a hyperthermal up to eV energy regime and a De Broglie wavelength of the order of the interatomic spacing. Diffraction patterns arise due to the interference produced by the periodic array of channels, modulated by that originated within a given channel [6,7]. GIFAD was first reported for light projectiles impinging on wide band-gap insulating surfaces [8,9], but has since been observed for a variety of systems including semiconductors [10,11], metals [12,13], adsorbate-covered metal surfaces [14], ultrathin films [15] and organic molecules on metal substrates [16].
From the theoretical standpoint, the extreme sensitivity of GIFAD poses a challenge for achieving an appropriate description of this phenomenon. The construction of a projectile-surface potential which includes the key features of the interaction and a scattering dynamics representation which retains the quantum character of the process are necessary ingredients for attaining good accord with experimental diffraction patterns.
Projectile-surface interaction potentials for GIFAD simulations are usually built from Density Functional Theory (DFT) calculations. Within a standard DFT approach, the exchange-correlation energy is described by means of local or semi-local functionals, thus leaving long-range dispersion forces aside. This level of approximation is appropriate for high-enough E ⊥ GIFAD, but it may not be sufficient for the low-E ⊥ regime. In this latter case, the projectile is scattered farther from the surface, in lower electron-density regions where Van der Waals (VdW) interactions should not be neglected.
Studies on the dynamic aspects of VdW interactions are still scarce [11,[17][18][19][20]. In GIFAD literature, Zugarramurdi et al. [4] modeled the interaction of He atoms with a graphene layer on 6H-SiC(0001) by means of a pairwise additive Lennard-Jones potential fitted from Helium Atom Scattering (HAS) data; Debiossac et al. [11] considered ad-hoc corrections to the attractive part of a DFT potential for He atoms impinging on the β 2 (2 × 4) reconstructed GaAs(001) surface; and Schüller et al. [20] included VdW contributions into a DFT interaction potential by means of Grimme's semi-empirical approach [21]. The latter authors addressed the role of VdW contributions for He atoms impinging on the insulating MgO(001) surface, reporting no significant effect due to VdW either on the rumpling, the interaction potential V (z), or the corrugation of the equipotential curves in the normal energy range relevant for GIFAD.
In this article, we report on VdW effects on GIFAD for H/LiF(001), one of the systems where this phenomenon was initially observed [8,9]. Despite the considerable amount of experimental data for this collision system [5,8,[22][23][24][25], theoretical research is scarce [26,27], and ab initio simulations that satisfactorily reproduce all the experiments are still lacking.
We describe the elastic scattering of H atoms off the LiF(001) surface within the Surface-Initial Value Repre-sentation (SIVR) approximation [28], which is a semiquantum method that affords a clear representation of the main physical mechanisms in terms of classical trajectories through the Feynman path integral formulation of quantum mechanics [29]. The SIVR method includes an approximate representation of classically forbidden transitions on the dark side of the rainbow angle, providing an appropriate description of GIFAD patterns along the whole angular range [28,30]. Another noteworthy point is that, with a relatively low computational cost, the SIVR approach takes into account the threedimensionality of the PES, without averaging it along the incidence direction. These features make SIVR a most attractive alternative to quantum wave packet propagation methods.
We show that a description of the H-Surface interaction based on PAW pseudopotentials and the semi-local GGA-PBE functional fails in reproducing the diffraction patterns for incidence along the 100 channel in the E ⊥ 0.5 eV energy range. However, upon inclusion of VdW interactions in the PES through the semi-empirical approach by Grimme [21], we achieved a much improved, almost quantitative accord between our simulations and the experiments. To our knowledge, this agreement provides the first indication of the non-negligible role of VdW interactions in GIFAD, obtained from a non adhoc potential.
The paper is organized as follows: The theoretical models used to describe the quantum scattering and the projectile-surface interaction are summarized in Sec. II. Results for incidence along the 110 and 100 channels are presented and discussed in Sec. III, with the focus on the influence of the VdW contribution. In Sec. IV we outline our conclusions.

II. THEORETICAL MODEL
Our theoretical description of GIFAD combines a semiquantum representation of the scattering process with an accurate projectile-surface interaction potential. They are both summarized in the following subsections.

A. Scattering process
We treat the scattering dynamics of H atoms grazingly colliding with the LiF(001) surface by means of the SIVR approximation [28,30], expressing all quantities in atomic units (a.u.). Within this approach the transition amplitude per unit of surface area S reads [28] A K o respectively being the starting (t = 0) position and momentum of the projectile. In Eq. (1), the starting position is expressed as is the component parallel to the surface plane, z is the normal to the surface and Z o is a reference distance for which the projectile is hardly affected by the surface interaction. In turn, the starting momentum where V P S represents the projectile-surface interaction, (3) is the SIVR phase at the time t, with − → P t = m P d − → R t /dt the classical projectile momentum. In Eq. (2) the Maslov function [31] is a Jacobian factor (a determinant) evaluated along the classical trajectory − → R t , with |J M (t)| the modulus of J M (t) and ν t an integer number that accounts for the sign of J M (t) at a given time t, increasing by 1 every time that J M (t) changes its sign along the trajectory.
The SIVR differential probability for elastic scattering with final momentum [28], with θ f the final polar angle, measured with respect to the surface, and ϕ f the azimuthal angle, measured with respect to the channel direction (see Fig. 1). In the present work, the transition amplitude A (SIV R) if is obtained from Eq. (1) by employing the MonteCarlo technique with more than 4 × 10 5 points in the − → R os and Ω o integrals. In such integrations, the random − → R os values are derived from a Gaussian distribution covering an area S equal to 2 or 3 reduced unit cells, while the Ω o values are obtained from a Gaussian distribution encompassing an angular region determined by the Heisenberg uncertainty relation. In this aspect, it should be mentioned that the − → R os and Ω o distributions are in principle defined by the profile of a coherent wave packet associated with the impinging particle, which depends on the collimation of the incident beam [30,32,33]. Here we have considered standard sizes of the − → R os and Ω o distributions because reported experimental data [5,8,[22][23][24][25] lack information about collimating parameters. In addition, the starting normal distance was chosen as Z o = 1.4 a (a is the lattice constant), to ensure a negligible projectile-surface interaction.

B. Projectile-surface potential
To obtain the H-LiF(001) potential V P S required in Eqs. (2) and (3), we make use of DFT, as implemented in the QUANTUM ESPRESSO code [34], to calculate the system's energy for a grid of (X i , Y i , Z i ) positions of the H atom over a relaxed LiF(001) surface. The grid is three-dimensional (3D) and is built out of a selection of 6 high-symmetry (X i , Y i ) configurations and 62 Z i values (Z = 0 falls on the topmost F layer). We then apply an interpolation technique, which combines the Corrugation Reducing Procedure (CRP) [35] with the cubic spline method, to obtain the potential energy for an arbitrary (X, Y, Z) position.
For the DFT calculations, we use projector augmentedwave (PAW) pseudopotentials [36,37] to describe the electron-core interaction, while for the exchangecorrelation functional we consider two different models: (a) the generalized gradient approximation (GGA), with the Perdew-Burke-Ernzerhof (PBE) functional [38] (henceforth referred to as PAW-PBE), or (b) the DFT-D2 approach of Grimme [21,39], which introduces a semiempirical correction to the GGA functional to account for long-range VdW interactions (henceforth PAW-PBE-VdW).
For both PAW-PBE and PAW-PBE-VdW models, we choose relevant parameters of the DFT calculations so that ab-initio energies differ from the converged result in less than 5 meV. The energy cutoff in the plane-wave expansion is 100 Ryd for the wave functions and 1000 Ryd for the charge density and potential; fractional occupancies are determined through a gaussian broadening approach with σ = 10 −6 Ryd.; and the Brioullin-zone integration is performed with a 4 × 4 × 1 Monkhorst-Pack grid of special k-points with an offset. The LiF lattice constant is a = 4.066Å for the PAW-PBE case and a = 4.063Å for the PAW-PBE-VdW case, both slightly higher than the experimental value of 4.02Å [40].
We represent the LiF(001) surface by means of the supercell-slab scheme. The supercell consists of a 2 × 2 surface cell, a five-layer slab and a vacuum distance of 5 a. The surface equilibrium geometry is reached by relaxing the two topmost LiF(001) planes from their bulk posi-tions. F and Li atoms initially in the same plane relax differently and the resulting surface geometry presents a rumpling, which we define as the distance between relaxed F and Li planes. For the topmost F and Li planes, we get rumplings of +0.070Å (PAW-PBE) and +0.088Å (PAW-PBE-VdW), with F atoms moving outward and Li atoms moving inward. These values are consistent with LEED experiments which yield a rumpling of 0.036 ± 0.1 A [41] and, particularly, our PAW-PBE result compares very well with Vogt's 0.068Å [42], also obtained from a GGA calculation. The relaxed surface is thereafter kept frozen both for the energy grid calculations, and during the scattering process.

III. RESULTS AND DISCUSSION
In this section we first discuss the general features of the PESs used in this work; then we proceed to a systematic theoretical study of GIFAD patterns for H/LiF(001) in terms of the incidence channel and the normal energy E ⊥ = E tot sin 2 θ i (E tot the total energy and θ i the incidence angle relative to the surface plane. See Fig. 1); and finally, we compare our theoretical simulations with experimental GIFAD distributions available in the literature [5,8,[22][23][24][25]27], in order to probe the PAW-PBE and PAW-PBE-VdW projectile-surface potentials.
The geometry of GIFAD for H/LiF(001) as well as the channeling directions 110 and 100 are illustrated in Fig. 1.

A. Analysis of the PES
The PES profile and corrugation near the reflection region are central in determining GIFAD patterns. We stress again that our PES is 3D and no dimension reduction is made during the dynamics. However, the fast motion of the projectile along the channel is in fact mainly sensitive to the average of V P S in this direction and thus we will discuss the PES role in these terms.
For our PAW-PBE and PAW-PBE-VdW PESs, in Figs. 2a and 2b we consider the energy averages respectively along the 110 and 100 channels, and depict equipotential contours across them. Across the 110 channel the equipotential curves have local maxima both at the border and at the middle of the channel, respectively corresponding to the rows of F and Li ions. Across the 100 channel, the equipotential curves have only one maximum at the border of the channel, corresponding to the F-Li rows. In Table I, we consider equipotential curves in the 0.3-0.9 eV range and show their Zrange (Z min , Z max ) and corrugation ∆Z = Z max − Z min for both PESs and channels. Note that, for projectiles that run parallel to the channel without suffering any azimuthal deflection (i.e. ϕ f = 0), these Z-ranges determine the reflection region. From Fig. 2 and Tab. I we obtain that i) the corrugation across the 100 channel is much higher than that across the 110 channel; ii) the inclusion of VdW interactions results in a less repulsive and higher corrugated PES, these effects growing stronger for lower energies; and iii) VdW interactions seem to play a more important role for the 100 channel than for the 110 one.
While the feature described in i) has already been observed both experimentally [22] and theoretically [26], and the one in ii) is related to the essentially attractive character of the VdW interaction, that dominates at long distances from the surface, the feature given in iii) deserves further discussion. Along the 110 channel, in Ref. [43] it was found that for He/LiF(001) the presence of cationic and anionic rows (see Fig. 1) induces polarization within the projectile, resulting in marked effects on GIFAD patterns. Polarization is implicitly included in a standard DFT calculation and therefore in our PAW-PBE PES. It mainly affects the projectile-surface interaction in the medium to large distance range, where it may compete with VdW interactions, probably overshadowing them. In contrast, along the 100 channel, neutrally charged rows of alternating F and Li ions (see Fig. 1) result in no projectile polarization effects, and hence in more visible VdW contributions.

B. Simulated GIFAD patterns
Simulated GIFAD patterns obtained for both incidence channels and PESs are displayed in Fig. 3, as a function of the final azimuthal angle ϕ f , considering a selection of E ⊥ values ranging from 0.4 to 0.8 eV.
GIFAD angular distributions present peaks associated with Bragg diffraction, which are produced by interference among equivalent trajectories whose starting positions − → R os lie on different parallel channels [5,6,44]. These peaks are situated at azimuthal angles that verify sin ϕ f = nλ/d, where λ = 2π/K i is the de Broglie wavelength of the incident atom, d is the width of the channel (d = δ for 110 and d = a/2 for 100 ), and n is an integer number that determines the Bragg order. Hence the positions of Bragg peaks, which depend on the total energy through λ, provide crystallographic information only. Their intensities however are determined by a unit-cell form factor that is originated from interference among trajectories with − → R os within the same reduced unit-cell. This unit-cell form factor acts as an oscillatory envelope function that can reduce or even suppress the contribution of a given Bragg order [28,44]. In GIFAD the intensities of the Bragg peaks are extremely sensitive to the shape of the PES across the incidence channel [2,45,46], being in most cases completely governed by E ⊥ [7,47]. For the present collision system we have in fact verified that the spectra of Fig. 3, obtained for a fixed incidence energy E tot = 1 keV, does not change appreciably upon setting E tot = 3 keV, while keeping the E ⊥ values unchanged. Finally, the number of observed Bragg orders is fully determined by the unit-cell form factor, which depends on E ⊥ [27].
In Fig. 3 we observe that the PESs features discussed in subsection III-A directly affect the intensity profiles. Increasing E ⊥ along a given channel results in reflection at more corrugated regions and thus in wider diffraction patterns while, for fixed E ⊥ , the lower corrugation across the 110 channel relative to that across the 100 one (See Tab. I), results in a narrower pattern for the former. Regarding the VdW contribution, it has a more visible role along the 100 channel than along the 110 one, which is expected due to the difference in the projectile polarization contribution. Also, for the 110 channel VdW effects grow stronger for lower E ⊥ values, that is, when projectiles probe regions farther from the surface. But for the 100 channel the most marked difference between PAW-PBE and PAW-PBE-VdW patterns is observed for E ⊥ = 0.6 eV. Noteworthy, we find that, for E ⊥ = 0.4 eV, VdW has a marked effect for both channels.

C. Comparison to Experiments
In this subsection we present simulated GIFAD patterns for different incidence conditions, corresponding to available experimental data. We compare to the experiments the results derived from both the PAW-PBE and the PAW-PBE-VdW PESs, with the aim of assessing their performance in terms of incidence channel and normal energy. In Tab. II we enumerate the various experimental settings considered in this work, and assign a code name (channel|E ⊥ × 100) to each of them, with E ⊥ expressed in eV. We will hereafter use these code names to refer to each particular setting.

Incidence along the 110 channel
We first present and discuss the simulations and experiments for incidence along the 110 channel, where weaker VdW effects are predicted.
In Fig. 4, the simulated intensity distributions for E ⊥ = 0.56 eV (case (110|56)) are contrasted with the experimental projected intensity profile reported by Rousseau et al. [22]. Theoretical distributions are here normalized to the central maximum (i.e., at ϕ f = 0), while the experimental one is normalized to the secondorder Bragg peak corresponding to the PAW-PBE-VdW potential. For (110|56), PAW-PBE and PAW-PBE-VdW potentials produce very similar patterns which satisfactorily reproduce the overall characteristics of the experimental distribution, showing a very intense central peak sided by two much lower maxima associated with the n = ±2 orders. The small differences between our PAW-PBE and PAW-PBE-VdW theoretical distributions can be explained from the values given in Tab. I, where the latter shows a slightly higher corrugation relative to the former, which results in increased intensities of its nevertheless low n = ±1, ±2 peaks. At the positions corresponding to the n = ±1 peaks, the intensity is higher in the experiment than in our calculations. However, this discrepancy might be attributed to experimental limitations that do not allow one to distinguish the different Bragg orders. At such angular positions the experimental profile seems to include contributions from the broad central maximum, which hinder a stringent comparison with the theoretical curves.
We continue on to address the higher normal energy (110|88) case, corresponding to the experiments for E ⊥ = 0.88 eV, reported by Schüller et al. [8]. In Fig. 5 we show two-dimensional angular distributions, as a function of the final polar and azimuthal angles (θ f , ϕ f ).
Since the θ f -length of GIFAD patterns is affected by the collimating conditions of the incident beam [30,33], not given in Ref. [8], the polar angle is here plotted in arbitrary units. Both the PAW-PBE and PAW-PBE-VdW calculations nicely reproduce the experimental pattern displaying five maxima with comparable intensities. The similarity of PAW-PBE and PAW-PBE-VdW GI-FAD distributions can be explained on inspecting the E ⊥ ∼ 0.88 eV data in Fig. 2a and Tab. I, which show that the corrugations and Z-ranges for PAW-PBE and PAW-PBE-VdW differ only slightly. The likeness of both simulated diffraction patterns can thus be traced to that of the averaged PESs near the reflection region. It is worth noting that, for this case, the H-Surface reflection distances are only a little larger than the HF molecule internuclear distance of 0.917Å.
For incidence along the 110 direction, we can then conclude that our calculations well reproduce the experiments, confirming that PAW-PBE performs quite satisfactorily while VdW plays a negligible role both in the PES and in the diffraction patterns for E ⊥ 0.5 eV. Lack of experiments for lower E ⊥ values prevents us however from exploring the interesting E ⊥ 0.4 eV regime, for which our theoretical study predicts a visible contribution of VdW interactions.

Incidence along the 100 channel
Following, we address the 100 channel which, according to our results from Subsection III-A, is more favorable for studying possible VdW effects, due to the absence of polarization.
In Fig. 6 we plot azimuthal angle spectra for the low normal energy case (100|29) (E ⊥ = 0.29 eV). Both our simulated patterns, as well as the experimental data, taken from Ref. [27], are normalized to the central peak.
Remarkably, PAW-PBE results in a rather poor agreement with the experimental pattern, which strikingly contrasts with the almost quantitative accord achieved by PAW-PBE-VdW. A key ingredient of this better performance is the increased corrugation of PAW-PBE-VdW near the reflection region which, as shown in Fig. 2b and Tab. I, approximately doubles that of PAW-PBE.
In an analogous fashion, the simulated azimuthal angle spectra for cases (100|45) and (100|51) (respectively E ⊥ = 0.45 eV and E ⊥ = 0.51 eV), are shown in Figs. 7 and 8, together with the corresponding experiments, taken also from Ref. [27]. For both cases, PAW-PBE performs rather unsatisfactorily, overestimating the n = ±1 Bragg orders, while PAW-PBE-VdW provides a much better description, reproducing the similar intensities of the n = 0 and n = ±1 orders and the much lower n = ±2 peaks.
Cases (100|29), (100|45) and (100|51) share the common features of unsatisfactory PAW-PBE patterns and very good descriptions provided by PAW-PBE-VdW along the complete angular range. In this low-tointermediate-E ⊥ regime the H atom scattering dynamics is restricted to regions where the semi-local PBE functional appears to provide an inadequate description of the H-surface interaction, of which PAW-PBE-VdW seems to give a much better representation. Concerning a possible influence of the functional choice on our PES (and hence on the simulated GIFAD patterns) we can mention that Muzas et al. [27] have recently reported GIFAD simulations for the low-to-intermediate E ⊥ cases along 100 . These authors modeled the projectile-surface interaction with a PAW-PW91 PES, where the PW91 exchangecorrelation functional was used, instead of the PBE. This PAW-PW91 PES was built with precision criteria comparable to our PAW-PBE PES and, noteworthy, leads to GIFAD patterns very similar to our PAW-PBE results. This fact contributes to support our claim that the reason behind PAW-PBE poor performance for this channel and normal energy regime is not our functional choice, but the neglect of VdW interactions.
In Fig. 9 we show two-dimensional angular distributions in (θ f , ϕ f ) for the (100|53) case (E ⊥ = 0.53 eV), corresponding to experiments by Winter et al. [5]. In this case our PAW-PBE calculation performs remarkably well, correctly reproducing both the outer low-intensity n = ±2 peaks as well as the higher intensity of orders n = ±1 relative to the central n = 0 peak. On the contrary, the pattern obtained with PAW-PBE-VdW is not as good, yielding an apparent overestimation of both the n = 0 and n = ±2 peaks.
On increasing the E ⊥ value, the VdW contribution is expected to gradually fade out. Therefore, it is intriguing that our PAW-PBE-VdW simulation yields worsebehaved patterns than PAW-PBE ones. This failure is also observed in Figs. 10 and 11 for the higher energies E ⊥ = 0.55 eV and E ⊥ = 0.64 eV, respectively, corresponding to the cases (100|55) and (100|64), and might probably be traced to a switch-off region issue. In a PES which includes VdW corrections, the relative importance of semi-local functionals and VdW terms varies with the H-Surface distance. The matching region, where VdW corrections are smoothly switched off, is where spurious effects of the VdW approach may arise. Further increasing the E ⊥ value results in a convergence of the PAW-PBE-VdW equipotential curves to the PAW-PBE ones, as shown in Fig. 2, and thus we expect both GIFAD simulations to eventually converge to a common pattern as well, as was the case for the 110 .
The very different experimental patterns reported for the (100|53) and (100|51) cases illustrate how a small variation in the normal energy may produce a substantial modification of the GIFAD spectrum. The unexpected patterns PAW-PBE yields for E ⊥ = 0.64 eV, shown in Fig. 11, are probably another example of this. Noteworthy, while PAW-PBE for 0.64 eV does not satisfactorily reproduce the experimental pattern by Busch et al. [24], we do find a nice accord between this experiment and our PAW-PBE simulation for 0.6 eV, shown in Fig. 3.
Summing up, for incidence along the 100 channel we find that, on the one hand, the agreement of PAW-PBE-VdW patterns with experiments for cases (100|29), (100|45) and (100|51) is almost quantitative and strongly suggests a non-negligible role of VdW interactions for low-to-intermediate-E ⊥ cases in H/LiF(001) GIFAD. On the other hand, PAW-PBE-VdW performs worse than PAW-PBE for cases (100|53), (100|55) and (100|64), a feature that might probably be related to spurious effects of the VdW approach arising when both the functional and VdW interactions have non-negligible contributions to the energy. These most interesting results strongly call for more experiments and theoretical work in the E ⊥ > 0.5 eV energy range.

IV. CONCLUSIONS
In this article, we have studied GIFAD on H/LiF(001), comparing the diffraction patterns obtained with a PAW-PBE PES with those obtained with a PAW-PBE-VdW PES, where VdW has been included following Grimme's approach. We have theoretically investigated the relevance of VdW interactions along the 110 and 100 channels in the 0.4-0.8 eV E ⊥ -range and have predicted a marked influence of VdW corrections for both channels in the low-E ⊥ (E ⊥ ≤ 0.4 eV) regime. Also we have found that VdW corrections affect more strongly the 100 than the 110 channel, due to the presence of polarization effects for the latter direction. These effects compete with VdW in the intermediate to large distance region.
From the comparison with available experiments, we have found that PAW-PBE gives an adequate description of GIFAD along the 110 channel, for E ⊥ > 0.55 eV cases and, along the 100 channel, for E ⊥ = 0.53 and 0.55 eV. GIFAD simulations along the 100 channel for E ⊥ ≤ 0.51 eV cases however result in a poor agreement with experiments unless VdW interactions are considered, through the PAW-PBE-VdW PES, in which case our simulated patterns remarkably achieve quantitative agreement with experiments. This certainly is the main finding of the present work and to our knowledge, it might be the first evidence of marked VdW effects in GIFAD, obtained with a DFT potential which includes VdW interactions in a non ad-hoc fashion, through Grimme's semiempirical approach [21].
Other worthmentioning points are: i) The prediction that, for incidence along the 110 channel, VdW contributions should become relevant for low impact energies (E ⊥ 0.4 eV), and ii) the quite intriguing patterns obtained for the 100 channel in the E ⊥ ∼ 0.5 − 0.65 eV range. More experiments for both these incidence conditions would be most desirable for a thorough analysis of the reliability of PES models and the influence of the VdW interaction.
Moreover, the present results open the way for many further studies on this topic. In particular, more experiments are important to study the performance of this semiempirical approach to VdW for E ⊥ ∼ 0.6 eV and incidence along the 100 channel. The high-E ⊥ regime is technologically appealing as, under it, GIFAD can reach topological resolution thus becoming a reciprocal space analog of a perfect tip AFM [48]. Low-E ⊥ as well as intermediate-E ⊥ GIFAD might prove just as relevant, providing a highly sensitive quality check for Potential Energy Surfaces, and thus contributing in the development of an accurate description of VdW interactions within DFT.    Code Name indicates how we will refer to the particular incidence setting in the text. ( * ) refers to cases for which the simulations were performed considering Etot = 1 keV, verifying the stability of the pattern obtained upon Etot variation (keeping E ⊥ fixed).