Confined hydrogen atom: endohedrals H@C₃₆ and H@C₆₀

In the majority of problems the time-independent Schrödinger equation, the central object in quantum mechanics does not admit exact solutions in terms of elementary or special functions. In such cases, it is necessary to solve this equation by means of numerical and approximate methods which inherently carry a certain degree of accuracy. For instance, some of the most common numerical methods are the following: the Numerov method, the spline-based method, the finite difference, the finite element, the Lagrange mesh method, and the variational one, to name a few. Interestingly, artificial intelligence algorithms have emerged as a promising tool for tackling eigenvalue problems of differential and integrodifferential operators [1, 2] with important applications on the solution of Schrödinger equation [3–7] and open quantum systems [8–10]. In particular, artificial neural networks have shown to be accurate and efficient algorithms for calculating eigenvalues and eigenfunctions, due mainly to their high capabilities to approximate complex functions from data [11–13].

Classical numerical schemes are widely used in the study of both spatially confined and unconfined quantum systems. In particular, the investigation of spatially confined quantum systems has gained much interest because some of their physical properties change abruptly with the size of the confining barrier. Furthermore, many physical phenomena can be modeled by a confined quantum system for example: atoms and molecules subject to high external pressures, atoms and molecules in fullerenes, inside cavities such as zeolite molecular sieves or in solvent environments, the specific heat of a crystalline solid under high pressure, etc. A complete list of applications can be found in the articles and reviews on the subject [14–26].

Eighty-five years ago, Michels et al [27] calculated the variation of the polarizability of hydrogen under high pressure. To this end, they proposed a simple model of the confined hydrogen atom where, to a first approximation, the nucleus is anchored in the center of a spherical box of radius r₀ and impenetrable walls. It is assumed that this infinite potential is due to the presence of neighboring negative electric charges. The corresponding wave function vanishes at the surface of the sphere, i.e. it must obey Dirichlet boundary conditions. Since then, this model has been successfully applied in the study of confined many-electron atoms and molecules (see [28, 29] and references therein). However, it takes into account the effects of repulsive forces only. A more realistic potential that embodies an attractive force considers a softer confinement in cavities of penetrable walls. The simplest penetrable confining potential of this type is that of a step potential [30, 31].

Another relevant example of a penetrable potential is realized by the logistic potential used in the analysis of the entropic properties of hydrogen atom [32].

Likewise, in the study of atoms inside fullerenes, the corresponding potential has been modeled effectively by the attractive spherical shell [18], the δ-potential [33] and a Gaussian spherical shell [34, 35]. Moreover, a spatially confined system inside an impenetrable spherical barrier, in addition to the inverted Gaussian potential [36–39], was employed by Lumb et al to describe the emission of two photons in the hydrogen atom [40].

It has been recognized as an advantage to study free systems, where bound states exist and the potential vanishes $V\rightarrow 0$ at large distances, using the methods applied to the spatially confined systems (see [41–43]). In such a procedure the original system with potential V is additionally enclosed inside an impenetrable box of radius r₀. Afterward, the wave functions as well as the eigenvalues are calculated in a sequence by increasing the value of this parameter r₀ until the desired accuracy is achieved. Mathematically, it has been demonstrated that this procedure converges to the exact results of the original system where the spatial confinement is absent (see [41–43]). Hence, the artificial spatial confinement gives highly accurate results with respect to the exact values, which can not be guaranteed by other methods [41]. Some examples to illustrate this point are the free hydrogen atom, an electron in the Yukawa potential, the Hulthen potential, and the singular oscillator (see [41–43]).

In this work, in the Born-Oppenheimer approximation, we calculate the energies and wave functions of the lowest states of a hydrogen atom confined inside a spherical penetrable barrier. This system has been used in [34, 35] (and references therein) to model the hydrogen atom confined by C₃₆ and C₆₀ fullerenes, where the corresponding confining potential is described by an inverted Gaussian function. Here, a more accurate and systematic study is carried out. To solve the Schrödinger equation we employ three different numerical methods: the Lagrange mesh, finite difference, and finite element. They are complemented by the training and use of an artificial neural network. For the sake of self-consistency, a brief description of the aforementioned numerical methods is presented in the next section with the corresponding references.

The goal of the present study is threefold. Firstly, we aim to improve previous results reported in the literature. Secondly, we make a comparison of the numerical results obtained with different methods and discuss the advantages of each individual approach. Finally, we show that the complementary use of an artificial neural network leads to efficient numerical interpolation.

The Schrödinger equation of hydrogen atom confined by a C₃₆ and C₆₀ endohedral cavities, in atomic units $(\hbar = m_e = e = 1)$, is of the form

$$\begin{equation} \bigg[\,-\frac{1}{2} \Delta \ - \ \frac{1}{r} \ + \ w(r)\,\bigg]\,\psi\ = \ E\, \psi , \end{equation} \tag{ 1 }$$

where Δ is the three-dimensional Laplacian, and the Gaussian spherical barrier w(r) is given by

$$\begin{equation} w(r)\ = \ -\,\omega_0\, \textrm{exp}[-(r-r_c)^2/\sigma^2], \end{equation} \tag{ 2 }$$

here, ω₀ is the well depth, r_c is the position of the center of the peak and σ is the width of the Gaussian, respectively. In concrete calculations, for the purposes of comparison with the existing results in the literature, from time to time the parameters r_c and σ will be presented in Angstroms (Å).

As is usual for any central potential, the angular momentum is conserved and the solutions of the Schrödinger equation (1) can be factorized in spherical coordinates ($r,\theta,\phi$). Explicitly, they read

