Assessing the Global Natural Orbital Functional Approximation on Model Systems with Strong Correlation

,


I. INTRODUCTION
Since the electron-electron interaction is described by a two-body operator in the non-relativistic Born-Oppenheimer approximation, the resulting energy is an exact functional of the second-order reduced density matrix (2RDM).Although the existence of the functional in terms of the first-order RDM (1RDM) was proven [1], its specific form remains unknown.The adopted strategy for electronic calculations involves the reconstruction of the 2RDM in terms of the 1RDM.However, the energy obtained through this process is not the exact functional of the 1RDM.Consequently, an approximate functional in this manner retains a hidden dependence on the 2RDM, leading to the consequential functional Nrepresentability problem [2,3].The question arises: How can we express the 2RDM in terms of the 1RDM?
In 1956, Löwdin and Shull demonstrated [4] that by constraining the 1RDM to the diagonal representation, we can express the two-electron wavefunction in the simplest manner.In this representation, the diagonal elements of the 1RDM and their corresponding eigenvectors are denoted as occupation numbers (ONs) and natural orbitals (NOs), respectively.The most straightforward generalization of two-electron functional form to N-particle systems is to divide the system into independent electron pairs that leads to the Piris natural orbital functional 5 (PNOF5) [5].An interesting relationship between geminal approaches and NOF theory emerges here, as PNOF5 has been demonstrated to be equivalent to a specific case of an antisymmetrized product of strongly orthogonal geminals (APSG) [6,7].Additionaly, PNOF7 [8,9] mantains PNOF5 intra-pair interactions but includes inter-pair correlation with an antisymmetrized geminal power (AGP) like term.
Despite recent advances obtained in the context of geminal wavefunction approaches [10] and Richardson-Gaudin states [11], including dynamic correlation effects continue being a challenge.The same arguments motivated the development of a perturbative correction to PNOF5 [12].In this line, hybrid approaches have been recently extensively explored in NOF theory [13][14][15][16][17].For instance, Wenna et al. [15] combined DFT and NOF with the aim to exploit the former in the equilibrium region and the latter for long-range interactions.More recently, Gibney and co-workers [17] introduced a damping factor to avoid double counting in a functional with terms from both DFT and NOF approximations.Nevertheless, the access to correlated density matrices is limited in this context and thereby it is preferable to have a fully correlated energy functional.A global NOF (GNOF) has been recently proposed [18] as a method to deal with global electronic structure problems.Despite it is an electron-pairing approach, GNOF is not an independent-pair approximation anymore.This functional introduces interactions between electron pairs to treat in a balanced manner the static and dynamic electron correlations [19][20][21][22][23]. Therefore, in contrast to its predecessors [8,9,12,24], GNOF does not require perturbative corrections to include dynamic correlation.
The motivation that brings us to this study is the assessment of GNOF.Although NOF approximations have been proposed for two decades now, benchmarking new functionals is essential to establish them among daily used methods.The set of systems selected for this purpose has been recently employed by Gaikwad and co-workers in Ref. [10].A few of these systems revealed in Ref. [19] the ability of GNOF to describe H model systems of strong correlation, which can be computed at the Full-Configuration-Interaction (FCI) level making them appropriate for comparisons.A more complete and challenge set is included in the present work, where we focus on the performance of the GNOF approximation on energy curves.This article is structured as follows: The first section (II) emphasizes the key aspects of the GNOF approximation for spin-compensated systems.The subsequent section (III) analyzes a set of model systems, specifically conducting a comparison between NOF calculations and Full Configuration Interaction (FCI).Finally, conclusions are drawn in section IV.

