Computer program ATOM-MOL-nonBO for performing calculations of ground and excited states of atoms and molecules without assuming the Born–Oppenheimer approximation using all-particle complex explicitly correlated Gaussian functions

In recent years, there has been increasing interest in quantum mechanical calculations of molecular bound ground and excited states without assuming the Born–Oppenheimer (BO) approximation.^1–3 There are at least two major reasons why such calculations are interesting and can provide useful information. The first reason is related to higher accuracy in calculating the molecular rovibrational and electronic spectra, as well as the electronic spectra of atoms. When the BO approximation is not assumed, the wave functions and the corresponding nonrelativistic energy levels of the atom or molecule explicitly include effects originating from the finite mass of the nuclei and from the coupling of the motions of the nuclei and electrons. If such non-BO wave functions are then used to calculate relativistic and quantum electrodynamics (QED) corrections, these corrections will directly include the finite-nuclear-mass (FNM) effects, i.e., the so-called recoil effects. We have showed that such an approach can produce results whose accuracy match the accuracy of the most accurate experimental measurements.³ The second reason for carrying out the non-BO atomic and molecular calculations stems from the interest in describing properties and structures of molecules and atoms by calculations where all particles forming the system are treated on equal footing. With that, the molecular structure, dipole moment, polarizabilities, etc., are obtained as expectation values of operators representing these properties, and for molecules containing identical nuclei, the indistinguishability principle leads to interesting effects not present when the BO approximation is assumed.

Despite the fact that both BO and non-BO approaches are based on the principles of quantum mechanics, there are significant differences between them when it comes to calculating the energies and the corresponding wave functions. In order to describe a molecular rovibrational spectrum within the BO approximation, one needs to perform separate calculations of the electronic wave functions and the corresponding energies at some selected configurations of the nuclei placed in different fixed positions in space. These calculations yield the so-called potential energy surface (PES), which is used in the subsequent calculation of bound rovibrational states of the molecule. In non-BO calculations, as the nuclei and the electrons forming the molecule are treated on an equal footing, the calculations yield the total energies and the corresponding total wave functions, which explicitly depend on the coordinates of both the nuclei and the electrons. Both the energies and wave functions directly include all effects (including the high-order ones) that originate from the coupling of the motion of the nuclei and electrons, i.e., adiabatic and non-adiabatic effects.

At the nonrelativistic level, the Hamiltonian used to calculate the internal states of an atomic or molecular system commutes with the operators representing the square of the total orbital angular momentum operator and its projection on some axis (usually chosen to be the z-axis). Therefore, the corresponding quantum numbers, which can be used to label the eigenstates of the Hamiltonian, are good quantum numbers. As the Hamiltonian is isotropic (i.e., invariant with respect to rotations about the center of the internal coordinate system described later in this work), the wave functions are atom-like [i.e., they transform according to the irreducible representations of the SO(3) rotation group]. Thus, if basis functions that are eigenfunctions of the square of the total orbital angular momentum and its projection on the z-axis are used in the calculations, the problem of calculating atomic or molecular bound states separates into independent calculations that are performed for groups of states, each corresponding to a particular set of the total orbital angular momentum quantum numbers. This applies equally to atoms and molecules, as the non-BO Hamiltonians for both types of systems have the same spherical symmetry.

The starting point of the approach we have developed for non-BO atomic and molecular calculations is the transformation of the nonrelativistic Hamiltonian, which depends on laboratory-frame Cartesian coordinates of all particles forming the system. This Hamiltonian describes both the relative motion of the particles around the center of mass of the system (here, we termed it the internal motion) and the translational motion of the center of mass in 3D space. The latter corresponds to a free motion of the system as a whole, which has continuous (i.e., non-discrete) energy spectrum, and thus has to be separated. After this separation, the calculation can only focus on system’s bound states, which correspond to the internal motion. The internal bound states are eigenstates of the Hamiltonian (we call it the internal Hamiltonian), which is obtained by separating out the operator representing the motion of the center of mass from the laboratory Hamiltonian. The internal Hamiltonian used in our atomic and molecular non-BO calculations is given in Sec. II.

One of the central issues in non-BO calculations is the selection of appropriate basis functions for expanding the spatial part of the wave function. In this selection, several points have to be considered. They are related to the description of the interparticle correlation effects and to the description of the radial and angular nodes in the wave functions of the ground and excited states.

If the nuclei and electrons in a non-BO calculation are treated on equal footing, the nucleus–nucleus (n–n), electron–electron (e–e), and nucleus–electron (n–e) correlation effects need to be represented in the wave function of system. A natural and very effective form of this representation, at least for the case of pairwise interaction between particles, involves expanding the wave function in terms of basis functions, which explicitly depend on the inter-particle distances. Such functions are called explicitly correlated functions. Various types of explicitly correlated functions have been adopted in quantum mechanical calculations of atomic and molecular systems within the BO approximation. Most commonly, these functions are taken in some form of multi-particle Gaussian functions.