$$\begin{equation} \psi_{nlm}(r,\theta,\phi) \ = \ R_{nl}(r) \,Y_{l,m}(\theta,\phi) , \end{equation} \tag{ 3 }$$

with $Y_{l,m}(\theta,\phi)$ being a spherical harmonic function. In the standard approach, it is convenient to introduce an auxiliary function $u_{nl}(r)$ defined by the relation $R_{nl}(r)\ = u_{nl}(r)/ r$. In this case, from (1) and (3) it follows that $u_{nl}(r)$ must obey the isospectral radial problem

$$\begin{equation} \bigg[\,-\frac{1}{2}\frac{d^2}{d r^2} \ + \ V_\textrm{eff}(r) \,\bigg] u_{nl}(r)\ = \ E\, u_{nl}(r), \end{equation} \tag{ 4 }$$

with the boundary condition $u_{nl}(r = 0) = 0$. The spectral problem (4), is defined in the domain (configuration space) $r\in [0,\,\infty)$. The effective potential appearing in (4) is given by

$$\begin{equation} V_\textrm{eff}(r)\ =\ -\frac{1}{r}\ + \ \frac{l(l+1)}{2\,r^{2}}\ + \ w(r), \end{equation} \tag{ 5 }$$

$l = 0,1,2,\ldots,(n-1)$ is the quantum number of angular momentum. Formally, (4) describes a one-dimensional particle of unit mass in the half positive line with an effective potential $V_\textrm{eff}(r)$. The shape of $V_\textrm{eff}$ is displayed in figure 1 for fixed l = 0, $\omega_0 = 1$ a.u., σ = 0.5 a.u. and four different values of r_c (the position of the center of the peak); there exists a critical value $r_c^\textrm{crit}\approx 1.05$ a.u. such that for $r_c\gt r_c^\textrm{crit}$ the effective potential develops a second minima. In turn, figure 2 shows the behavior of $V_\textrm{eff}$ (5) as a function of the width σ of the Gaussian for fixed ω = 0.5 a.u. and $r_c = 3.54$ Å.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to solve the radial Schrödinger equation (4), three different methods will be employed: (i) The Lagrange-mesh method (LMM), (ii) finite difference method (FDM), and (iii) finite element method (FEM). By combining our accurate results with an artificial neural network, we construct an efficient numerical interpolation for the energies.

2.1. The Lagrange-Mesh method

In the context of the LMM [44–46], a set of N Lagrange functions $f_i\,(x)$ defined over the domain of the radial variable is associated with N mesh points x_i which correspond to the zeros of Laguerre polynomials of degree N, i.e. $L_N(x_i) = 0$. The Lagrange-Laguerre functions $f_i\,(x)$ which satisfy the Lagrange conditions

$$\begin{equation} f_i(x_j)=\lambda_i^{-1/2}\delta_{ij}, \end{equation} \tag{ 6 }$$

at the N mesh points are given by

$$\begin{equation} f_i(x)=(-1)^{i} \frac{x}{x_i^{1/2}}\frac{L_N(x)}{(x-x_i)}e^{-x/2}. \end{equation} \tag{ 7 }$$

The coefficients λ_i are the weights associated with a Gauss quadrature

$$\begin{equation} \int_0^{\infty} G(x) dx \approx \sum_{k=1}^{N} \lambda_k\, G(x_k). \end{equation} \tag{ 8 }$$

In terms of the N Lagrange functions $f_i(x)$ (7), the solution of the Schrödinger equation (4) is expressed as

$$\begin{equation} \psi(r) =\sum_{i=1}^{N}c_i f_i(r). \end{equation} \tag{ 9 }$$

The trial function (9), together with the Gauss quadrature (8) and the Lagrange conditions (6) leads to the system of variational equations

$$\begin{equation} \sum_{j=1}^N \left[h^{-2}\,T_{ij}+\left(\frac{l(l+1)}{2\, h^2\,x_i^2}-\frac{1}{h\,x_i}+\omega(h\,x_i)\right)\delta_{ij}\right]\,c_j = E\,c_i, \end{equation} \tag{ 10 }$$

where $\omega(x_i)$ is the Gaussian potential (2) evaluated at the mesh points x_i and T_ij are the kinetic-energy matrix elements whose explicit expression is found in [46]. h is a scaling factor that allows to adjust the mesh to the system in consideration. By solving the system (10), not only the energies E are obtained but also the eigenvectors c_i from which the approximation to the wave function (9) is obtained.

Inside the LMM approach, the expectation value of the radial coordinate r is easily calculated. Given the approximation to the wave function (9) together with Gauss quadrature (8) and the Lagrange condition (6) leads to

$$\begin{equation} \langle r \rangle = \sum_{i=1}^{N} c_i^2\, x_i, \end{equation} \tag{ 11 }$$

where c_i are the eigenvectors resulting from solving (9) and x_i are the mesh points.

2.2. Finite difference method

The FDM is a numerical approach easy to implement on a computer, and it is used to solve ordinary as well as partial differential equations in an approximate way. The method is based on the discretization of the Hamiltonian on a spatial grid. Consider the system defined in a domain $[0,L]$, the Schrödinger equation is then solved on a uniform grid [47] defined by the set of discrete points which are the nodal points $\left\{r_{j}\right\}$. The domain is split into N subintervals of equal length h:

$$\begin{equation} 0=r_{1}\lt r_{2}\lt \ldots \lt r_{N+1}=L. \end{equation} \tag{ 12 }$$

where

$$\begin{equation} h=r_{j+1}-r_{j}. \end{equation} \tag{ 13 }$$

In the present work, we use a fourth-order centered difference approximation to the second derivative, namely