II. NATURAL ORBITAL FUNCTIONAL APPROXIMATIONS
In this section, we provide a brief description of GNOF for singlet systems.A comprehensive review of recent developments in approximate NOFs has just been published [25].GNOF was proposed as a versatile method for addressing various electronic structure problems [18].It has demonstrated accurate comparisons with experiments and state-ofthe-art calculations for energy differences between ground states and low-lying excited states with different spin states of atoms ranging from H to Ne, ionization potentials of first-row transition-metal atoms (Sc-Zn), total energies of a selected set of 55 molecular systems in different spin states, homolytic dissociation of selected diatomic molecules in different spin states, and the rotation barrier of ethylene.More recently, GNOF has been successfully employed to reduce the delocalization error in NOF approximations [21], conduct ab-initio molecular dynamics simulations [26], develop excited state methods by coupling GNOF with extended random phase approximation [27], and, most importantly for our purposes, investigate problems with high nondynamic correlation such as the aromaticity of the Al − 3 anion [23] and Iron(II) Porphyrin [20], which poses a challenge due to its computational demands.Hereafter, nondynamic and static terms are used interchangeably to denote strong correlation.
The most common strategy for developing NOFs is to approximate the 2RDM in terms of the ONs and NOs.This reconstruction is not applicable in another representation of the 1RDM matrix but in its diagonal form.Therefore, it is more precise to refer to NOF theory rather than 1RDM theory.The exact non-relativistic energy functional of the 2RDM is employed, leading to the solution being established by optimizing the energy with respect to the ONs and the NOs separately.As a result, orbitals vary throughout the optimization process until the most favorable orbital interactions are attained.
Hereafter we focus on spin compensated systems (S = 0).Consequently, the spin-restricted theory can be applied, ensuring that all spatial orbitals are doubly occupied with equal occupancies for particles with α and β spins: It is important to note that our NOFs preserve the total spin [28] regardless the value of S. Electron pair models require to divide the orbital space Ω.We divide the latter into N Ω = N/2 mutually disjoint subspaces Ω g , each of which contains one orbital |g with g ≤ N Ω , and N g orbitals |p with p > N Ω , namely, Taking into account the spin, the total occupancy for a given subspace Ω g is 2, which is reflected in the following sum rule: Here, the notation p ∈ Ω g represents all the indexes of |p orbitals belonging to Ω g .In general, N g may be different for each subspace.In this work, N g is equal to a fixed number for all subspaces Ω g ∈ Ω.We adopt the maximum possible value of N g which is determined by the number of basis functions (N B ) used in calculations.Indeed, variable N g will play a key role in the Be + H 2 reaction studied in subsection III C. Taking into account Eq. ( 3), the trace of the 1RDM is verified equal to the number of electrons: Using ensemble N-representability conditions [29], we can generate a reconstruction functional for the 2RDM in terms of the ONs that leads to GNOF [18]: For singlet states, the intra-pair component is formed by the sum of the energies of the pairs of electrons with opposite spins, namely ) where Π(n p , n r ) = c(n p )c(n r ) and c(n p ) is defined by the square root of the ONs according to the following rule: that is, the phase factor of c p is chosen to be +1 for a strongly occupied orbital, and −1 otherwise.H pp are the diagonal oneelectron matrix elements of the kinetic energy and external potential operators and L pq = pp|qq are the exchange-timeinversion integrals [30].The inter-pair Hartree-Fock (HF) term is where J pq = pq|pq and K pq = pq|qp are the Coulomb and exchange integrals, respectively.The prime in the summation indicates that only the inter-subspace terms are taking into account (p ∈ Ω f , q ∈ Ω g , f = g).The inter-pair static component is written as where Φ p = n p h p with the hole h p = 1 − n p .Note that Φ p has significant values only when the ON n p differs substantially from 1 and 0. In Eq. ( 9), Ω b denotes the subspace composed of orbitals below the level N Ω (p ≤ N Ω ), so interactions between orbitals belonging to Ω b are excluded from E inter sta .Finally, the inter-pair dynamic energy can be conveniently expressed as The dynamic part of the ON is defined as with h c = 0.02 √ 2 [18].The maximum value of n d p is around 0.012 in accordance with the Pulay's criterion that establishes an occupancy deviation of approximately 0.01 with respect to 1 or 0 for a NO to contribute to the dynamic correlation.Considering real spatial orbitals (L pq = K pq ) and n p ≈ n d p , it is not difficult to verify that the terms proportional to the product of the ONs will cancel out, so that only those terms proportional to Π will contribute significantly to the energy.