Atomic and molecular calculations where the BO approximation is not assumed represent a more recent development in the field of quantum molecular science. Such calculations require some modifications of the explicitly correlated basis functions, as the correlation effects, i.e., the e–e, n–e, and n–n (n–n is only relevant for molecules) correlation effects, that need to be described are more intricate than in the calculations concerning purely electronic states. As electrons are lighter and more delocalized in space, their wave functions (if we could talk about independent single-particle orbitals) overlap more strongly than the wave functions of nuclei. Thus, the correlation in the motion of electrons, on the relative scale, is smaller than the correlation in the motion of nuclei.

As the e–e, n–e, and n–n correlation effects are due to the electrostatic interactions between the particles, they depend on the charges of the particles. However, they also depend on particles’ masses. The mass dependence appears due to the following effect. For the light electrons, there is non-negligible probability of finding two of them (with opposite spins) at the same point in space. On the other hand, due to significantly larger masses of the nuclei, the internuclear correlation is stronger because the nuclei in a molecule avoid each other almost completely, even when their spins are opposite. The third type of correlation, the n–e correlation, can be called anti-correlation, as it is arises due to attraction and describes an effect of electrons, particularly the core electrons, following the nuclei very closely. The direct dependence of the basis functions on the electron–nucleus distances is key to accurately describing the latter correlation.

The explicitly correlated Gaussian functions that have the dependency on the interparticle distance only in the exponential factor are usually quite effective in describing the n–e and e–e correlation effects. However, as the probability of finding two nuclei at the same spatial point is several orders of magnitude smaller than finding two electrons in the same point, additional correlation factors need to be included in the explicitly correlated Gaussians for non-BO calculations. These additional factors are also needed to describe the radial and angular nodes of the wave function that appear due to ro-vibrational (and also electronic) excitations in molecular non-BO wave functions and due to electronic excitations in atomic non-BO wave functions. The need for such factors is due to Gaussians at coalescent points (r_ij → 0) having maxima, which are natural and desirable in describing the n–e correlation, but not desirable in describing the n–n and, partially, the e–e correlations. The factors can be chosen in different forms, but their role is to lower (or even zeroing) the amplitude of the wave function at r_ij = 0. The factors should also allow us to generate radial oscillations of the wave functions for “vibrationally excited” states. The angular oscillations of the wave functions, which are also needed to describe atomic and molecular excited states, can be achieved by multiplication of the Gaussians by appropriate Cartesian spherical harmonics. Finally, it is desirable to use a common set of basis functions in both atomic and molecular non-BO calculations, as the general structures and the symmetries of the Hamiltonians for these two types of systems are the same. Such a common basis set that can be used in very accurate atomic and molecular non-BO calculations is implemented in the computer code, which is described in this work.

Non-BO calculations of ground and excited states of atomic and molecular systems are interesting and useful for the following reasons: (i) treating nuclei and electrons of the system on equal footing is approximately equivalent to computing the entire molecular potential energy surface and determining rovibrational states of the system in a single calculation, (ii) the non-BO approach enables a direct account of the coupling of the vibrational and electronic degrees of freedom and can provide a very accurate description of this effect, (iii) non-BO calculations when performed with extended well optimized basis sets and with the inclusion of leading relativistic and quantum electrodynamics (QED) effects are capable of predicting spectral transitions with accuracy that matches or even exceeds the accuracy of the most accurate high-resolution experiments, and (iv) non-BO calculations provide a unique description of chemical and physical properties of atoms and molecules outside the realm of the BO approximation.

The non-BO effects can be determined directly by an approach such as the one presented in this work or by calculating these effects using the perturbation theory, where the BO solution is taken to be the zeroth-order approximation. A question that one can asks is which of the two approaches is more accurate and convenient to apply. An example that sheds some light on this point can be found in the paper by Pachucki and Komasa⁴ where one finds a comparison of the non-adiabatic corrections to the adiabatic energy levels of the H₂ molecule obtained in direct non-BO calculations by Wolniewicz, by our group, and by a perturbation-theory method derived in that work (see Table I in Ref. 4) based on the BO approach. While the agreement between our results and the results of Wolniewicz is very good, it is slightly worse when comparing with the results of Pachucki and Komasa. It seems that higher order corrections, which are more cumbersome to calculate, are needed to obtain a better agreement of the results obtained with the direct non-BO method and the perturbation theory results.

Let us consider an isolated atom or molecule formed by N particles with masses {M_i} and charges {Q_i}. We start with the positions of the particles first described using the Cartesian coordinates in the laboratory frame, $Ri$⁠. The laboratory coordinates and the corresponding linear momenta of the particles are