$$\begin{align} \begin{array}{l} {{\tilde u}^\prime }^\prime = \;\frac{{ - \tilde u({r_{j + 2}}) + 16\tilde u({r_{j + 1}}) - 30\tilde u({r_j})}}{{12{h^2}}}\quad \\ \,\,\,\,\,\,\,\,\,\, + \frac{{16\tilde u({r_{j - 1}}) - \tilde u({r_{j - 2}})}}{{12{h^2}}} + \;O({h^4}). \end{array} \end{align} \tag{ 14 }$$

Substituting this definition on the Schrödinger equation (4) we find that there is one equation for each j value. For j = 1 and j = N we impose the boundary condition $\psi (r_{1}) = \psi (r_{N+1} = L) = 0$. Hence, these systems of equations can be written as an eigenvalue problem $H\,C\, = \,E\,C$ where the matrix H is a pentagonal matrix, whereas the vector C contains the values of $\tilde{u}$ evaluated on the grid, i.e. $C_{i} = \tilde{u}(r_{i}).$ To solve the differential equation, we select the parameters of the potential and angular momentum l of the states, we propose a value of L much higher than r_c , and solve the eigenvalue problem. We increase the value of L and solve the eigenvalue problem again, compare the energy values and repeat this process until the desired accuracy is achieved. We call $L_{\max }(l)$ to the value of L at which the desired accuracy is achieved. Note that the boundary condition at $\psi (r_{N+1} = L) = 0$ is equivalent to the existence of an impenetrable barrier at that precise point, this fact establishes the relationship with the confined systems. In this case, the Schrödinger equation to be solved can be written as (cf (4)):

$$\begin{equation} -\frac{1}{2}\frac{d^{2}\tilde{u}_{nl}}{dr^{2}}\,+\,V_{\mathrm{eff}}(r)\, \tilde{u}_{nl}\,+\,V_{c}\,\tilde{u}_{nl}\,=\,E_{\max }\,\tilde{u}_{nl}, \end{equation} \tag{ 15 }$$

where $V_{\mathrm{eff}}(r)$ is the effective potential (5) and the confinement potential V_c is defined as follows

$$\begin{equation} V_{c}(r)=\left\{ \begin{array}{l} 0, \qquad \thinspace\ if \ r\leqslant L\thinspace , \\ \infty , \qquad if \ r\gt L\thinspace . \end{array} \right. \end{equation} \tag{ 16 }$$

The above spectral problem (15) is defined in the reduced domain $r\in \lbrack 0,\,L]$. The Schrödinger equation (15) was studied by Lumb et al [40] in connection with the two-photon emission in confined hydrogen atom. In the region r < L the Schrödinger equation becomes

$$\begin{equation} -\frac{1}{2}\frac{d^{2}\tilde{u}_{nl}}{dr^{2}}\ +\ V_{\mathrm{eff}}(r)\ \tilde{u}_{nl}\ =\ E_{\max }\,\tilde{u}_{nl}, \end{equation} \tag{ 17 }$$

cf (4). The function $\tilde{u}$ satisfies the Dirichlet boundary conditions:

$$\begin{align} \tilde{u}_{nl}(0)\ & =\ 0,\nonumber\\ \tilde{u}_{nl}(L)\ & =\ 0. \end{align} \tag{ 18 }$$

The original energy eigenvalues and eigenfunctions of the free (spatially unbounded) system (4) are recovered in the limit $L\rightarrow \infty$,

$$\begin{equation} \tilde{u}_{nl}(r)\,\rightarrow \,u_{nl}(r),\qquad E_{\max }\,\rightarrow \,E. \end{equation} \tag{ 19 }$$

2.3. Finite element method

The FEM combines the discretization of the space in elements (intervals) and the use of polynomial interpolating functions on each element. This method is widely used in engineering, classical physics, and quantum mechanics problems, among others. The discretization is based on the reformulation of the relevant differential equation as an equivalent variational problem. The Garlekin methods are then employed in the corresponding minimization procedure [48]. In general, one can identify the following steps to solve the associated differential equation: (i) to present the problem in a variational formulation, (ii) a discretization of the domain using FEM, and finally, (iii) to find the solution of the discrete problem, which may consist of the solution of a system of simultaneous equations or an eigenvalue problem.

In Quantum Mechanics the FEM was used from a few years ago [47, 49–54]. An excellent introduction to FEM in Quantum Mechanics can be found in the Ram-Mohan's book [52]. Only the key points of the method will be presented here.

The time independent Schrödinger equation for a particle of mass m subjected to a potential energy $V({\mathbf r})$ is given by:

$$\begin{equation} \bigg[\,-\frac{\hbar ^{2}}{2\,m}\Delta\ +\ V({\mathbf r})\,\bigg]\,\psi ({\mathbf r})\ =\ E\,\psi ({\mathbf r}) . \end{equation} \tag{ 20 }$$

This equation can be obtained by minimizing the following action integral I:

$$\begin{equation} I[\psi,\,\psi ^{\ast }]=\int d^{3}r\left[ \frac{\hbar ^{2}}{2m}\nabla \psi ^{\ast }\cdot \nabla \psi +\psi ^{\ast }\left( \,V\,-\,E\,\right) \psi \right] , \end{equation} \tag{ 21 }$$

where $\psi ^{\ast }$ and ψ are considered as two independent 'fields'. It is assumed that ψ is continuous up to its second derivative. By varying the action I with respect to $\psi ^{\ast }$ we obtain the Schrödinger equation (20). For a spherical symmetric potential $V = V(r)$, only the dynamics of the r coordinate is not trivial. The problem reduces to the one-dimensional case in the r-space, the variable r varies in the interval $\left[ 0,L\right]$. Hence, the action integral I (in atomic units) becomes:

$$\begin{equation} I=\int_{0}^{L}dr\left\{\frac{1}{2}\frac{d\psi ^{\ast }}{dr}\frac{ d\psi }{dr}+\psi ^{\ast }\left[ V_\textrm{eff}(r)-E\right] \psi \right\} , \end{equation} \tag{ 22 }$$

where $\psi = \psi(r)$, and $V_\textrm{eff}(r)$ is an effective potential $V_\textrm{eff}(r) = V(r) + \frac{l\,(l+1)}{2\,r^2}$, cf (5).

Now, the interval $\left[ 0,L\right]$ is divided into small subintervals called elements. The action integral (22) can be decomposed as the sum of the action computed in each element,

$$\begin{equation} I\ = \ \sum_{j=1}^{n}I^{(j)}, \end{equation} \tag{ 23 }$$

where n is the number of elements and $I^{(j)}$ is the action integral evaluated on the jth element. Explicitly, the wave function ψ_j defined in the jth element is expanded as a linear combination

$$\begin{equation} \psi _{j}(r)\ = \ \sum_{j=1}^{n} c_{j}\,N_{j}(r), \end{equation} \tag{ 24 }$$

where c_j , $j = 1,2,\ldots,n$, are unknown coefficients to be determined whereas $N_{j}(r)$ are interpolating polynomials. These are defined for the jth element and they are identically zero out of this element.

The basic idea is to make the variation of the action integral with respect to the coefficients $c_{i}^{\ast }$,

$$\begin{equation} \frac{\delta I}{\delta c_{i}^{\ast }}\ =\ 0, \quad i=1,2,\ldots,n , \end{equation} \tag{ 25 }$$

and by solving these equations to obtain the optimal energies and eigenfunctions.

In particular, Guimaraes and Prudente [53] developed an alternative version of FEM called the p-finite element method (pFEM) to study the confined hydrogen atom, and Nascimento et al [34] employed a pFEM version to study the electron structure of endohedrally confined atoms using an attractive Gaussian potential to model atoms inside fullerenes.

As L becomes very large the energies of a confined system approach to those of the confinement-free system. In practice, a large value of both L and the parameter n, and a polynomial degree for the $N_{j}(r)$ appearing in (24) are chosen and then, the generalized eigenvalue problem is solved. For a fixed L, by increasing either the value of n or the polynomial degree, or both, higher precision in the results can be achieved as we will explain in the next section.

2.4. Artificial neural networks

Artificial intelligence has emerged as a collection of computational techniques which seek to mimic the human brain in order to complete tasks for which standard algorithms lead to partially satisfactory results or are costly to implement [55–58]. Particularly, neural networks are artificial intelligence algorithms inspired by the workings of neurons in the human brain. These algorithms have demonstrated the capacity of pinpointing relevant specific pieces of information 'buried' in huge data sets and unveiling complex non-linear relationships between the inputs and target, which would be all but impossible to accomplish through a standard visual inspection [59–61].

In this work, we implement a neural network to estimate the eigenvalues of the radial Schrödinger equation (4) for different positions of the peak of the Gaussian, r_c . This neural network consists of a two-layer feed-forward architecture with ten neurons in the hidden layer and a linear neuron in the output layer. Note that the input node of our neural network corresponds to the center of the inverted Gaussian potential, r_c . Figure 3 shows the architecture of the regression neural network. In general, the output of each neuron before the activation function reads,

$$\begin{equation} z=\sum_{i=1}^N \omega_i x_i, \end{equation} \tag{ 26 }$$

where ω_i are the synaptic weights, x_i are the inputs, and N is the number of inputs. Importantly, all the neurons in the output layer contain linear activation functions whereas the neurons in the hidden layer have sigmoid functions given by

$$\begin{equation} \alpha(z)=\frac{1}{1+e^{-z}}. \end{equation} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

Synaptic weights of the neural network are optimized with the Levenberg-Marquardt backpropagation method in a direction that minimizes the mean squared error [62, 63]. This method approaches second-order training speed without having to compute the Hessian matrix. Because the performance function is given by a sum of squares then the Hessian matrix can be approximated as $\mathbf{H} = \mathbf{J}^T\mathbf{J}$, where J is the Jacobian matrix. Using this approximation, the synaptic weights can be updated by the following expression,

$$\begin{equation} \mathbf{w}_{k+1}=\mathbf{w}_{k}-[\mathbf{H}+\mu I]^{-1}\mathbf{J}^{T}\mathbf{e} \end{equation} \tag{ 28 }$$

where k denotes the kth iteration, µ is the learning rate and e is the vector of networks errors. To train the neural network, we use a subset of the dataset that contains eigenvalues of the radial Schrödinger equation (4) for different positions r_c of the peak of the Gaussian. These eigenvalues are calculated by the LMM. After the training stage, the neural network is able to predict the eigenvalues of the whole dataset with a coincidence at six significant digits at least.

In order to describe the confinement of the hydrogen atom inside C₃₆ and C₆₀ fullerene cages, an attractive spherical Gaussian potential [34] is considered (1). Energies, eigenfunctions and expectation values of r for the first six states of $s, p$ and d symmetries are accurately calculated. In order to be able to compare with previous results, the values of the parameters of the Gaussian potential (1) that are investigated in detail are: $r_c = 2.5$ Å (for H@C₃₆) and $r_c = 3.54$ Å (for H@C₆₀), $\omega_0 = 0.5$ a.u. and $\sigma = 0.26, 0.57$ and 1.59 Å. Previous work has shown [64–66] that the Lagrange-mesh (LMM) method provides very accurate values for the energy values. The corresponding results obtained using the LMM are shown in tables 1–3 with 12 significant figures; overall within this precision, as a result of explicit computations it turns out that the LMM is the most practical, accurate and easier to implement an approach for the present case (see discussion below). We emphasize that, in principle, the three distinct numerical approaches (i) LMM, (ii) fourth order finite differences, and (iii) the FEM, can lead to extremely highly accurate results by choosing more refined values of the parameters of the methods with the corresponding growth of the computation time.