III. RESULTS
The systems investigated in this study can be categorized into three groups.The first and second groups both involve the H 4 model, while the last group includes the symmetric dissociation of H 8 in one and three dimensions and the insertion of a Be atom into an H-H bond.All calculations were performed using the DoNOF code [31] with the STO-6G basis set [32], unless otherwise specified.Full Configuration Interaction (FCI) results were obtained from [10], except for those related to the Be + H 2 reaction, for which the Be(3s2p)/H(2s) basis was extracted from Ref. [33].The planar H 4 system was previously investigated in the context of NOF theory by Ramos et al. [34].They concluded that nondynamic inter-pair correlation is crucial for accurately describing the D 2h − D 4h transition.Consequently, the independent-pair model PNOF5 inaccurately exhibited a cusp along this transition.In the following sections, we examine whether GNOF can rectify the deficiencies of PNOF5 in the rectangular H 4 system.Additionally, for comprehensive analysis, we include additional geometries and variants, referred to as planar systems, as depicted in Fig. 1.The notation follows that of Ref. [10], although these isomers of H 4 were originally introduced by Paldus and co-workers [35].
GNOF absolute energies can be observed in the top plot of Fig. 2 for the planar D4, S4, and P4 configurations.Qualitatively, all the energy curves exhibit similar behavior: energies decrease as the distance between two hydrogen molecules increases in the case of D4 and P4, and similarly for S4, but now with respect to inter-nuclear distances.However, it is worth noting a slight shoulder in the green curve corresponding to the rectangular P4 system.This corresponds to α = 2 Bohr, the point at which P4 converges to the S4 structure.
The bottom plot in Fig. 2 illustrates GNOF energy errors for the D4, S4, and P4 model systems.The former corresponds to a one-dimensional problem, where two hydrogen molecules are dissociated.Since GNOF is a quasi-exact FIG. 3. GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for the planar H4 system with different a distances (in Bohr).method for two-electron systems, the error approaches zero as the correlation between two molecules decreases.According to Fig. 2, the maximum error occurs slightly above α = 1.0 Bohr, where one orbital is practically doubly occupied and the other strongly occupied orbital yields an occupancy of 1.7, which represents a typical intermediate case between dynamical and static electron correlation effects.However, the maximum GNOF error remains below 0.01 Hartree.It is noteworthy to observe the alteration in bonding patterns as D4 transitions from α < a to α > a.This change is discerned through GNOF in the occupation scheme.Beyond α = 2.0 Bohr, the preference shifts towards two strongly occupied quasidegenerate NOs and two weakly occupied quasi-degenerate NOs.
As observed in the case of the Hubbard model [36], the error increases as we transition from linear to two-dimensional electronic systems.In both the square and rectangular H 4 systems, the agreement with FCI is remarkable beyond a = α = 4.0 Bohr.Note that both geometries converge to the same system at a = α = 2.0 Bohr.However, GNOF exhibits larger errors for the square S4 compared to the rectangular P4 geometry at both intermediate and short distances between hydrogen molecules.Interestingly, in the case of S4, GNOF yields two spatial orbitals that are singly occupied (2n 2 = 2n 3 = 1) along the entire dissociation curve.This may be associated with an overestimation of the static correlation at short and intermediate distances.In this system, a minimal basis set that can only provide a description of bonding-antibonding orbitals appears insufficient to capture the correlation effects adequately.In-cluding more weakly occupied orbitals by using a larger basis set seems to be a better option, as we will demonstrate in subsection III C for the BeH 2 system, where we encounter a similar problem.
On the other hand, the rectangular P4 system is more analogous to the H 4 model previously studied in the context of NOF approximations [34].GNOF performs quantitatively accurately here compared to FCI at very short distances and α ≥ 3.0 Bohr.At the dissociation limit, two ONs of strongly and weakly occupied orbitals degenerate, specifically with 2n 1 = 2n 2 = 1.925 and 2n 3 = 2n 4 = 0.075, respectively, at α = 8 Bohr.However, the electronic configuration changes to four distinct doubly occupied orbitals as α shortens.For instance, at α = 1 Bohr, the occupancies are equal to 1.993, 1.980, 0.020, and 0.007.The intermediate region where both configurations compete is also well-handled by GNOF, with the error not exceeding 0.045 Hartree.It is worth noting that P4 is a typical example wherein GNOF can find multiple solutions, making it challenging to consistently adhere to the minimum energy curve.As pointed out in Ref. [11], it is in this region around α = 2 Bohr where the dominant Slater determinant in the FCI expansion changes due to a switching of the long and short sides of the rectangle.
To extend our analysis of planar model systems, we investigate the trapezoidal H 4 configuration with distances between successive atoms fixed to a = 1.2, a = 1.6, and a = 2.0 Bohr.As illustrated in the geometry shown in Fig. 1, when δ = 0, the S4 geometry is retrieved, and the D4 structure is obtained when δ = π/2.The top and bottom plots in Fig. 3 display absolute and error energies corresponding to this system.Overall, the absolute energy curves resemble those obtained for the planar D4, S4, and P4 systems shown in Fig. 2. It is important to note, however, that Fig. 3 presents energies with respect to the parameter α, which causes the structure to vary between S4 and D4 as described above.Hence, for any a distance between hydrogen molecules, the minimum energy is obtained for the linear geometry and the maximum energy for the square configuration.In comparison with the planar models studied above, the performance of GNOF remains consistent for trapezoidal systems, as evident from the difference energy curves shown in FIG. 5. GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for the non-planar V4 system.The distance a is fixed to 2.0 a.u. in both plots.the bottom plot of Fig. 3. Johnson and DePrince pointed out [11] that trapezoidal H 4 systems require a description of weak correlation effects not present in the aforementioned S4 and P4 systems.Despite posing a challenge for the orbitaloptimized doubly-occupied configuration interaction method, these systems are well-described by GNOF.
We observe small deviations in energy below FCI at α ≥ 0.1 for the a = 1.2 Bohr case (bottom plot of Fig. 2).Although details regarding N-representability are omitted in Section II, it is worth noting that the reconstructed 2RDM that leads to GNOF only satisfies some necessary Nrepresentability conditions.Consequently, violations of functional N-representability [2,3] may occur, leading GNOF to result in energies below FCI values.On the other hand, the 1RDM is an ensemble N-representable matrix since the ONs belong to the interval [0,1].Additionally, it also satisfies the pairing conditions (3) linked to pure N-representability [3], allowing the avoidance of obtaining fractional charges in homolytic dissociations [37,38].

