Nonadiabatic localization of H2 in the field of two external positive tip charges

J. Vac. Sci. Technol. A 39, 053206 (2021); https://doi.org/10.1116/6.0001138 39, 053206 © 2021 Author(s). Nonadiabatic localization of H2 in the field of two external positive tip charges Cite as: J. Vac. Sci. Technol. A 39, 053206 (2021); https://doi.org/10.1116/6.0001138 Submitted: 11 May 2021 • Accepted: 15 July 2021 • Published Online: 04 August 2021 W. Schattke, M. A. Van Hove and R. Díez Muiño COLLECTIONS Paper published as part of the special topic on Commemorating the Career of Charles S. Fadley ARTICLES YOU MAY BE INTERESTED IN X-ray spectroscopic identification of strain and structure-based resonances in a series of saturated carbon-cage molecules: Adamantane, twistane, octahedrane, and cubane Journal of Vacuum Science & Technology A 39, 053208 (2021); https:// doi.org/10.1116/6.0001150 Green CVD—Toward a sustainable philosophy for thin film deposition by chemical vapor deposition Journal of Vacuum Science & Technology A 39, 051001 (2021); https:// doi.org/10.1116/6.0001125 Epitaxial film growth by thermal laser evaporation Journal of Vacuum Science & Technology A 39, 053406 (2021); https:// doi.org/10.1116/6.0001177 Nonadiabatic localization of H 2 in the field of two external positive tip charges Cite as: J. Vac. Sci. Technol. A 39, 053206 (2021); doi: 10.1116/6.0001138 View Online Export Citation CrossMar k Submitted: 11 May 2021 · Accepted: 15 July 2021 · Published Online: 4 August 2021 W. Schattke, 1,2 M. A. Van Hove, 3 and R. Díez Muiño 1,4 AFFILIATIONS 1 Donostia International Physics Center DIPC, P. Manuel de Lardizabal 4, 20018 San Sebastián, Spain 2 Institut für Theoretische Physik und Astrophysik der Christian-Albrechts-Universität, Leibnizstraße 15, 24118 Kiel, Germany 3 Department of Physics and Institute of Computational and Theoretical Studies (ICTS), Hong Kong Baptist University, 224 Waterloo Road, Kowloon, Hong Kong 4 Centro de Física de Materiales CFM (CSIC-UPV/EHU)—Materials Physics Center MPC, P. Manuel de Lardizabal 5, 20018 San Sebastián, Spain Note: This paper is a part of the Special Collection Commemorating the Career of Charles S. Fadley. ABSTRACT For two external spherical tips with equal positive charges, the ground state of a hydrogen molecule is variationally determined within the quantum Monte Carlo scheme. For finiteness, the system is enclosed in a spherical container with randomly reflecting unstructured walls. The 12-dimensional system of two protons and two electrons is investigated ab initio without any adiabatic restrictions. We focus on the hydrogen center of mass (COM) distribution in continuation of previous work on ground state and vibration modes at a COM that was fixed in space. Our general purpose is to control the molecule’s position and orientation by external means, such as by two charged tips. To this end, knowledge is needed in order to find specific COM regions with a desired molecular orientation in the container at significant probability density of the molecule. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1116/6.0001138 I. INTRODUCTION There is continuous interest in the hydrogen molecule and in its technical applications in the nanoscale and subnanoscale domains. This involves the hydrogen nuclei as well as their electrons where the latter usually deserve a quantum mechanical consideration. The nuclei are because of their larger mass on the borderline to classical mechanics. Nevertheless, quantum mechanics influences their motion, too. For instance, clamping a hydrogen molecule as a bridge between tunneling tips can be a successful but tedious task. It is a first step toward external control of position and orientation of an object that is small enough to penetrate tiny openings and heavy enough to transport tiny vibrational energy quanta. An important research direction has already developed as supramolecular chemistry 1 from molecular switches to machines 2–6 where larger objects are manipulated by combining nanomechanics with nanochemistry. Nevertheless, the very basic compound chosen here may still offer interesting insights if considered from the full quantum mechanical ab initio level. Here, we address experiments such as investigating the vibrational energies of hydrogen bridging two tunneling tips 7–10 or measuring the hydrogen-induced forces when trapped in tunneling nanocavities 9,11–13 as well as calculations on the stability of molecular configurations involving a hydrogen molecule as a mediating glue. 14 Those findings ask for the investigation of simple systems that might lead to a basic understanding of more complicated situations. Based on complete quantum mechanical ab initio calculations, one has to include the proton degrees of freedom on the same level as the electrons without allowing an adiabatic assumption that more or less ignores the slow proton motion in comparison with the fast electrons. In terms of energy, the electrons’kinetic energy is larger by about a factor 2000 than that of the protons, suggesting a decoupling of the degrees of freedom or a perturbational treatment of the protons. Our previous investigation 15 on the lowest vibrational excitation of H2showed essential deviations from the adiabatic ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-1 ©Author(s)2021 approximation for a system with a fixed center of mass. Thus, we will avoid all adiabatic approximations for the nuclear COM degrees of freedom, even though they are used in the vast majority of calculations. Fundamental work on this molecule ranges from highprecision calculations 16,17 to high-precision measurement as, e.g., highly accurate vibrational determinations. 18 It extends as well into the challenge of examining the generally used Born–Oppenheimer assumption. 19–21 The latter is relevant to this work, because the light hydrogen molecule is highly sensitive to nonadiabatic corrections, i.e., the quantum mechanical nuclear motion coherently taken into account within the QMC electronic calculation. The nuclear degrees of freedom are connected with a spatial scale that is much smaller than the electronic scale by three orders of magnitude. Much larger in contrast is the scale of macroscopic objects such as the vacuum vessel containing the probe equipment and the H2gas. Our theory can handle the first two scales, while we must limit the macroscopic scale to a practical size. The confinement of simple molecules in a box as a model for a cavity has been widely investigated, which can serve to estimate how an enclosure will affect our system. 22–24 The results of such calculations may shed some light on the classical view of electrostatics and diffusion in the picture of a compact hydrogen molecule with an induced dipole subject to an electrostatic force. The properties of the isolated molecule, such as its ground state energy and polarizability, are well known by ab initio investigations, 21,25,26 where Tubman et al. utilize QMC methods beyond the adiabatic assumption too. The diffusion analogy to the Schrödinger equation is also an illustrative aspect of quantum mechanics. We expect to learn what can be deduced from such simple imagination or what replaces this picture in order to manage hydrogen in the nanoscale. In Sec. II, we describe some details of the calculational method. The results address the calculations of the ground state energy and wavefunction dependent on the distance between the tunneling tips, which is presented in Sec. III.Section IV summarizes the findings. II. SYSTEM AND METHOD A. Geometric arrangement The external configuration in which the H2molecule is embedded consists of a spherical container and two tips on the SX axis symmetrical to and near its center as depicted in Fig. 1. We first tried to treat the free H2system without external confinement, e.g., without container, expecting that the system could be limited by the presence of tips alone as an attractive region. But large energy fluctuations arose when optimizing the proton parameters and could not be handled in a satisfactory way. Thus, we decided to enclose the total system in a container which proved to be a controllable approach as it led to appreciably smaller fluctuations. The walls of the container act as a boundary condition that sets the outside wavefunction to zero. This procedure corresponds to continuing the protons’random walk of the Monte Carlo process after reaching the wall with subsequent random steps until a step into the interior occurs, which is called here “diffuse reflection.”Thus, we disregard the kinetic energy contribution which would necessarily arise with a finite boundary potential. Eventually, we checked this assumption which yields a controllable error in the limit of large cut-off radius as expected by the literature. 24 Both tips shall carry equal positive charges where we assume two equal extended spherical bodies that are homogeneously charged spheres with radius fixed here to the radius 2.55 a.u., say, of a spherically symmetric state of a tip molecule not further specified. In so doing, we simulate the finite extension of a charged tip and simultaneously avoid the singular point charge behavior. The tip-tip distance Dis a system parameter that is varied over a few atomic units in the regime of physisorption. The tip charge is fixed at Q¼þ0:5 a.u. with only a short nonoptimized extrapolation shown at the end. B. Calculation scheme and optimization parameters The variational quantum Monte Carlo calculations are performed by extending the scheme previously presented. 15 We use computer clusters and random number seeds that are distributed over 24–48 parallel CPUs with run lengths of 2 109Monte Carlo steps, which take about 12 h and use 100 GB of memory. Several sweeps over the more than 50 optimization parameters are manually controlled by inspecting their importance and sensitivity to reduce their number during the course of optimization sweeps. Regarding the sampling, we obtain in the final sweeps an acceptance ratio of the proton moves between 50% and 70% which usually is considered as reasonable. Nonergodic runs were FIG. 1. Setup for the H2center of mass (COM) configuration plotted in the (SX,SY) plane of COM coordinates (SX,SY,SZ). An outer spherical container functions as a system boundary. Two positively charged spherical tips are external tools to manage the probability distribution of the H2molecule. The objects are drawn in the actual scale given by the axes, e.g., 13.1 a.u. for the tip-tip distance Dand 20 a.u. for the radius of the container. The angles between the H2 axis and the three COM coordinates (SX,SY,SZ) are denoted by (ΘX,ΘY,ΘZ) and sketched in the drawing. ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-2 ©Author(s)2021 occasionally obtained when the parameter values force the MC run into very negative potential-energy regions (small electron-tip distance) trapping the run with an acceptance down to 0% to a single configuration at a single energy and variance. We had to adjust the proton maximum step width in the beginning of these investigations which also led to similar sampling problems. We are confident about a sufficient control of the MC runs by energy and variance and other expectation values such as hydrogen bond length and various densities. In the optimization procedure approximating the true ground state, we rely on the variance of the local energy which was suppressed to a similar low value as in the former calculations with pinned COM. The possibility of variationally missing the true minimum because of an incomplete Hilbert space cannot be ruled out, of course, and is present in most published calculations. The main extension with respect to our original work 15 is associated with the three COM degrees of freedom. The Monte Carlo step width for a move of the COM position is found to be reduced to 1/8 of that of the electrons in order to produce stable Monte Carlo (MC) runs over the full optimization procedure, especially avoiding the MC runs getting trapped in very low potential points. In addition, a direction scaling is employed in the step width by increasing it along the tip-tip axis by a factor of 10 for the protons’COM move and a factor of 3 for the protons’relative move. The wavefunction was extended by Gaussian type functions with different decay lengths parallel and perpendicular to the tip-tip direction. Additionally, regarding the influence caused by the tip positions double-zeta Gaussians with centers displaced from the origin had to be introduced. All decay lengths and centers are optimization parameters. Further, Jastrow-type functions have been added for the protons in order to account for stronger fluctuationsnearthetipcharges. The wavefunction is written as 15,27 ϒ(ri,Rj)¼Ψel(ri,Rj;αk)Φnu(Rj;ζl), i,j¼1, 2, (1) with an increased number of optimization parameters αand ζas follows below. It yields the energy E¼hϒj ^ Hjϒi¼ðϒ*(ri,Rj)Hϒ(ri,Rj)d6rd6R, (2) for an ab initio Hamiltonian ^ Hcomprising the 12-dimensional position space of two protons and two electrons in the field of two electrostatic charges enclosed in a finite spherical volume. This notation is fully general even though we factorize the wavefunction into a mixed electronic-nucleonic part Ψel(ri,Rj;αk) containing the electronic variables with some contributions from the proton variables and a nucleon-only part Φnu(Rj;ζl). The former part is mainly identical with the previous work, 15 while the latter has some important extensions that refer to the center of mass degrees of freedom not contained in the former paper with the COM fixed at the origin, i.e., at the center between both tips. Thus, the full ground state ansatz consists of a singlet one-component determinant for both electrons combined with a Jastrow factor and contains (1) two-body terms compensating the electron-electron, electronproton, and proton-proton Coulomb singularities; (2) one-body terms compensating the electron-tip and proton-tip Coulomb singularities, all written within the Jastrow exponent, and additionally; (3) two general geminal sums 28 for orbital-orbital electron coupling, one of which replaces a simple product for the determinant while the other in the Jastrow exponent reflects electron-electron correlation. In more detail, the ground state nuclear wavefunction is given by Gaussian distributions and Jastrow-type functions, written here partially in relative R:¼R2R1¼(X,Y,Z) coordinates, partially using directly the single proton coordinates Ri¼(Xi,Yi,Zi), i¼1, 2. The Gaussian part describing the relative motion is explicitly Φ(0) nu (R1,R2)¼exp ζ0ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (X2þY2þZ2) pR0  2 ζ0x(jXjR0x)2ζ0yz ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (Y2þZ2) pR0yz  2: (3) The three parameters R0,R0x,andR0yz are used as centers of the gaussian ansatz for the spherical, axial, and planar radial dimensions, respectively, instead of a single hydrogen bond length, and are associated with three corresponding parameters for the widths. They represent a variational freedom of in total six nuclear parameters which have to be optimized, in addition to the electronic parameters for the total energy minimum. The above variational functions reflect the actual symmetry. The nuclear positions are dynamically influenced by the electrons’positions and the tips’force field. In view of the latter, we need additional functions for their COM degrees of freedom. A gaussian for the nuclear center of mass coordinates S:¼(R1þR2)=2¼(SX,SY,SZ) reads as Φ(COM) nu (R1,R2)¼exp ζCOM,x(jSXjRCOM,x)2ζCOM,yz(S2 YþS2 Z) hi , (4) adding another set of three parameters where RCOM,xis related but not necessarily equal to the position of the tips. However, we think it is useful to expand the function set by superimposing a second gaussian with center shifted along the SXaxis in order to take special care of the tips. With a linear weight cξ, this introduces four further parameters so that Φ(COM,ξ) nu (R1,R2)¼exp ζCOM,ξ(jSXjRCOM,ξ)2ζCOM,ρ(S2 YþS2 Z) hi : (5) We considered only the above double-zeta gaussian with axial symmetry and no angle dependence. Also, an extra term for an explicit polarization by differing COM positions of electrons and protons did not improve the optimization. The polarization of the electron-proton charge cloud enters the wavefunction through the ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-3 ©Author(s)2021 double-zeta orbitals of dx2and dyz type involved in the electronic part which seems to yield sufficient accuracy. A kind of proton-tip Jastrow factor was introduced to deal with the Coulomb field of the positive tip charge, especially for close tips and small tip curvature. Led by the electron-nucleus functional form of the Jastrow factor, we analogously write for this contribution Φ(COM,J) nu (R1,R2)¼exp C FX k¼1,2 X i¼l,r 1 1þFjRkTij "# , (6) with Tias left and right tip positions. Two further optimization parameters C and F appear here. Altogether, we obtain Φnu ¼Φ(0) nu Φ(COM) nu þcξΦ(COM,ξ) nu  Φ(COM,J) nu :(7) III. RESULTS The optimization toward minimum energy and variance with precedence of energy was carried out for several tip-tip distances D by a series of MC runs over different parameters as well as their values. Of course, the variance had to be kept below a reasonable bound given by the fluctuations with parameter-value change. Larger variance levels usually were connected with extremely low potential values unbalanced by respective kinetic energy, indicating either still unadapted parameters or unsuitable functional forms as a whole. Misfit of the parameter values was usually smoothed in the course of subsequent sweeps, while an unsuitable functional dependence required enlargement of the wavefunction set. However, single statistical events outside expectation naturally always occur and are not a matter of concern. Besides energy and variance, other quantum mechanical observables are also sampled, namely, electron density, nuclear COM density, or the angle between the molecule axis and one of the coordinates (x,y,z) as a function of position. Figure 2 shows the COM distribution for an energy close to optimization at several tip-tip distances, namely, D¼7:1, 9.1, 13.1 a.u. In the case of the largest distance, the distribution seems to follow the lines of the electric field between two positive charges; these emerge in the (SX,SY) plane radially from both charge centers, run downhill in potential, meet each other at a potential’s saddle point at the origin, and run further downhill with denser flow along a channel around the perpendicular bisector. The hydrogen’s probability concentrates along a narrow middle disk there decreasing from the center toward the container wall, however, with a descent much smaller than the steep descent across the disk. In the large remaining part outside the disk, the density almost vanishes. The situation appears quite different for smaller tip-tip distances. The closely spaced positive tip charges seem to repel the molecule into the outside, primarily behind the two tips. Such a behavior would correspond to a diffusion of the hydrogen into the outer space because not enough space is left for the density between the tips. Of course, this is not an explanation but at most a key to remember. A more quantitative picture of the density landscape is obtained by plotting the density vs SXalong horizontal cuts averaged over stripes as shown in Fig. 3 referring to two graphs of Fig. 2. Considering the angle between the H2axis and the COM coordinate axes SX,SY,Sz, its distribution is rather simple. An almost homogeneous average angle, see Fig. 4, appears with the hydrogen axis oriented parallel to the tip-tip connection over the full available space. The spatial distribution of the angle shows in all depicted regions the same qualitative behavior, namely, a maximum at 90for the angle with reference to the SYor SZaxis and close to 0with reference to the SXaxis, keeping in mind the sin(/(H2,SX)) factor for the phase volume of the solid angle at the center S¼0. Several sections of the SX,ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S2 YþS2 Z p  -plane demonstrate in connection with the density representation of Fig. 2,a homogeneous unique orientation of the angle’s quantum mechanical expectation value all over the accessible space. FIG. 2. Probability density distribution of the H2molecule plotted on the (SX,SY)-plane scaled in color by tenths of maximum density of ρmax 2:86 1051/(a.u.)3. The color sections are defined by (n1)=10 ,ρ=ρmax n=10 for n3 with lowest interval 1=100 ,ρ=ρmax 2=10 and keys showing the values of n. The positive charge Q¼0:5 a.u. of the tips extends over a radius of 2:55 a.u. as sketched in the figure. ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-4 ©Author(s)2021 The direction of the molecule axis was at first surprising as we initially supposed that it might spatially fluctuate more and correspond to the orientation with the center of mass pinned at the origin. Nevertheless, the kinetic energy tries to smooth the spatial fluctuations and the center of mass pinned at the center of the container represents only an extremely small volume measure of the available position space. These properties dominate the distribution and reduce the importance of the pinned COM position allowing for the found structure of the density and the averaged alignment of the H2orientation with the tip-tip direction. The shape of the container seems to be of minor influence as the boundary condition drives the wavefunction exponentially to zero within an interval small compared to the container radius. Thus, the energy expectation value is essentially determined via the probability distribution by the bulk of the position space between the origin and the container surface where it is forced by the strong long-ranged Coulomb field of the tips. Thus, the bulk probability distribution is influenced less than might be expected by the values at the origin and by the boundary. Preliminary checks on the orientation of H2 in a system with the same boundary condition but without tips yield rather homogeneous distributions over all angles and very similar average shapes on spherical domains where the COM is averaged. The ground state energy is depicted in Fig. 5 as a function of the tip-tip distance in comparison with the previous result for the fixed center of nuclear masses. Its scatter, as shown by the error bars, might be expected to dominate the scatter for fixed COM because of the increased number of degrees of freedom. Nevertheless, its size has appreciably decreased in the present optimization with a quantum mechanical COM variable and supports the validity of the present procedure. In contrast to the clamped COM, the new degrees of freedom drastically smooth the energy variation because the COM probability distribution can avoid strong density fluctuations by depleting regions with large potential differences. Single checks for D13 a.u. have shown the stability of the energy value against smoothing the wall’s hard cut-off by a repulsive potential as long as it remains below 0.1 a.u. in depth. A small elevation in energy is apparent that separates the region of larger tip-tip distance Dwith main density accumulated within the tip-tip disk from the small D region where the density is expelled to the outside of that disk. Comparing the amount of change of energy vs tip-tip distance, a reduction similar to the pinned COM 15 case is observed FIG. 3. Average probability density in sections 2(n1) a.u. ρ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y2þZ2 p ,2na.u. vs SXCOM coordinate for (a) D¼7:1 and (b) 13.1 a.u and for various values of n; compare with horizontal cuts in Fig. 2 at respective Dvalues. FIG. 4. Distribution of the angles (ΘX,ΘY,ΘZ), see sketch in (a) and definition in Fig. 1, between H2axis and (SX,SY,SZ)-coordinate axes from highest density areas (a) to lowest density (d) for D= 9.1 a.u. The keys denote the (SX,Sρ) areas plotted in Fig. 2 with Sρ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S2 YþS2 Z qequal to the distance from the symmetry axis. Note the different ordinate scales. ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-5 ©Author(s)2021 when degrees of freedom are freed. In the latter case, 15 the ground state energy of the potential-energy surface changes by 0.1 a.u. along a tip-tip distance change of 10 a.u., Fig. 2 of Ref. 15,as opposed to 0.005 a.u. energy change in the relaxed nonadiabatic case, Fig. 3 of Ref. 15. Here in Fig. 5, the further relaxation by the COM degrees of freedom yields a 0.0003 a.u. variation of energy. Nevertheless, one observes in Fig. 5 around 9 a.u. tip-tip distance a still rather abrupt change of slope which is, of course, more pronounced in the previous curve. This distance seems to represent the break below which the probability density is squeezed out. The dependence of the ground state energy on tip charge Q has been calculated with the parameter values optimized at Q¼0:5 a.u., i.e., without a full optimization sweep, in order to give a first impression about the charge’s influence. In Fig. 6, a linear behavior is displayed, as it would be seen by an electrostatic test charge in the tip’s field. For very small Q, the energy should deviate from the linear behavior to approach the free molecule value. IV. CONCLUSIONS The center of mass probability distribution of an H2molecule shows a strong anisotropy in decaying from its peak around the center to the limiting walls. A smooth decline perpendicular to the tip-tip axis along one single ridge contrasts with a rather steep variation parallel to the tip-tip axis for D13 a.u. tip-tip distance. With decreasing the tip-tip distance, the ridge splits into two ridges appearing as disks with axial symmetry near but clearly apart from the tip charges. Both increasingly separate and the volume occupied by the center of mass spreads over the volume between both tips. Below a tip-tip spacing of D9 a.u., those disks jump outside the region between the two tips and fill the outer volume, indicating that H2is squeezed out. The H2axis is directed on the quantum mechanical average all over the space parallel to the tip-tip axis. This is not explained by electrostatic arguments reasoning with compact H2dipoles anisotropically induced by the tips’field. The quantum mechanical probability spreading of course plays a role, though it does not explain the depletion of the perpendicular disks. Also, the orientation of the molecule axis deviates from the equatorial orientation found for a center of mass fixed at the origin. Simple arguments either based on the classical view or using the suppression of three degrees of freedom are not conclusive. We may view the H2protons as reacting to the external field more slowly than the electrons being expelled in the long term from the regions connected with the tip charge positions. The spatial H2distribution will be relevant when discussing the vibrational modes inelastically excited by tunneling currents. It might also be important for the container design to enclose a low density gas of that molecule, as the present probability reflects the density in the noninteracting case. H2can be held between two suitably separated positively charged tips on disks perpendicular to the tip-tip axis, while H2is expelled behind the tips when the tip-tip distance shrinks. Within the here considered range of external parameters, the orientation of the H2axis direction quantum mechanically fluctuates around its average parallel to the tip-tip axis and will influence via the transition matrix the strength of the coupling to the tips as well as the amount and the frequency of the vibrational excitations. We expect that the present results can be continuously scaled up. By increasing the container radius to an experimentally manageable size, this configuration presents a possibility of enclosing a few hydrogen molecules within a nanoscale region being subject to controlled external probing. ACKNOWLEDGMENTS This work has been supported in part by the Basque Departamento de Educación, Universidades e Investigación, the University of the Basque Country UPV/EHU (Grant No. IT1246-19) and the Spanish Ministerio de Ciencia e Innovación (Grant No. PID2019-107396GB-I00/AEI/10.13039/501100011033). The assistance of the DIPC computer center is thankfully acknowledged. W.S. is indebted to the Institute of Theoretical Physics and Astrophysics of CAU Kiel as well as to the university’s computer center. M.A.V.H. acknowledges financial support from the Collaborative Research Fund of the Research Grants Council of Hong Kong (Grant No. C2014-15G). ICTS is supported by the FIG. 5. Total energy plotted by black stars vs the distance between the tips as compared to the fixed COM result (circles). 15 The error bars denote the statistical uncertainty. FIG. 6. Energy vs tip charge for D¼5:1 and 13.1 a.u. tip-tip distance. ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-6 ©Author(s)2021 HKBU Institute of Creativity, which is sponsored by the Hung Hin Shiu Charitable Foundation. DATA AVAILABILITY Data sharing is not applicable to this article as no new data were created or analyzed in this study. REFERENCES 1 J.-M. Lehn, Supramolecular Chemistry: Concepts and Perspectives (Wiley-VCH, Weinheim, 1995). 2 R. A. Bissell, E. Córdova, A. E. Kaifer, and J. F. Stoddart, Nature 369, 133 (1994). 3 C. Cheng, P. R. McGonigal, S. T. Schneebeli, H. Li, N. A. Vermeulen, C. Ke, and J. F. Stoddart, Nat. Nanotechnol. 10, 547 (2015). 4 L. Grill, K.-H. Rieder, F. Moresco, G. Rapenne, S. Stojkovic, X. Bouju, and C. Joachim, Nat. Nanotechnol. 2, 95 (2007). 5 T. Kudernac, N. Ruangsupapichat, M. Parschau, B. Macia, N. Katsonis, S. R. Harutyunyan, K.-H. Ernst, and B. L. Feringa, Nature 479, 208 (2011). 6 G. Rapenne and C. Joachim, Nat. Rev. Mater. 2, 11 (2017). 7 R. H. M. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. van Hemert, and J. M. van Ruitenbeek, Nature 419, 906 (2002). 8 D. Djukic, K. S. Thygesen, C. Untiedt, R. H. M. Smit, K. W. Jacobsen, and J. M. Van Ruitenbeek, Phys. Rev. B 71, 161402(R) (2005). 9 A. Halbritter, P. Makk, S. Csonka, and G. Mihaly, Phys. Rev. B 77, 075402 (2008). 10 M. L. Trouwborst, E. H. Huisman, S. J. van der Molen, and B. J. van Wees, Phys. Rev. B 80, 081407 (2009). 11 H. Wang, S. Li, H. He, A. Yu, F. Toledo, Z. Han, W. Ho, and R. Wu, J. Phys. Chem. Lett. 6, 3453 (2015). 12 C. Lotze, M. Corso, K. J. Franke, F. von Oppen, and J. I. Pascual, Science 338, 779 (2012). 13 E. Carbonell-Sanromà, M. Corso, J. Li, A. X. Brión-Ríos, D. Sánchez-Portal, and J. I. Pascual, Surf. Sci. 678, 189 (2018). 14 M. Sobczyk, I. Anusiewicz, and P. Skurski, J. Chem. Phys. 118, 7297 (2003). 15 W. Schattke, T. Frederiksen, M. A. Van Hove, and R. Díez Muiño, Phys. Rev. B 101, 195436 (2020). 16 W. Kolos and L. Wolniewicz, J. Chem. Phys. 49, 404 (1968). 17 K. Pachucki and J. Komasa, Phys. Chem. Chem. Phys. 20, 26297 (2018). 18 G. D. Dickenson, M. L. Niu, E. J. Salumbides, J. Komasa, K. S. E. Eikema, K. Pachucki, and W. Ubachs, Phys. Rev. Lett. 110, 193601 (2013). 19 D. R. Yarkony, J. Phys. Chem. 100, 18612 (1996). 20 B. T. Sutcliffe and R. G. Woolley, Phys. Chem. Chem. Phys. 7, 3664 (2005). 21 N. M. Tubman, I. Kylänpää, S. Hammes-Schiffer, and D. M. Ceperley, Phys. Rev. A 90, 042507 (2014). 22 E. Ley-Koo, Rev. Mex. Fís. 64, 326 (2018). 23 T. Pang, Phys. Rev. A 49, 1709 (1994). 24 S. B. Doma, F. N. El-Gammal, and A. A. Amer, arXiv:1509.02567 (2015). 25 W. Kolos and L. Wolniewicz, J. Chem. Phys. 46, 1426 (1967). 26 R. L. Wilkins and S. Taylor, J. Chem. Phys. 48, 4934 (1968). 27 W. Schattke and R. Díez Muiño, Quantum Monte-Carlo Programming: For Atoms, Molecules, Clusters, and Solids, Monograph (Wiley-VCH, New York, 2013). 28 M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys. 121, 7110 (2004). ARTICLE avs.scitation.org/journal/jva J. Vac. Sci. Technol. A 39(5) Sep/Oct 2021; doi: 10.1116/6.0001138 39, 053206-7 ©Author(s)2021