Table 1. Energy E and expectation value $\langle r \rangle$ of the hydrogen atom of the $1s \ldots 6s$ states in presence of a inverted Gaussian potential (2) centered at $r_c=2.50$ and 3.54 Å with $w_0=0.5$ a.u. as a function of σ. The free case corresponds to σ = 0. For $r_c=3.54$ Å comparison is done with results presented in [35]. The values for E and $\langle r \rangle$ are displayed with 12 significant figures, and they correspond to the results obtained with the Lagrange Mesh Method. The coincidence with the FEM results is in 7–8 digits whilst for the FDM, the agreement with the results of the LMM is in all figures.

		$r_c=2.50$ Å		$r_c=3.54$ Å
ns	σ Å	E a.u.	$\langle r \rangle$ a.u.	E a.u.	[35]	$\langle r \rangle$ a.u.
1s	0.00	−0.500 000 000 000	1.500 000 0000	−0.500 000 000 000		1.500 000 000 00
	0.26	−0.505 803 094 144	1.625 016 8416	−0.500 226 076 582	−0.500 2261	1.510 216 498 28
	0.57	−0.528 322 517 980	2.102 338 6012	−0.501 274 477 556	−0.501 2745	1.562 697 296 13
	1.59	−0.700 338 868 882	2.292 977 5088	−0.558 460 325 443	−0.558 4603	3.053 971 9446
2s	0.00	−0.125 000 000 000	6.000 000 0000	−0.125 000 000 000		6.000 000 000 00
	0.26	−0.240 727 300 618	4.641 295 9696	−0.222 678 640 702	−0.222 6786	6.334 049 117 02
	0.57	−0.352 513 053 170	4.066 550 9890	−0.341 831 613 201	−0.341 8316	6.421 677 897 48
	1.59	−0.493 011 698 186	4.218 131 7661	−0.489 180 987 884	−0.489 1810	5.057 577 1424
3s	0.00	−0.055 555 555 556	13.500 000 000	−0.055 555 555 556		13.500 000 0000
	0.26	−0.064 455 619 315	12.272 530 044	−0.056 490 840 224	−0.056 490 84	13.650 810 7402
	0.57	−0.069 257 879 934	11.329 276 155	−0.063 868 006 446	−0.063 868 01	11.172 615 4458
	1.59	−0.204 794 618 028	5.969 868 8054	−0.247 983 103 635	−0.247 9831	6.542 671 5036
4s	0.00	−0.031 250 000 000	24.000 000 000	−0.031 250 000 000		24.000 000 0000
	0.26	−0.034 143 357 000	22.414 917 645	−0.031 553 410 884	−0.031 553 41	23.564 684 5778
	0.57	−0.036 157 620 392	21.040 808 457	−0.036 248 303 973	−0.036 248 30	19.451 723 6193
	1.59	−0.056 183 917 194	13.817 270 956	−0.070 803 944 469	−0.070 803 95	10.888 054 490
5s	0.00	−0.020 000 000 000	37.500 000 000	−0.020 000 000 000		37.500 000 0000
	0.26	−0.021 319 794 124	35.527 535 710	−0.020 260 703 774		36.568 599 620
	0.57	−0.022 355 880 943	33.759 026 704	−0.022 999 812 535		31.563 158 8038
	1.59	−0.030 981 035 591	24.577 416 804	−0.034 734 278 117		22.431 757 286
6s	0.00	−0.013 888 888 889	54.000 000 000	−0.013 888 888 889		54.000 000 0000
	0.26	−0.014 606 403 251	51.633 586 493	−0.014 097 009 471		52.672 305 540
	0.57	−0.015 206 956 729	49.482 686 203	−0.015 733 494 004		46.870 636 8157
	1.59	−0.019 751 915 407	38.280 917 640	−0.021 384 692 368		35.795 142 275

Table 2. Energy E and expectation value $\langle r \rangle$ of the hydrogen atom of the $2p, \ldots, 7p$ states in presence of a inverted Gaussian potential (2) centered at $r_c=2.50$ and 3.54 Å with $w_0=0.5$ a.u. as a function of σ. The free case corresponds to σ = 0. For $r_c=3.54$ Å comparison is done with results presented in [35]. The values for E and $\langle r \rangle$ are displayed with 12 significant figures, and they correspond to the results obtained with the Lagrange Mesh Method. The coincidence with the FEM results is in $7-8$ digits whilst for the FDM, the agreement with the results of the LMM is in all figures.

		$r_c=2.50$ Å		$r_c=3.54$ Å
np	σ Å	E a.u.	$\langle r \rangle$ a.u.	E a.u.	[35]	$\langle r \rangle$ a.u.
2p	0.00	−0.125 000 000 000	5.000 000 0000	−0.125 000 000 000		5.000 000 0000
	0.26	−0.235 017 656 232	4.607 321 7318	−0.205 773 905 331	−0.205 7739	6.194 941 0697
	0.57	−0.358 937 721 074	4.600 980 5888	−0.321 662 542 152	−0.321 6625	6.497 295 5252
	1.59	−0.524 721 262 609	4.498 847 2199	−0.487 790 459 579	−0.487 7905	6.400 068 2529