The laboratory frame nonrelativistic Hamiltonian of the system is

Next, the 3N-dimensional problem represented by the above Hamiltonian is separated into two problems. The first is a 3D problem of the motion of the center of mass of the system in space. The second problem is an $3N−3$-dimensional internal problem of the motion of the particles forming the atom or a molecule with respect to each other. The separation of the two problems is rigorous and is achieved by transforming Hamiltonian (2) to a new coordinate system. The first three coordinates, r₀, in this transformation are the coordinates of the center of mass in the laboratory coordinate frame and the remaining 3N − 3 coordinates are the internal Cartesian coordinates. These coordinates, denoted as r_i (i = 1, …, N − 1), are the position vectors of particles 2 to N with respect to particle 1. Particle 1 is called the reference particle. While any particle in the system can be chosen to be the reference particle, it is natural to assign this role to the heaviest one. The approach is a generalization of the textbook approach applied to solve the time-independent Schrödinger equation for the hydrogen atom. Usually in an atomic calculation, the reference particle is the atomic nucleus, and in a molecular calculation, it is the heaviest nucleus. If we denote n ≡ N − 1, then, in the new coordinate system, Hamiltonian (2) becomes

where q_i = Q_i+1 (i = 0, …, n) are the charges, $μi=m0mim0+mi$ are the reduced masses, M_tot is the total mass of the system, m₀ is the mass of the reference particle, and m_i = M_i+1. $∇ri$ is the gradient vector expressed in terms of the x_i, y_i, and z_i coordinates of vector r_i, r_ij = |r_i − r_j| = |R_i+1 − R_j+1|, and r_0i ≡ r_i = |r_i| = |R_i+1 − R₁|. The prime symbol is used to denote the vector/matrix transposition.

In the new coordinate system, the laboratory frame Hamiltonian (3) becomes a sum of the operators representing the kinetic energy of the center-of-mass motion, $Hnrcm(r0)$⁠, and the internal Hamiltonian, $Hnrint(r)$⁠,

where r is a 3n-component column vector. Its first three components are coordinates r₁, the next three are r₂, etc. In the present work, we are only concerned with the internal bound states of the system. Thus, only the eigenvalues and eigenfunctions of the internal Hamiltonian are calculated. The variational calculation in the framework of the Rayleigh–Ritz approach involves solving the generalized eigenvalue problem (finding eigenvalues ϵ and eigenvectors c),

with the (internal) Hamiltonian and overlap matrices H and S. The elements of the Hamiltonian and overlap matrices, $Hkl=⟨ϕk|Hnrint|ϕl⟩$ and S_kl = ⟨ϕ_k|ϕ_l⟩, are calculated in the chosen set of explicitly correlated basis functions, {ϕ_k}.

The computational time for calculating a matrix element of each of the two matrices is, in general, proportional to n³, as it involves matrix operations such as inversion and multiplication. The time is also proportional to the product of the factorials of the numbers of particles in the groups of identical particles (for example, in the case of the $H3+$ ion, this product of factorials is 3! × 2! as there are three identical nuclei and two identical electrons). This later dependency results from the number of permutation symmetry operators that need to be applied to the wave function to impose the proper permutational symmetry. When all matrix elements are computed, solving the generalized eigenvalue problem (5) also makes a sizable contribution to the overall cost of the calculation despite using a quick solver involving finding only a single eigenvalue and the corresponding eigenfunction. It should be noted that both the Hamiltonian and overlap matrices are complex Hermitian matrices, while the overlap matrix is also positive definite.

Upon examining the form of the internal Hamiltonian, one can see that it describes a system of n pseudoparticles with the masses equal to reduced masses μ_i and charges q_i (i = 1, …, n) moving in the central field of the charge of the reference particle, q₀. The pseudoparticles interact with each other via the Coulomb potential. In addition, their motions are coupled (correlated) through the mass-polarization terms, $−12∑i≠jn1m0∇ri′∇rj$⁠. As mentioned, (3) has the symmetry of an atomic Hamiltonian. Thus, the spatial basis functions used to expand the wave function of such a Hamiltonian need to be one-center functions that provide bases for irreducible representations of the SO(3) group of rotations.

As nuclei (nucleus in an atomic calculation) and electrons in the calculation are treated on equal footing, the basis functions need to properly represent this egalitarian approach. In this work, we use all-particle explicitly correlated Gaussian functions as basis functions for expanding the spatial component of the non-BO wave function. There are various forms of explicitly correlated one-center n-particle Gaussian functions (ECGs). The Gaussian exponent of these functions can be represented in the following form that shows its explicit dependence on the inter-particle distances, r_ij:

where r_i is the distance between particle i and the center of the Gaussian and α_ik and β_ij,k are the non-linear parameters of the Gaussian. The above function can be represented in the following alternative compact form:

where $A¯k$ is a 3n × 3n real symmetric matrix of exponential parameters. $A¯k$ can be written as $A¯k=Ak⊗I3$⁠, where I₃ is the 3 × 3 unit matrix, ⊗ denotes the Kronecker product, and A_k is a n × n symmetric matrix. To ensure square integrability of function ϕ_k(r), matrix A_k must be positive definite. This happens automatically if A_k is represented in the Cholesky factored form as A_k = L_k$Lk′$⁠, where L_k is an n × n, rank n, lower triangular matrix. ϕ_k(r) is automatically square-integrable for the L_k matrix elements being any real numbers.

As mentioned, the description of the n–n correlation and the nodes of the wave function of excited states requires including additional factors in the basis functions. It was shown that, in the case of a diatomic molecule, an effective factor can be a pre-exponential multiplier in the form of a non-negative power of the inter-nuclear distance.^3,5,6 For a molecule with more than two nuclei, the factor should be a product of non-negative powers of all intermolecular distances. Let us call the Gaussians with such multipliers “power Gaussians”. For a diatomic system, the power Gaussian has the following form:

where r₁ is the distance between the reference nuclei and particle i and 2p_k is its integer non-negative power (ranging from 0 to 250 in our calculations; we only use even powers as this considerably simplifies the calculation of the Hamiltonian integrals but has almost no effect on the accuracy of the calculation). For a triatomic molecule, the power Gaussian has the following form:⁷ 

where r₁, r₂, and r₁₂ are internuclear distances and $2p1,k$⁠, $2p2,k$⁠, and $2p12,k$ are their respective non-negative even-integer powers.

The larger is the 2p_k power value at fixed A_k, the more the two nuclei are separated from each other [i.e., the maximum of the Gaussian, which is located, for example, for Gaussian (8) on a sphere, shifts to a sphere with a larger radius]. Gaussians with the zero value of the powers need to be included to assure that the probability of finding two or more nuclei in a single point in space may not be exactly zero. In a calculation of an excited vibrational state, Gaussians (8) with different values of the 2p_k powers are combined to generate the radial oscillations and the nodes in the wave function. As the non-BO wave function for an excited ro-vibrational state of, for example, a diatomic system oscillates in terms of r₁, it also oscillates in terms of the distances between the electrons and the nuclei because the maxima of the electron densities appear in the same regions where there are maxima of the densities of the nuclei. In the non-BO wave function of a diatomic molecule, where the internal coordinate r₁ describes the relative position of the second nucleus with respect to the first (reference) nucleus, the internal coordinates r₂, r₃, etc., describe the relative positions of the electrons with respect to the reference nucleus. Thus, the oscillations of the wave functions in terms of r₁ are accompanied by similar oscillations of the wave function in terms of r₂, r₃, etc. However, as the wave functions do not have nodes in terms of these latter coordinates, there is less need to include in the basis functions pre-exponential multipliers being powers of the r₂, r₃, etc. distances. The dependency of the exponents of the Gaussians on these distances usually suffices in this case. This explains the effectiveness of the basis set of Gaussians (8) in the non-BO calculations of diatomic molecules.^3,5,6 However, this effectiveness can be increased by making the pre-exponential multipliers also dependent on powers of r₂, r₃, etc. This point is discussed further when the basis functions used in the present work are introduced.

In recent years, we have developed and tested an alternative form of ECGs for molecular non-BO calculations.^8–10 These alternative basis functions are all-particle explicitly correlated Gaussians with complex exponential parameters (complex ECGs or CECGs, for short). The general form of such functions for states with the rotation quantum number of zero (L = 0) is⁸ 

where $A¯k$⁠, as in (7), $B¯k$ is a real symmetric matrix of variational parameters, and $i=−1$⁠. $A¯k$ and $B¯k$ can be written as $A¯k=Ak⊗I3$ and $B¯k=Bk⊗I3$⁠. For L = 1 (M_L = 0) states, the appropriate form of the CECGs becomes¹⁰ 

where $zik$ is the z internal coordinate of particle i_k. Note that the value of the integer parameter i_k can be varied from 1 to n for each basis function and, thus, it has the corresponding index k. Basis functions (10) and (11) have been implemented.^8–10 Work on implementation of L = 2 (M_L = 0) CECGs,

is currently in progress. The limitation of L = 0, 1, or 2 in the non-BO atomic and molecular calculations means that only the states where the electronic excitations are limited to Σ, Π, and Δ states, states where the total rotational quantum number does not exceed three, and states of the mixed nature where the total angular momentum quantum number does not exceed three can be considered. Future developments may go beyond these constraints.