B. Non-Planar H 4 benchmark set
The amount of correlation energy is strongly related to the dimensionality of electronic systems.This phenomenon was demonstrated in the context of NOF approximations by studying the Hubbard model in one [39] and two dimensions [36].Here, we extend the benchmark conducted in the previous section to three-dimensional H 4 configurations.The geometries FIG. 6. GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for the non-planar T4 system.The distance a is fixed to 2.0 a.u. in both plots.R values are given in a.u..  4, where the notation is again taken from Ref. [10].
As shown in Fig. 4, the H 4 systems studied in this section can be obtained by twisting the rectangular P4 model shown in the previous section.In the case of T4, the angle θ varies in the range [0, π/2] for many fixed distances between the two centers of hydrogen molecules, whereas for V4, the angle between two bonds is set to π/2 and it is the distance which increases from α = 1.0 to α = 10.0 Bohr.For the latter, GNOF gives energies similar to the rectangular P4 shown in the bottom panel of Fig. 2. Fig. 5 shows GNOF absolute energies and errors with respect to FCI for V4 in terms of the distance between two hydrogens specified above.While at intermediate distance (α ≈ 1.5 Bohr) the error approaches 0.04 Hartree, it rapidly decreases at both shorter and longer distances.As in many other cases, the orbital scheme is different in these regions.GNOF gives two singly occupied spatial orbitals with half ONs (2n 3 = 2n 4 ≈ 1.0) at short distances, but then it changes smoothly to an HF-like configuration with two markedly occupied orbitals and corresponding weakly occupied orbitals at long distances.Regarding absolute energies, the curve in the top plot of Fig. 5 seems similar to the error since, after a slight increase in the energy from α ≈ 1.1 to α ≈ 1.5 Bohr, it rapidly decreases to a constant beyond α ≈ 4 Bohr.
As illustrated by Fig. 6, for T4, the situation changes dramatically depending on the fixed distance between two hydrogen bonds.When the two hydrogen molecules are brought as close as R = 1.1428Bohr, the GNOF error stays under 0.04 Hartree for small angles up to approximately 0.45π.However, the performance rapidly deteriorates from this point to θ = π/2, where we obtain a difference E GNOF − E FCI close to 0.1 Hartree.According to the top plot of Fig. 6 corresponding to absolute energies, as the energy increases and the red curve ascends, the error with respect to FCI also increases, indicating that GNOF is not able to capture a significant portion of dynamic electron correlation in this region.Natural ONs are close to integer values along the full red curve, so a larger basis set allowing a larger N g value is necessary to improve GNOF performance, as happened for BeH 2 in subsection III C.
Similar to the behavior shown by geminal wavefunctions developed by Gaikwad and co-workers in Ref. [10], after reaching some distance close to R = 2 Bohr, the difference with respect to FCI remains nearly constant for any angle between two hydrogen bonds, as illustrated by the solid curve in Fig. 5 corresponding to R = 4 Bohr.There are no significant changes in the ONs and NOs when θ varies between 0 and π/2, so GNOF exhibits two quasi-degenerate orbitals with ONs close to 0.98 and corresponding weakly occupied orbitals with ONs ≈ 0.005.This suggests that the electronic problem is dominated by dynamic correlation.
Results corresponding to the T4 geometry with distances R = √ 2, 2, 3, 7 a.u.can be found in the supplementary material.