3p	0.00	−0.055 555 555 555	12.500 000 000	−0.055 555 555 555		12.500 000 000
	0.26	−0.059 201 294 702	12.295 298 412	−0.058 921 847 548	−0.058 921 85	9.714 367 4959
	0.57	−0.064 787 981 905	10.822 752 863	−0.072 949 413 182	−0.072 949 41	6.510 334 8451
	1.59	−0.248 825 019 321	5.325 475 7370	−0.253 036 633 049	−0.253 036 63	6.150 340 6694
4p	0.00	−0.031 250 000 000	23.000 000 000	−0.031 250 000 000		23.000 000 000
	0.26	−0.032 157 720 549	22.700 493 507	−0.036 553 597 541	−0.036 553 60	17.230 927 430
	0.57	−0.034 940 841 538	20.400 019 294	−0.042 230 125 318	−0.042 230 13	16.420 194 663
	1.59	−0.063 945 693 111	11.042 617 509	−0.084 718 765 323	−0.084 718 77	8.806 674 0094
5p	0.0	−0.020 000 000 000	36.500 000 000	−0.020 000 000 000		36.500 000 000
	0.26	−0.020 360 792 830	36.096 483 134	−0.023 349 792 551		30.203 357 619
	0.57	−0.021 903 580 743	33.091 901 133	−0.025 258 633 309		28.971 093 917
	1.59	−0.033 489 056 655	21.958 625 107	−0.035 932 061 579		21.062 981 799
6p	0.00	−0.013 888 888 889	53.000 000 000	−0.013 888 888 889		53.000 000 000
	0.26	−0.014 070 782 545	52.498 386 898	−0.015 883 759 185		45.834 836 839
	0.57	−0.014 999 966 315	48.833 976 319	−0.016 767 326 295		44.081 940 239
	1.59	−0.020 903 914 137	35.390 733 125	−0.021 816 312 170		34.305 718 085
7p	0.00	−0.010 204 081 633	72.500 000 00	−0.010 204 081 633		72.500 000 000
	0.26	−0.010 309 583 975	71.904 566 27	−0.011 458 478 666		64.278 040 559
	0.57	−0.010 907 611 689	67.596 640 10	−0.011 950 076 576		62.104 185 020
	1.59	−0.014 343 119 913	51.737 977 258	−0.014 800 187 457		50.435 203 259

Table 3. Energy E and expectation value $\langle r \rangle$ of the hydrogen atom of the $3d, \ldots, 8d$ states in presence of a inverted Gaussian potential (2) centered at $r_c=2.50$ and 3.54 Å with $w_0=0.5$ a.u. as a function of σ. The free case corresponds to σ = 0. Here the values for E and $\langle r \rangle$ are displayed with 12 significant figures, and they correspond to the results obtained with the Lagrange Mesh Method. The coincidence with the FEM results is in 7–8 digits whilst for the FDM, the agreement with the results of the LMM is in all figures.

		$r_c=2.50$ Å		$r_c=3.54$ Å
nd	σÅ	E a.u.	$\langle r \rangle$ a.u.	E a.u.	$\langle r \rangle$ a.u.
3d	0.00	−0.055 555 555 555	10.500 000 000	−0.055 555 555 555	10.500 000 000
	0.26	−0.128 471 359 504	5.449 216 8184	−0.146 974 205 177	6.894 926 1986
	0.57	−0.251 571 681 618	4.986 228 7243	−0.269 235 578 841	6.727 277 2114
	1.59	−0.411 853 426 848	5.127 318 4744	−0.432 795 177 458	6.729 277 3865
4d	0.00	−0.031 250 000 000	21.000 000 000	−0.031 250 000 000	21.000 000 000
	0.26	−0.041 811 532 197	15.885 456 275	−0.038 490 481 788	17.852 091 651
	0.57	−0.044 189 228 465	14.896 944 804	−0.040 240 022 616	16.883 522 700
	1.59	−0.127 836 599 792	6.881 247 2011	−0.172 635 611 475	7.440 790 0868
5d	0.00	−0.020 000 000 000	34.500 000 000	−0.020 000 000 000	34.500 000 000
	0.26	−0.024 569 869 745	28.278 727 512	−0.022 880 682 620	30.857 778 724
	0.57	−0.025 611 985 197	27.021 828 122	−0.023 721 575 291	29.575 002 055
	1.59	−0.037 461 999 111	18.248 819 233	−0.039 108 470 145	17.049 986 326
6d	0.00	−0.020 000 000 000	34.500 000 000	−0.013 888 888 889	51.000 000 000
	0.26	−0.016 327 055 430	43.560 989 985	−0.015 360 981 734	46.717 692 906
	0.57	−0.016 886 116 520	42.025 765 284	−0.015 837 354 470	45.132 637 791
	1.59	−0.022 489 169 307	31.285 420 901	−0.023 341 455 972	29.854 962 236
7d	0.00	−0.010 204 081 633	70.500 000 00	−0.010 204 081 633	70.500 000 00
	0.26	−0.011 667 991 818	61.812 020 47	−0.011 065 988 156	65.536 819 33
	0.57	−0.012 004 198 432	59.997 674 38	−0.011 362 746 680	63.657 452 22
	1.59	−0.015 168 052 460	47.209 127 992	−0.015 654 825 275	45.536 405 971
8d	0.00	−0.007 812 500 000	93.000 000 00	−0.007 812 500 000	93.000 000 00
	0.26	−0.008 763 164 361	83.050 904 42	−0.008 363 446 047	87.340 123 04
	0.57	−0.008 981 421 027	80.957 563 96	−0.008 560 845 234	85.170 841 06
	1.59	−0.010 955 981 002	66.101 951 82	−0.011 258 294 889	64.170 188 267