As mentioned, to ensure the square integrability of ϕ_k(r), matrix A_k in (7) is represented in the Cholesky factored form as A_k = L_k$Lk′$⁠. The L_k matrix elements, which are variational parameters in the calculation, can take any real values, positive or negative. In our recent works, we showed that functions (10) can be contracted to form pre-exponential multipliers being products of sin and cos functions dependent on squares of the inter-particle distances.¹¹ Such products multiplied by the real part of the exponent of the Gaussian, which “dumps” the sin/cos oscillations at longer distances, are capable of very effectively describing the e–e, n–n, and n–e correlation effects. They can also describe the radial oscillations and the radial nodes of the wave functions of excited states as effectively as the power Gaussians for both atomic and molecular systems. Finally, the CECG basis functions are universal and can be used in non-BO calculations of atoms and diatomic molecules, as well as in calculations of molecules with more than two nuclei.

The exact wave functions of states of atoms and molecules corresponding to the projection of the total orbital angular momentum vector on the z axis equal to zero are functions that are not complex. Even though complex Gaussians can be contracted into sin/cos linear combinations that are real functions, there is some value in using CECGs in the uncontracted form in the non-BO calculations. If the complex Gaussians are uncontracted and both A_k (or L_k) and B_k matrix elements are variationally optimized as independent parameters, it gives the basis more flexibility and may open some additional trajectories in the process of optimizing the nonlinear variational parameters, which in the end can lead to more compact expansion of the wave functions in terms of the basis functions. At the same time, the complex-space optimization approach may produce a wave function with a small imaginary component. This is acceptable as the function is variational.

The angular excitations of the system, for both the atomic and molecular cases, are described by including angular functions as pre-exponential Gaussian multipliers. As mentioned, for L = 1, the Gaussian should be multiplied by the $zik$ factor, where i_k = 1, …, n. The value of i_k in our calculations is treated as an additional variational parameter and is partially optimized. While in the atomic non-BO calculations, i is an electron index, in molecular calculations, it may correspond to either an electron label or a label of a nucleus. Thus, in a calculation of an excited rovibrational state, along with the main contributions to the non-BO wave function that describe angular excitations of the nuclei, there might be some minor contributions corresponding to angular excitations of the electrons. We showed that such contributions are indeed small, yet not negligible.¹² For rovibrational states where more than one particle is angularly excited (we calculated such states for atoms¹³), there can be contributions to the wave function from Gaussians where both nuclei (nucleus for an atom) and electrons are simultaneously angularly excited. The CECGs with appropriate mixed Cartesian angular factors can effectively describe such excitations.

The algorithms for performing non-BO calculations with CECGs are implemented in a computer program written in FORTRAN and employ MPI (Message Passing Interface) protocol for parallelization on distributed memory systems. The implementation includes algorithms for calculating the matrix elements of the internal-Hamiltonian matrix and the overlap matrix. At present, the calculations can be performed with basis functions (10) and (11). The implementation also includes the matrix elements of the analytical energy gradient determined with respect to the elements of L_k and B_k matrices, i.e., the Gaussian non-linear parameters. The procedure for calculating the gradient is described in Refs. 8 and 14. The implementation is general and can be applied to an arbitrary number and types of particles. As the calculations of matrix elements requires very little communication between different MPI processes, this part of the code is highly scalable.

The scheme implemented in the computer program is fairly straightforward. The major goal of the calculation is to generate a CECG basis set for expanding the wave function of the desired state of the system and to calculate its energy and other properties. The basis set generation is based on variational minimization of the total nonrelativistic internal non-BO energy of the considered state of the system. The calculation of the energy involves constructing the Hamiltonian and overlap matrices and solving the generalized eigenvalue problem. The variational energy minimization employs the energy gradient to navigate the search for an energy minimum (or at least a low enough point in the space of nonlinear variational parameters). The gradient is dependent on the derivatives of the Hamiltonian and overlap matrix elements determined with respect to the matrix elements of the L_k and B_k matrices.

The basis set is generated and tuned independently for each considered state of the system. The basis set is grown from a small set of Gaussians that can be obtained by a very thorough optimization, constructing it from a set of single-particle Gaussian orbitals, or taken from an already calculated neighboring state of the system. Using a small set of Gaussians whose parameters are obtained from a random number generator is also a possibility. The growing process involves adding basis functions one by one or in small groups of functions. The initial guess for a new added basis function is obtained by randomly perturbing the non-linear parameters of a subset of some function already included in the basis set and then by selecting the candidate that lowers the energy the most. For that function, the particle indices in the angular pre-exponential multiplier [for example, the i_k index in the $zik$ pre-exponential angular multiplier in basis function (11) or the i_k and j_k indices in function (12)] are optimized first. This is followed by an optimization of the function’s L_k and B_k parameters. At this point, the overlap integrals between the optimized function and the functions already included in the basis set are calculated to check for linear dependencies. If any linear dependencies exceeding a predefined threshold are detected, the function is discarded and a new function is selected and optimized. It is undesirable to have linearly dependent functions in the basis set, as this may lead to severe numerical instabilities and loss of precision in the calculation. After addition of a certain number of new functions to the basis set, the whole basis set is reoptimized by cycling over all functions one by one and reoptimizing their L_k and B_k exponential parameters (no reoptimization of the i index is done at this point). The growing of the basis set continues until the desired energy convergence is reached.

