The unbiased Diffusion Monte Carlo: a versatile tool for two-electron systems confined in different geometries

Computational codes based on the Diffusion Monte Carlo method can be used to determine the quantum state of two-electron systems confined by external potentials of various nature and geometry. In this work, we show how the application of this technique in its simplest form, which does not employ complex analytic guess functions, allows to obtain satisfactory results and, at the same time, to write programs that are readily adaptable from one type of confinement to another. This adaptability allows an easy exploration of the many possibilities in terms of both geometry and structure of the system. To illustrate this, we show results in the case of two-electron hydrogen-based species (H$_2$ and H$_3^+$) and two different types of confinement, nanotube-like and octahedral crystal-field.


Introduction
The fundamental relevance and concrete applications of the confined quantum systems has been established since the early days of quantum physics [Michels, 1937;Sommerfeld, 1938;De Groot, 1946], up to the present day [Sabin, 2009;Sen, 2009;Sen, 2014;Ley-Koo, 2018]. Simple confined physical systems have a very long tradition in physics because of their importance as basic theoretical problems [Aquino, 2007;Burrows, 2006;Micca Longo, 2015a;Micca Longo, 2015b;Micca Longo, 2015c;Longo, 2018], and as prototypes to gain insights in semiconductor physics (quantum dots [Chuu, 1992;Porras-Montenegro, 1992], quantum wires and quantum wells [Harrison, 2005;Munjal, 2018]). Atoms imprisoned in zeolite traps, clusters and fullerene cages provide good examples of confined systems in atomic physics and inorganic chemistry [Belosludov, 2003;Connerade, 1999;Hernandez-Rojas, 1996;Soullard, 2004]. Many of these studies on the physical-chemical properties of confined systems focus on systems with a small number of electrons: for example, species with small molecules containing hydrogen, helium and lithium. Furthermore, confinement is produced by means of an infinite potential well with different geometries, in order to produce analogies with different real systems without considering their details. In recent years, our group has developed a very simple and versatile methodology to study this kind of systems [Micca Longo, 2015a;Micca Longo, 2015b;Micca Longo, 2015c;Longo, 2018]. It is a modification of the Diffusion Monte Carlo (DMC) method: since the number of electrons is low, no variance reduction technique requesting a trial wavefunction is employed; instead, efficient methods of convergence of the eigenvalue are employed, inspired by adaptive control theory. Thanks to this method, it is possible to study different systems, with different orientations, and even quickly change the geometry of the potential well, from spherical to elliptical, from cylindrical to cubical, using Cartesian coordinates and implementing few modifications to the starting code. Our previous works [Micca Longo, 2015a;Micca Longo, 2015b;Micca Longo, 2015c;Longo, 2018] focused on the effects of various confinement geometries on the excited states of monoelectronic systems, such as H and H2 + . In this paper, we show the application of the method to the ground state of two-electron systems (H2 and H3 + ) with non-trivial confinement geometries. An exhaustive report on the many-electron confined system can be found in the work by Jaskólski [Jaskólski, 1996]: different methods of analysis and description of spatial confinement effects are reviewed, with a particular attention to their importance in semiconductor structures. Free H2 is one of the simplest molecular systems, and its electronic properties represent a key-point in understanding larger molecular systems; therefore, its confinement has received much attention during the last decades. The molecular hydrogen confined in a rigid spheroidal box, with fixed and not fixed nuclear positions, was studied with variational calculations and quantum Monte Carlo techniques [LeSar, 1981;LeSar, 1983;Pang, 1994;Colín-Rodríguez, 2010;de Oliveira Batel, 2018;Doma, 2018], also with Monte Carlo methods beyond the Born-Oppenheimer approximation [Sarsa, 2013]. Many studies concerning the ground state and the potential energy surfaces of the H3 + can be found in the literature [Conroy, 1964;Röhse, 1994;Jin, 2000;Pavanello, 2009;Adamowicz, 2012;Mizus, 2018; to cite a few], but almost no information is available about the quantum confinement of this interesting molecular system. Quantum calculations based on diffusion for the study of the free H3 + in equilateral triangle configuration trace back to the works by Anderson [Anderson, 1975;Anderson, 1992]. In this paper, we present the first results concerning the confined H2 inside a six-negative-charge cage and the collinear H3 + enclosed in a nanotube. All results are obtained by diffusion Monte Carlo method.