Before discussing the concrete results let us briefly mention some details about the LMM. In the system (10), there are two free parameters: the size of the mesh N and the scaling factor h. The optimal values of these two parameters (N, h) depend on the considered state of the system as well as on the value of the parameters occurring in the Gaussian potential. In general, the results presented in tables 1–3 for states s, p and d, respectively, are obtained with a mesh of at least N = 250 points and h in some interval between $h\in [0.02-1.0]$. For states with small quantum numbers, even meshes of size n = 120 are adequate. For example, if we consider the state 1s with $w_0 = 0.5$, $r_c = 2.5$ Å and σ = 2.26 Å the energy is obtained by using $(N,h) = (120,0.04-0.12)$ while for the state 8d, for the same values of the parameters of the potential, $(N,h) = (250,0.26-0.30)$.

The convergence of the LMM is determined by the stability of the results with respect to an increase in the size of the base N and the variation of h. Table 1 presents the results for the $1s, \ldots, 6s$ states. For $r_c = 3.54$ Å a comparison with [35] is possible for the levels 1s, 2s, 3s, 4s with $\sigma = 1.59, 0.57$ and 0.26 Å where it can be seen a complete agreement in 7 decimal digits (for 1s and 2s states) and 8 decimal digits (for 3s and 4s). For completeness, the case σ = 0 (the free hydrogen atom) is also presented. For all these s-states, the energy decreases by increasing the value of σ as can be seen in figure 4(c). On the other hand, the presence of the Gaussian potential has an important effect on the expectation value of r as indicated in the seventh column of table 1 (see also figure 4(d)): at fixed $r_c = 3.54$ Å (i) as a function of $\sigma \in[0.00,1.59]$ Å, $\langle r \rangle$ increases for the ground state, whilst (ii) for the states $4s, \ldots, 6s$, $\langle r \rangle$ decreases. Columns 3 and 4 of table 1 (see also figures 4(a) and (b)) present the results of the energy and the expectation value $\langle r \rangle$ when the center of the Gaussian potential is $r_c = 2.50$ Å. As well as for $r_c = 3.54$ Å, the energy decreases as a function of σ. The behavior of the expectation value $\langle r \rangle$ is depicted in figure 4(b). This effect on $\langle r \rangle$ reflects how the electronic charge is attracted to the Gaussian part of the potential.

Zoom In Zoom Out Reset image size

States with l = 1 ($2p, \ldots, 7p$) are presented in table 2 for $r_c = 2.50$ and 3.54 Å. In both cases, with respect to the free configuration (σ = 0), the system gets more bound as the value of σ increases (see figures 5(a) and (c)). For $r_c = 3.54$ Å a comparison with [35] is also possible: we see a complete agreement in 7 decimal digits for the 2p state and 8 decimal digits for the 3p and 4p states for the three values of $\sigma = 0.26, 0.57$ and 1.59 Å. The expectation value $\langle r \rangle$ is shown in figures 5(b) and (d).

Zoom In Zoom Out Reset image size

Results of the energy and the expectation value $\langle r \rangle$ for l = 2 ($3d, \ldots, 8d$) are displayed in table 3 for $r_c = 2.50$ and 3.54 Å. For both values of r_c , by increasing σ the energy becomes more negative compared to the unconfined system (σ = 0). The expectation value $\langle r \rangle$ exhibits a decreasing behavior with the increasing of σ for all states except the lowest 3d state, which decreases and eventually increases.

3.1. Comparison of LMM with the finite difference and FEMs

The FDM was implemented in the second and fourth order in Matlab. Calculations using second-order finite differences require a very large number n of nodal points, which implies to diagonalize very large matrices, and therefore a significant computational time is involved. For this reason, we decided to use fourth-order finite differences, in which case the error in the solution of the eigensystem is of order h⁴, here h being the distance between two consecutive nodal points. For $\omega_0 = 0.5$ a.u., σ = 0.26 Å and $r_c = 3.54$ Å, by using h = 0.01, we obtained an accuracy of 9 decimal places in energy eigenvalues when they are compared with the results calculated with the LMM. For l = 0 and l = 1, $L_\textrm{max} = 160$, whereas for $l = 2, L_\textrm{max} = 200$. Analogous results are found for other values of $\omega_0, \sigma$, and $L_\textrm{max}$ and different values of l. For H@C₆₀ the step size h = 0.01 provides results such that the rounding error does not affect the first 11 decimal digits. Another quantity that can affect the accuracy of the eigenvalues is the quantum state to be calculated. For example, at fixed values $r_c = 3.54$ Å, $\omega_0 = 0.5$ a.u. and σ = 0.26 Å, in order to obtain eigenvalues with 8 decimal places of accuracy it is necessary that $L_\textrm{max} = 150$ a.u. and a grid with $n = 15\,000$ for states with l = 0, whereas for states with l = 1, $L_\textrm{max} = 170$ a.u. and $n = 17\,000$, and for l = 2 the optimal parameters are $L_\textrm{max} = 200$ a.u. and $n = 20\,000$. This shows that to calculate excited states, while maintaining the same accuracy of the previous energy levels, it is necessary to increase the value of $L_\textrm{max}$ as well as n. Additionally, by increasing the value of $L_\textrm{max}$ and n the agreement with the results of the LMM is in all digits.