It should be noted that in our code, there is also an alternative option of reoptimization of all parameters of all basis functions simultaneously. In this case, in order to prevent the optimization process from yielding a set of highly linearly dependent functions, we add a penalty function as described in Ref. 15. In the approach, a positive-valued penalty term,

is added to the total variational energy of the considered state whenever two or more basis functions have the normalized overlap, S_ij, higher than the predefined threshold, t. After the addition, instead of minimizing the energy, the sum of the energy and the penalty term is minimized. β in (13) is a scaling factor. The penalty term depends on the overlap integrals between the basis functions and increases as these overlaps become closer to one (the basis functions are normalized). In the calculation, the penalty threshold is set to a value slightly smaller than one (e.g., t = 0.99). This ensures that during the optimization no two basis functions become almost linearly dependent.

The code, which contains a makefile that can be used to compile it, is placed in a GitHub repository and can be downloaded freely.¹⁶ A description of how to generate an input file is also supplied (some input parameters and variables used in the code are described in the comments of the Fortran source). We also provide examples of the input file (see inout.txt) and the corresponding output files. The examples concern calculations of the ground and the first excited states of the HD⁺ ion. The structure of the code reflects the purpose of the calculation, which achieves the lowest possible total nonrelativistic energy of the considered state of the system. The code contains subroutine DRMNG, modified by us, which realizes unconstrained optimization and is implemented via reverse communication.¹⁷ During the minimization process, DRMNG calls the procedure that calculates the energy and the energy gradient, which in turn calls the procedure to calculate the Hamiltonian and overlap matrix elements and their derivatives with respect to the non-linear parameters of the bra and ket basis functions. Thus, a significant amount of the computational time goes into calculating the matrix elements. The part of the code that does that is fully parallelized. To illustrate the efficiency of the parallelization, calculations are performed where the Hamiltonian and overlap matrices are generated for the basis set of 1500 CECGs for HD⁺ using different numbers of MPI processes (i.e., 28, 14, 7, 3, and 1). The timings for the calculations of matrices are shown in Table I. They were obtained on a server with Intel Xeon E5-2695v3 central processing unit (CPU) (14 cores/28 threads and 2.3 GHz base frequency). The data show that, until the number of MPI processes reaches 14, the computational time scales nearly linearly with that number. At 28 MPI processes, the scaling becomes less linear, which is partially due to the fact that hyper-threading in many modern CPUs is rarely capable of “doubling” the actual performance (gains due to enabling hyper-threading in most float-point intense applications are relatively small). In the table, we also include timings for calculations performed using 14 MPI processes using different numbers of basis functions. The numbers increase by a factor of two and the timings for calculating the matrix elements approximately increases, as expected, by a factor of four.

Finally, in order to illustrate how much time a typical calculation of a state of a molecular system takes, we performed a series of non-BO calculations for the ground state of the HD⁺ ion in which the CECG basis set is incrementally increased from 100 CECGs to 300 CECGs. Each time, the increment of the increase is equal to 50 CECGs. After each incremental increase of the basis set, 20 optimization cycles are performed for the basis functions, where in each cycle, all basis functions in the set are variationally reoptimized one by one using the gradient-aided optimization approach. The calculations were performed using seven MPI processes. The energies and the timings corresponding to each incremental increase of the basis set size are shown in Table II.

All three illustrative examples presented in this work concern molecular systems. For atoms, we used Gaussians (10), (11), and (12) without the imaginary components in the exponents (see our recent benchmark calculations of the four lowest ¹S states of the beryllium atom¹⁸ where we showed that the use of these functions can produce extremely accurate total energies for all Be isotopes). As real Gaussians, such as those used in the above-mentioned work on the ¹S spectrum of beryllium, are a partial case of the complex Gaussian basis set (10), the inclusion of the imaginary components in the Gaussian exponents would make these basis functions even better. However, for molecules, the effectiveness of the use of CECGs in ro-vibrational calculations has not been as well documented as for atoms. Thus, in the test calculations described below, we only consider molecules.