Method
Quantum Monte Carlo (QMC) is a class of computer algorithm that are able to simulate quantum systems and to compute the electronic ground state of atoms, molecules and solids. Among several Quantum Monte Carlo methods available [Foulkes, 2001], the diffusion Monte Carlo (DMC) method offers the possibility to consider confining walls of any geometry while retaining the use of cartesian coordinates and to include any distribution of positive charge within the cavity. Basically, DMC is a stochastic projector method that makes use of the similarity between the imaginary time Schrödinger equation and a generalized diffusion equation, which can be solved using stochastic calculus and simulating a random walk. The diffusion process of numerous "walkers" in the phase space is simulated. A walker is a mobile mathematical point in a set { # } of phase coordinates R in a 3N-dimensional space, where N is the number of physical electrons in the system, different from the number of walkers M, which is an integer numerical variable, ∼ 10 ) . An evolution is performed by repeated application of the short-time diffusionreaction propagator where t is a numerical step, and ? is the so-called energy offset that controls the walker total population. The kinetic energy term of the wave equation is connected with the diffusion process (the first exponential term of equation (1), including the normalizing constant), while the potential energy term (the second exponential term) controls the so called "birth-death algorithm", that is the destruction or the multiplication of walkers [Foulkes, 2001]. According to this second factor, a walker may randomly disappear, or duplicate itself, in a numerical time step . During the simulation, ? is constantly adjusted so that the average number of walkers is approximately constant [Thijssen, 2007]: where is a small positive parameter, and ? @ , ? @BC , # and #4I are the energy values and number of walkers at time step i and i-1.
At the numerical steady-state, ET becomes the estimate of the energy eigenvalue.
The method always provides the energy of the ground state of the system, with a very limited number of numerical parameters to be optimized (typically: energy offset ? , time step , number of walkers M ). We assume equal to to reduce the number of numerical parameters. A computational code has been developed by the authors, and it is carefully described in the previous works [Micca Longo, 2015a;Micca Longo, 2015b;Micca Longo, 2015c;Longo, 2018]. In order to reduce energy fluctuations and statistical deviations (typical of the DMC method), we implemented a sort of restart of the walker cloud, for each point of the potential energy curve. In a first simulation, the walker population is roughly initialized in a hypercube phase space: it diffuses according to equation (1), and the energy offset is adjusted according equation (2). At the end of this simulation, a new and more realistic configuration of the walker cloud is obtained. This new configuration represents a sort of "guiding wavefunction", and it is used as a starting point for a second simulation, which then gives a more precise energy eigenvalue.
In contrast with recent Monte Carlo calculations [Sarsa, 2012], our programs do not employ an analytic expression of a guess wavefunctions to reduce variance. Instead, it performs long imaginary time evolutions, in the course of which averages are collected. This is the reason we name it "unbiased": it is essentially the most primitive form of the algorithm [Anderson, 1975] with some corrections. A first solution consists of a simple technique that tries to reduce the oscillations of the energy selfvalue during the adaptation process algorithm: occasionally, this eigenvalue is replaced with an average calculated in the previous phases of the calculation at, this is done at a moment of the "oscillation" of ET, when it crosses a momentary mobile average.
The second solution, as mentioned before, consists in using the walker cloud distribution obtained in a previous calculation and the related eigenvalue estimate as starting points for a calculation with a slightly different value of the nuclear positions: in this way, the previous calculation works as an estimate of the assumption to be obtained in the new calculation. While this technique leads to somewhat larger oscillations and fluctuations compared to the employ of analytical guesses, it has the advantage to be applicable to different confinement geometries with minimal variations, without having to change the expression of the guess wavefunction, since no such guess is used. DMC can handle only positive distributions, so many problems arise when we want to deal with a many-electron wave function. The DMC method can address low excited states (whose wavefunctions are not always positive) with distinct symmetries as well, by adding nodal surfaces that reproduces the excited state symmetry. However, in the case of two electrons systems, an unrestricted, six-dimensional wavefunction ( , ) for the ground state can be calculated directly. This means that the electron correlation is accounted exactly, with a remarkably low analytical fatigue: no expansion of the wavefunction is needed, no function database need be included, at least, in the unbiased approach. In this work, we are dealing with such two-electron wavefunctions. Our native codes are therefore updated to simulate two-electron, confined systems thanks to the aforementioned features. All equations that govern the diffusion and branching algorithms can be changed accordingly with to necessary to abandon cartesian coordinates. As a consequence, we can apply the DMC method to the confined H2 and H3 + with great freedom concerning confinement geometry, by simulating the six-dimensional system directly with the additional constraint that no six-dimensional walker can cross the confining surfaces.

Cylindrical confined H3 +
The confinement by cylindrically symmetrical potential surfaces represents a good model for the chemical confinement of small molecules within nanotubes, such as carbon nanotubes. A quantum computation allows us to determine the potential energy surface of these species in the confined state. The potential energy surface can then be used to make structural studies, in particular to determine the equilibrium position of the nuclei, and dynamic studies, that lead to predictions of infrared and Raman spectroscopy. In light of the increasing importance of studies related to the confinement of hydrogen-based species, in various kind of materials specially designed to store large quantities of hydrogen (i.e., carbon nanotubes) [Coppola, 2019], it is important to undertake some first confinement studies of hydrogen-based species in nanotubes. In this section of the work, we show results in the case of H3 + : this ion is the simplest triatomic species with only two electrons, so it is of a considerable interest from a more fundamental point of view. Since these are the first calculations for this confined system, we have considerably reduced the complexity of the problem: in particular, we are assuming that the molecule-ion is collinear, with the molecular axis placed on the nanotube axis. We also assume that the two internuclear distances, from the central nucleus to the two terminal nuclei, are equal to each other and are indicated by d: the full system keeps therefore a D∞h symmetry (Figure 1). This limitation is assumed only to reduce the number of independent parameters: the code can treat any other symmetry without modification. The cylinder radius, c, is the confinement characteristic length. The results in figures 2 and 3 show what can be obtained using the DMC method in its primitive form, without the use of test functions. In figure 2, the comparison with the best available results, in the case of the unconfined system, shows how a very simple DMC code algorithm is able to obtain a satisfactory prediction with small deviations from the reference points in most cases. By inspecting the figure, it is also clear that the trend of the calculated solution is satisfactory: therefore, in a future application to ab-initio molecular dynamics, the method may be exploited to produce a reliable interpolated potential energy surface. Statistical methods can also be used to further improve the results. At this point, thanks to the versatility of the algorithm, it is possible to introduce readily the cylindrical infinite potential barrier illustrated in Figure 1. Figure 3 shows the results of the compression of the electron cloud of the ion molecule by a cylindrical barrier with decreasing radius: the reduction of the available volume leads to a general lifting of the potential curve. The most important result is the shift towards smaller values of the equilibrium distance: this means that the molecule contracts due to the increase in the density of the electron cloud. Furthermore, a greater curvature of the potential surface is obtained at the minimum, which corresponds to an increase in the frequency of the Raman line of symmetrical stretching, that in principle could be experimentally measurable.

Crystal-field confined H2
In the previous section, we have shown the effect of a confinement by an infinite potential surface and with a geometric section: this is essentially a mathematical abstraction of more complex chemical problems. Actually, this approach (potential well confinement) was followed by the vast majority of the literature on species. With our algorithm, however, it is also possible to study a somewhat more realistic confinement mechanism, closer to the chemical reality of the systems. The mechanism, inspired by the premises of Bethe's crystalline field theory [Basolo, 1986], is based on the idea of using point electric charges in realistic geometric configurations of a crystal lattice, and which therefore produce a static potential surface. This acts on the charged species inside the molecule, in particular on the electrons, producing also in this case a confinement effect and an increase in energy. Furthermore, we expect a modulation of the potential energy surface, when the species included in the crystalline field are molecular species. Figure 4 shows a sketch of the H2 molecule confined in a six-negative-charge cage: the hydrogen molecule has a molecular axis directed towards two of the opposite charges of the octahedron, and its center of mass coincides with the center of mass of the octahedron. We therefore have a D4h symmetry.  Figure 4. H2 confined by a system of six negative point charges, with the hydrogen molecule-point charges system in a configuration of D4h symmetry (although the octahedron is assumed to remain geometrically perfect, the symmetry is broken by the molecule). In Figure 5, the energies of a confinement as from Figure 4 are reported for several distances l of the point-charges from the center of the octahedron. Each point is the result of the average of ten runs; increasing l brings to a larger standard deviation, which means that the methods works even better for the strongly confined systems which are the target of our studies. Figure 5 shows some of the results that can be obtained with our code: as the point charges approach, even in this case an effect of general energy increase is observed, due to the Coulomb repulsion, and the trend of the potential energy surface is significantly modified. This is a promising beginning for a geometry that deserves further investigation, also by considering other possible symmetries, for example D3d.

Conclusions
In conclusion, one may wonder why reconsider, in the context of a new emerging class of problems, a Monte Carlo calculation technique in its primitive form, that dates back to the 1970s. The reason lies in the contrast between analytical fatigue and computational cost. The variance reduction techniques used for Monte Carlo methods require the formulation of fully correlated tentative wave functions: this formulation is not immediate in the case of confined systems, and it is not immediately transportable from one confinement geometry to another. On the other hand, the accuracy of the unbiased Monte Carlo method can be increased by accurate statistical analysis and simple solutions. Nowadays, the large availability of fast computing resources, even portable ones, might lead a researcher to favor an approach that allows to develop native, short and readable programs, and, at the same time, allows to explore different confinement geometries with minimal modifications. Even with the limits due to statistical fluctuations, such algorithm takes the correlations between electrons into full account. Furthermore, our technique, despite being in the spirit of the primitive one, includes improvements of straightforward implementation, such as the averaging of the energy oscillations and the method of restarting the walker cloud. Ultimately, our experience suggests that a simpler technique, employing no analytic functions of significant complexity, can be a resource in the exploration of the wide range of possibilities that has been opened up by the recent reconsideration of confined quantum systems.