When the FEM is applied, the agreement with the results of the LMM is in 8 decimal digits. It should be mentioned that the MATHEMATICA 12.3 software package allows us to calculate the solutions of the eigenvalue problem (4) as well. It can be easily done using the NDEigensystem command (based on the FEM) which, in general, provides the smallest eigenvalues and eigenfunctions of the involved linear differential operator on a certain finite region, $L_\textrm{max} \sim 150$. For the lowest states, in table 4 we display the relative difference $\frac{|E_M - E_{LM}| }{|E_{LM}|}$ between the energy E_M obtained with MATHEMATICA 12.3 and the corresponding value E_LM computed in the Lagrange-Mesh approach. Specifically, at fixed values of s-states (l = 0) with parameters $r_c = 3.54$, $\omega_0 = 0.5$, and σ = 1.59, by using a maximum cell measure of 0.001 the optimal value of $L_\textrm{max}$ increases monotonically from $L_\textrm{max} = 30$ at n = 1 ³ (ground state) up to $L_\textrm{max} = 120$ for n = 6. The calculations were run in MATHEMATICA 12.3 on a personal laptop.

Table 4. Relative error $|E_M - E_{LM}|/|E_{LM}|$ as a function of the quantum numbers n and l at $\omega_0=0.5$ a.u. Here E_M is the energy obtained with MATHEMATICA (finite element) whereas E_LM corresponds to the result of the LMM. Results are given by the number followed by the power of 10.

		$r_c=3.54$ Å			$r_c=2.5$ Å
n	σ Å	l = 0	l = 1	l = 2	l = 0	l = 1	l = 2
1	1.59	$1.4($−10)	$4.6($−10)	$5.7($−10)	$4.9($−10)	$2.6($−10)	$4.4($−10)
	0.57	$2.6($−10)	$7.8($−10)	$9.8($−10)	$2.5($−10)	$3.9($−10)	$6.4($−10)
	0.26	$2.4($−10)	$1($−9)	$1.7($−9)	$3.8($−10)	$6.6($−10)	$1.4($−9)
2	1.59	$5.7($−10)	$6.9($−10)	$1($−9)	$7.2($−10)	$7.4($−10)	$1.5($−9)
	0.57	$7.3($−10)	$1.8($−9)	$3.5($−9)	$4.5($−10)	$2.8($−9)	$3.5($−9)
	0.26	$1($−9)	$2.4($−9)	$3.5($−9)	$6.5($−10)	$2.9($−9)	$3.6($−9)
3	1.59	$5.5($−10)	$1.8($−9)	$3.5($−9)	$1.4($−9)	$2.6($−9)	$3.3($−9)
	0.57	$2.7($−9)	$3 ($−9)	$2.5($−9)	$2.7($−9)	$3.2($−9)	$2.7($−9)
	0.26	$2.9($−9)	$3 ($−9)	$2.7($−9)	$2.8($−9)	$2.9($−9)	$2.5($−9)
4	1.59	$2.6($−9)	$2.9($−9)	$2.6($−9)	$2.9($−9)	$2.9($−9)	$2.7($−9)
	0.57	$2.9($−9)	$2.5($−9)	$6.2($−9)	$2.9($−9)	$3.6($−9)	$5.7($−9)
	0.26	$2.5($−9)	$3($−9)	$6.8($−9)	$2.8($−9)	$4.3($−9)	$6($−9)

3.2. Numerical interpolation via an artificial neural network

In the Born-Oppenheimer approximation, for the lowest states with angular momentum $l = 0,1,2$ the energies and eigenfunctions of the hydrogen atom inside a C₃₆ and C₆₀ fullerenes are presented. The confining barrier was modeled by an inverted Gaussian function $w(r) = -\,\omega_0\, \textrm{exp}[-(r-r_c)^2/\sigma^2]$. The approximate solutions of the corresponding Schrödinger equation were determined by three different numerical methods: (i) the LMM, (ii) the (fourth-order) finite difference and (iii) the FEM. As a complementary tool, we use an artificial neural network to interpolate/extrapolate the results. In general, the confined hydrogen atom becomes more bound with the increase of the width σ and the depth ω₀ of the Gaussian potential, for the two considered cases of r_c .

Using the LMM accurate energies with not less than 11 significant figures were obtained. The optimal values of the size of the mesh N and the scaling factor h depend on the state being studied as well as on the parameters of the confining potential w(r). Since this method is not completely based on a variational principle, we must point out that the computed energies are not necessarily greater than or equal to the (unknown) exact ones. However, in all known cases where the Lagrange mesh results are stable with respect to the variation of N and h, like in the present study, it turns out that it converges rapidly and generates simple highly accurate solutions.

The FDM is a suitable method for solving the radial Schrödinger equation (14). In the present work, we used a fourth-degree approximation for the kinetic energy operator (18). The accuracy of the eigenvalues depends strongly on the corresponding step size h (equation (19)). In particular, using a fourth-degree approximation the error in the kinetic energy operator (equation (18)) is of order h⁴. The calculations were performed in a code written with double precision in MATLAB. It should also be noted that this method, like the LMM, is not based on the variational principle.

In the case of the FEM, precise energies (always from above the exact ones) can be calculated by simultaneously increasing $L_\textrm{max}$, the number of elements in which the interval $[0,\,L_\textrm{max}]$ is divided and the degree of the interpolating polynomials. Moreover, this method allows us to deal with problems in higher spatial dimensions with regular and irregular boundaries, which is not so easy to implement in the Lagrange-mesh and FDMs.

Finally, the artificial neural network is computationally a faster efficient tool to compute the spectra. Nevertheless, for the training stage, it requires to know in advance accurate results in several points on the space of parameters. In combination with the Lagrange-mesh or the FDM it significantly reduces the overall computational time, although with less accuracy. By means of the methods used in the present study, energies were obtained with higher accuracy than those reported in the literature. In future works, we plan to establish in a more systematic manner the implementation of artificial neural networks in the study of other physically relevant atomic and molecular systems.

No new data were created or analysed in this study.

The authors are grateful to S A Cruz for their interest in the work and useful discussions. M A Q-J would like to thank the support from DGAPA-UNAM under Project UNAM-PAPIIT TA101023.