The first test concerns the complete vibrational spectrum of the HD⁺ ion corresponding to zero total orbital angular momentum. The spectrum consists of 23 bound states (⁠$v$ = 0, …, 22). The total energies of the ro-vibrational states obtained in the calculations performed with CECGs are shown in Table III. The results are compared with the energies obtained in the calculations performed with power Gaussians.¹⁹ The number of CECGs for each state is significantly smaller than the number of power Gaussians. For example, for the lowest state, the largest number of CECGs is 1500, but there were 4000 power Gaussians in the largest basis set of this kind. However, as one can see, already with 1500 CECGs, the total energy of −0.597 897 968 588 a.u. is slightly lower than the best power Gaussian result of −0.597 897 968 577 a.u. For the first excited state, the 1500 CECG result is slightly higher than the 5000 power-Gaussian result. This trend continues as the size of the power-Gaussian basis set increases for higher states to reach 7000 for the top $v$ = 22 state, while the number of CECGs increases to only 1700. The reason we are only able to show CECG results up to 1700 functions are due to a yet unresolved issue with developing a fast and stable (for highly excited states) eigenvalue solver that can update the solution of the complex generalized eigenvalue problem at the cost of $O(K2)$ arithmetic operations (K is the size of the basis). Therefore, at large K values that exceed ∼1000, solving the generalized eigenvalue problem becomes a bottleneck in the current implementation. However, the analysis of the energy convergence indicates that is it likely that the CECG energies for all 23 states will eventually become lower than the power-Gaussian results as the number of CECGs keeps increasing.

As the level of vibrational excitation increases and the wave function acquires more radial nodes, more basis functions are needed to describe the oscillations of the wave function. The primary vibration oscillations appear in the HD⁺ non-BO wave function in terms of the coordinate representing the internuclear distance (the length of vector r₁ in the present calculations). However, the vibration quantum number $v$ is only approximately a good quantum number, as the vibrational motion is coupled in HD⁺ with the motion of the electron. The coupling increases with the vibrational excitation. Also, as more oscillation appears in the wave function in terms of r₁ with vibrational excitation, oscillations also appear in terms of the r₂ coordinate, which represents the distance of the electron from the reference nucleus. As mentioned, this happens because the wave-function maxima in terms of r₁ coincide with the maxima in terms of r₂ as the electron follows the moving nucleus (pseudonucleus). In the power Gaussians, there is no term that can directly represent the r₂ oscillations and they have to be described “indirectly” through the exponential correlation terms. On the other hand, appropriate linear combination of CECGs generates real Gaussians with sin/cos preexponential multipliers dependent on $r22$⁠, which can directly describe the needed oscillation of the wave function. This explains the effectiveness of CECGs in molecular non-BO ro-vibrational calculations. This also suggests that CECGs should be equally effective in non-BO ro-vibrational calculations of molecules with more than two nuclei as they are in the calculations of diatomics, such as the HD⁺ molecular ion.

The HD⁺ system is a good example where the Born–Oppenheimer approximation leads to unsatisfactory results as it fails to describe the asymmetry of the charge distribution as the vibrational excitation increases. In BO nonrelativistic electronic calculations, the average proton–electron and deuteron–electron distances are the same regardless of the level of the vibrational excitation (if no adiabatic corrections to the wave functions are included). However, when the non-BO approach is used in the calculations, the densities of the particles forming the system show charge asymmetry that increases with the level of the vibrational excitation. The values of the average interparticle distances calculated using the non-BO wave functions generated with the code described in this work presented in Table IV show that, while for lower vibrational states, the asymmetry (as measured by the difference between the proton–electron and deuteron–electron distances) is small, for the last two states before dissociation, the electron is almost entirely localized at the deuteron. Simultaneously, the character of the bond in the ion changes from being basically covalent at lower vibrational states to ionic at the top states.

The second test concerns a triatomic molecule, the $D3+$ ion. It is chosen because, due to its nuclei being bosons, the permutational symmetry of the nuclei is easier to implement in the wave function than in the case of $H3+$ with fermionic nuclei (⁠$H3+$ non-BO calculations are now possible with the newly implemented algorithm¹⁰ and will be pursued in the near future). When Gaussians (10) with B_k set to zero (real Gaussians) are used in a molecular calculation, for example, for H₃⁺, the natural occurrence of linear dependencies between the basis functions prevents the calculation from converging. The reason for this behavior can be explained by the optimization of the basis functions trying to describe the value of the wave function becoming equal to zero as r₁ → 0, r₂ → 0, or r₁₂ → 0. In order to achieve that for real Gaussians [i.e., Gaussians (10) with B_k = 0], the optimization has to effectively generate pre-exponential factors in the Gaussians that are products of powers of r₁, r₂, and r₁₂. This can be achieved by “differentiating” the real Gaussians with respect to appropriate matrix elements of L_k. During the optimization, this differentiation effectively occurs by making two or more Gaussians almost identical and having them contribute to the wave function with a linear coefficient with equal absolute values, but opposite signs. This effectively generates a “numerical” derivative that approximates the analytical one. This behavior explains the tendency of the optimization to create the linear dependencies, which make the optimization ineffective and even numerically unstable. For example, we observe such behavior in the optimization of real Gaussians for vibrational states of the $H2+$ molecule. By adding the imaginary components to the real Gaussians [i.e., by making B_k ≠ 0 in Gaussians (10)], thus effectively allowing the sin/cos pre-exponential factors to appear in the Gaussians, the zeroing of the Gaussians at r₁ → 0, r₂ → 0, or r₁₂ → 0 is achieved. Then, there is no need for the linear dependencies to appear in the calculation. This feature is tested for $D3+$ in calculations of its ground and first-excited vibrational states corresponding to the total rotation quantum number of zero. In the test, we use 400 CECGs for each state. The basis set is separately optimized for each state using the gradient-based variational energy minimization. The all-parameter optimization approach is employed. The linear dependencies appear not to be an issue in the calculation. The energies of the lowest two states obtained in the calculations are −1.328 377 6 a.u. and −1.316 858 8 a.u. These two results are preliminary. To obtain good-quality results, the basis sets for the two states have to be significantly increased. Also, more ro-vibrational states need to be calculated to fully assess the utility of the computer code presented in this work.