C. Symmetric dissociation of linear and cubic H 8 and the Be+H 2 reaction
Following the work by Gaikwad et al. [10], we conclude the assessment of GNOF with three paradigmatic cases: the dissociation of a 3D H 8 cube, a linear chain of eight hydrogen atoms, and the insertion of a beryllium atom into a hydrogen molecule.Similar systems to the former were studied in Ref. [19].Despite the qualitatively correct dissociation curve provided by GNOF in Ref. [19] for H 6 4 in three dimensions, the lack of a reliable reference method motivates us to include a similar system such as the cubic H 8 in the present work.Geometries and notation corresponding to the H 8 systems are shown in Fig. 7.Note that symmetric dissociation refers to stretching all bonds simultaneously to obtain, in the end, 8 isolated hydrogen atoms.When the correlation regime is increased, a metal-to-insulator transition is observed in both systems.Accordingly, they are paradigmatic models for strongly correlated Mott insulators.
Figs. 8 and 9 show absolute energies and energy differences with respect to FCI given by GNOF for the cubic and linear H 8 systems, respectively.As expected, at large distances in the strong correlation regime, the error rapidly approaches zero in both model systems.However, as illustrated in the bottom plot of Fig. 8, the performance at intermediate and short distances breaks down in the three-dimensional case.Let us focus on the latter, since here the error increases up to one order of magnitude larger than at intermediate distances.This system resembles the S4 model studied in Fig.  nately, in these systems, GNOF struggles to retrieve dynamic correlation at very short distances with a minimal basis set, so it is necessary to include more weakly occupied orbitals to improve GNOF energies in this region.
Regarding H 8 dissociation in one dimension, Fig. 9 confirms the conclusions obtained in Ref. [39] for GNOF predecessors.At large distances, GNOF retrieves FCI energies.However, GNOF misses some correlation energy at intermediate H-H distances; for instance, there is a difference of approximately 0.03 Hartree at R H-H = 2.5 a.u., where two strongly occupied orbitals yield 2n 1 = 2n 2 ≈ 1.80.It is worth noting that GNOF goes below FCI when hydrogen atoms are brought very close to each other.Overall, plots corresponding to absolute energies show qualitatively good behavior for both systems and resemble typical dissociation curves in molecular systems.A comparison between absolute energy and en-ergy error plots reveals that while equilibrium and dissociation limit regions are accurately described by GNOF, intermediate and short bond distance regions are especially challenging.
BeH 2 is employed here to represent a chemical reaction involving strong correlation.Recently, Ammar and co-workers [33] have studied the ability of the GW approximation to describe the multireference region of the Be+H 2 reaction.In fact, when a Be atom is inserted into H 2 following a C 2v pathway as illustrated in Fig. 10, the FCI wavefunction of the whole system is dominated by two different electronic configurations depending on the distance between beryllium and H 2 center.Hence, the switching region is plagued by multireference effects.It is worth noting that different basis sets were employed in Refs.[33] and [10] to study the strong correlation of the Be+H 2 reaction.While STO-6G was used in the latter, the former utilized the Be(3s2p)/H(2s) configuration.In Fig. 11, we show GNOF absolute energies and energy differences with respect to FCI along the beryllium insertion into H 2 , respectively in the top and bottom panels.In Ref. [33], Ammar et al. point out that in the region 2.5a 0 < x < 3a 0 , the FCI wavefunction switches from one electronic configuration to another.Therefore, it is at x = 2.75a 0 where both dominant configurations have the same weight, and maximum strong correlation effects might occur.Looking at absolute energies in the top plot of Fig. 11, there is an increase in energy in this region, which is significantly weaker if the basis set is large enough, as revealed by the red curve in comparison with the green and blue curves.In other words, GNOF has no problem describing this region as long as the orbital subspace containing an electron pair is large FIG.10.Structure corresponding to the Be+H 2 reaction based on Fig. 1 from Ref. [33].Note beryllium is placed at the coordinate origin, whereas hydrogens are at ±y = 2.54 − 0.46x varying x from 0 to 4 Bohr.
Be H H x y enough, as revealed by the energy errors in the bottom panel of Fig. 11.This issue is evident in the bottom plot of Fig. 11 when we focus on the red curve, which keeps the difference with respect to FCI below 0.01 hartree.
It is worth recalling here the meaning of the parameter N g = 1, which corresponds to a GNOF calculation with a minimal basis set, i.e., two orbitals per electron pair.According to the description given in section II, this means having one strongly occupied orbital (n p ≥ 0.5) coupled with a weakly occupied one (n p ≤ 0.5).Looking at the ONs obtained at x = 2.75a 0 , we have two orbitals almost fully occupied (n 1 ≈ 0.9998, n 2 ≈ 0.9825), two orbitals half-occupied (n 3 ≈ 0.5, n 4 ≈ 0.5), and two almost empty orbitals (n 5 ≈ 0.035, n 6 ≈ 0.0004).Nevertheless, when the orbital subspace is extended to four orbitals and thereby each strongly occupied NO is coupled with three weakly occupied NOs (N g = 3), we have a similar distribution of the ONs, except for the weakly occupied orbitals coupled to orbital number 2, which present a nonnegligible distribution of the ONs.This makes the difference with respect to FCI nearly zero.More details about the ONs corresponding to this system can be found in the supplementary material.