The third example of non-BO calculations concerns the largest system we have ever calculated with the non-BO approach using all-particle explicitly correlated Gaussian functions. It is the eight-particle (two nuclei plus six electrons) boronmonohydride molecule, BH. The calculations are performed for the lowest energy state of the system. As the above-described HD⁺ calculations showed, CECGs can be equally effective in non-BO molecular calculations as explicitly correlated power Gaussians,¹⁹ which are products of CECG Gaussian exponents with $B¯k=0$ and even non-negative powers of the internuclear distance. Thus, in the BH calculations, which started long time before the present work was initiated (the calculations are still running), we use the power Gaussians. At present, the basis set size reached over 4000 functions. This required about two years of continuous calculations on a 32-core computer. As for the memory requirements, for relatively large basis sets, the bulk of allocated RAM is used for storing two square matrices of size K. When double precision is used in the calculation, it translates into ∼16 × K² bytes. Only a small amount of disk space is used. The majority of CPU time goes into calculating the matrix elements of the Hamiltonian and overlap matrices and those of the energy gradient determined with respect to the Gaussian exponential parameters. Because there are six identical particles (electrons) in the BH molecule, each matrix element requires the calculation of 6! = 720 terms. As each Gaussian in the BH calculations contains 28 independent parameters, for 4000 Gaussians, the total number of parameters that are optimized is 112 000.

The results of the BH calculations are shown in Table V. They include the energy values obtained at different stages of the basis-set growing process. The reported energy values correspond to 500, 1000, 2000, 3000, and 4000 Gaussians. As one can see, the energy convergence is good, but with 4000 basis functions in the basis set, the energy is converged to about the fifth significant figure only. The BH example provides a good illustration of the effort involved in non-BO molecular calculations of larger molecules.

In this work, we describe our new computer code ATOM-MOL-nonBO for performing non-Born–Oppenheimer calculations of atoms and molecules with an arbitrary number of electrons and for molecules with an arbitrary number of nuclei. The nuclei and electrons are treated on equal footing. The basis sets of all-particle explicitly correlated complex Gaussians are employed in the calculations. The basis set is universal and can be used to calculate ground and excited electronic states of atoms and ro-vibrational states of molecules. The code is tested in the calculations of the complete vibrational spectrum of the HD⁺ ion in the ground rotational state and the lowest two vibrational states of the $D3+$ ion. Some results are also shown from non-BO calculations of the ground state of the eight-particle BH molecule, which is the largest system we have ever calculated within the non-BO approach.

The code continues to be developed and new features are being added. Currently, it allows us to calculate states corresponding to the total rotation quantum number of zero and one. One of the most interesting applications of the code, which will be performed in the future, is calculations of the ro-vibrational spectra of $H3+$ and its isotopologues. To increase the accuracy of the calculations, we are planning to add algorithms for computing the leading relativistic corrections (the algorithm for calculating these corrections for the ground rotational state are already implemented²⁰).

See supplementary material for the source code used in the present calculations can be downloaded/cloned from a git repository.¹⁶ It contains the FORTRAN source (files with extensions .f90 and .f), a UNIX makefile, and several example input files, which set up calculations for the ground and first excited states of the HD⁺ ion with 500, 1000, and 1500 basis functions. Example output files are also provided. Users need to make sure that a working installation of the MPI library is available (e.g., Open MPI²¹) and adjust the makefile according to their local configuration.

This work was supported by a grant from the National Science Foundation (Grant No. 1856702). S.B. also acknowledges funding from MES RK state-targeted Program (No. BR05236454) and NU faculty development (Grant No. 090118FD5345).