IV. CONCLUSIONS
This paper examines the performance of the Global Natural Orbital Functional (GNOF) approximation when applied to hydrogen atom clusters in one, two, and three dimensions, comparing its outcomes with those obtained from the Full Configuration Interaction (FCI) method.The GNOF energy expression, relying solely on two-index electron integrals, minimizes computational scaling with system size while maintaining the capacity to describe diverse types of electron correlation.This feature renders it highly attractive as an approximation method.
In this study, GNOF has demonstrated its effectiveness in capturing strong correlations across various scenarios.Three distinct blocks were distinguished, with the first and second corresponding to known groups of H 4 atoms.In these cases, GNOF errors remain below 0.05 Hartree, with the intermediate region showing maximum differences when dynamic and static electron correlation effects compete.Notably, challenges arose in the T 4 system, where two hydrogen molecules rotate with respect to each other, particularly when the distance between H 2 is small.
The third block focused on more intricate systems, including the symmetric dissociation of two H 8 systems.In one dimension, GNOF yielded some energies below those of FCI at very short internuclear distances, yet the overall accuracy remained high.However, during the dissociation of the H 8 cube, significant deviations between GNOF and FCI were observed.While agreement with FCI was accurate at bond H-H distances exceeding 1 Bohr, the error rapidly increased to approximately 0.55 Hartree as the internuclear H-H distance decreased to 0.5 Bohr.The assessment culminated with a typical reaction demonstrating strong multireference effects: the insertion of a beryllium atom into a hydrogen molecule.This reaction underscored the importance of expanding the number of orbitals within each orbital subspace.Indeed, GNOF achieved greater precision with an augmented basis set, transitioning from STO-6G to Be(3s2p)/H(2s).
In summary, the challenges encountered by GNOF in the T4 system, the dissociation of the H 8 cube, and the Be+H 2 reaction underscore the difficulty of capturing dynamical correlation at extremely short distances using minimal basis sets.This highlights the necessity of incorporating more weakly occupied orbitals to enhance GNOF energies in such scenarios.Coupled with recent publications endorsing the efficacy of GNOF across various atomic and molecular systems, this study reinforces GNOF as a well-rounded method for addressing global electronic structure problems.

FIG. 4 .
FIG.4.Structure of the non-planar H 4 model systems studied in this work.Labels set according to Ref.[10].
involved are illustrated in Fig.

FIG. 9 .
FIG.9.GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for the linear H 8 .

FIG. 11 .
FIG.11.GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for Be+H 2 reaction.Note legend to distinguish between STO-6G and Be(3s2p)/H(2s) basis sets, as well as N g values.
[10]1.Structure of the planar H 4 model systems studied in this work.Labels are set according to Ref.[10].FIG.2.GNOF absolute energies in a.u.(top) and GNOF energy differences with respect to FCI in a.u.(bottom) for the planar D4, S4 and P4 systems.Note that hydrogen bonds are fixed to a = 2.0 Bohr in both D4 and P4.