Solving One-Electron Systems in a Novel Gaussian-Sinc Mixed Basis Set
Page 1
Solving One-Electron Systems in a Novel Gaussian-Sinc Mixed Basis Set
Jonathan L. Jerke, Young Lee and C. J. Tymczak
Department of Physics, Texas Southern University, Houston, TX 77004, USA
A novel Gaussian-
Sinc mixed basis set for the calculation of the one-electron electronic structure within a uniform magnetic field in three dimensions is presented. The one-electron system is used to demonstrate the utility of this new methodology and is a first step in laying the foundation for further development of many-electron atomic and molecular methodology. It is shown in this manuscript how to effectively calculate all basis set integrals, which includes the mixed Gaussian-
Sinc integrals, with a fast and accurate method. The
Sinc basis is invariant to the choice of the position of the Coulomb potential, as opposed to traditional grid based methods. This invariance guarantees that the choice of the grids origin has no effect on the electronic structure calculation. This is because the Coulomb potential is treated properly in this methodology, as opposed to DVR methodologies. The off-diagonal terms are sparse but very important around the Coulomb singularity. In general, five to six significant digits of accuracy on all converged results without the linear dependency problems of the Gaussian methodologies are achievable. This methodology is applied to calculate the ground state energy of H atom, H ion and H ion in magnetic fields up to a magnetic field strength of 2.35x10 G (10,000 au). From these calculations it is shown that H ion is unstable without relativistic considerations.
I. INTRODUCTION
A fundamental problem with electronic structure theory is the accuracy of the basis sets that are utilized [1]. Electronic structure basis set methodologies are usually focused on accurately representing either atomic valence states or atomic core states. Core focused methodologies can achieve unprecedented accuracy with the atomic core states with localized basis elements, e.g. Gaussian Type Orbitals (GTO) [2], but fail with valance states. In contrast, valence methodologies that use broad semi- to delocalized basis functions to capture the intricate inter-atom bonding orbitals will achieve high accuracy with bonding orbitals, e.g. Plane Waves [3], but usually fail to capture the highly localized core orbitals. Ideally one would like to be able to move between these two basis themes in a unified construction. In general, core bases sets require high angular symmetry basis functions in order to span the molecular arrangements of electrons found in molecular and condensed matter systems. However, these high angular symmetry basis functions have two difficulties: a) bonding orbitals can have a symmetry not reflected in the angular symmetries of the core focused bias set used leading to an ill-conditioned expansion; and b) calculation of angular symmetry type basis functions scale as O( L ) limiting their utilization because of computational complexity. However, it should be noted that these methods have been very successful in solving first row elements, e.g. Crystal03 [4], FreeON [5], or Gaussian [6] for example, but this can be attributed to the low angular symmetry of the constituent atoms. The electronic valence structure has been adequately modeled with plane-waves basis functions for the primary purpose of spanning the molecular bonding orbitals [3,7-9]. A plane-wave Page 2 density functional theory (DFT) like ONETEP has found many successes [10]. However, we have chosen to not proceed with a DFT approximation. Furthermore, one can attain the accuracy of plane waves but with A semi-local basis set. A function does exist that is universal enough to reconstruct both plane-waves and be local enough to define properly all interaction tensor terms. The Sinc function is a semi-local function that is both universal and local. Information theory sets a foundation for the
Sinc function. We will show in this paper that this replacement to plane-wave is local enough to define all interaction tensor elements perfectly, but with the accuracy and utility of the plane wave methods The valence structure can be spanned with a regular lattice of
Sinc functions. The
Sinc is a name given by Shannon [11] and coworkers for the cardinal sine defined in antiquity. Information theory states that band-limited functions can be perfectly represented using a
Sinc basis. Therefore, compact valence structure can be determined with a finite basis. There is a long history of the
Sinc function in computational science. ONETEP also uses a periodic version of the
Sinc function they call pSincs.
However, we do not use a Fast Fourier Transform (FFT) to define the kinetic energy, instead we use the matrix representation of the
Sinc derivative to directly represent the kinetic energy of the wave function without loss of generality. The kinetic energy used is found in the literature as the Discrete Variable Representation (DVR) methodology [9,12-14]. However, DVR fails to achieve a solution for the singular Coulomb potential. One property that the DVR method relies upon to obtain accuracy is that all off-diagonal matrix elements approach zero (in relation to the on-diagonal) as the lattice spacing approaches zero. For most cases this is true with the exception of the singular Coulomb potential, which is the source of the previously mentioned problem. Our new method overcomes this problem and in contrast to the DVR methodologies can solve the Coulomb potential with a variational basis set approach in full three-dimensional space. For a review of the
Sinc representation as used in computational science see Cattoni [15]. For additional accuracy and to correctly represent the atomic core region, this methodology can mix GTO and
Sinc bases together to achieve a highly accurate representation. This highly flexible mixed basis set allows the choice to either focus on valance, core, or both at the same time in a unified framework. A new methodology of computing the integrals of the
Sinc and Gaussian functions with the Coulomb potential is developed, which is described in detail in Section III. In this methodology, we are able to compute all integrals within numerical precision necessary to achieve the variational bound, including all the Gaussian-
Sinc mixed integrals. In what follows we demonstrate the utility of the
Sinc function as a variational basis of electronic structure. The
Sinc function is semi-local enough to define all interaction and matrix elements. Coupled with a Gaussian basis near the cores of atoms, we can achieve a near perfect construction of valence and core orbitals. We generally achieved five to six significant digits of accuracy on all converged results without the linear dependency problems of the Gaussian methodologies [4-6], which will eventually lead to more powerful linear scaling techniques. This methodology was then applied to calculate the ground state energy of H, H ion and H ion in magnetic fields completely up to a magnetic field strength of 2.35x10 G (10,000 au). From an extrapolation from low field calculations there is little evidence for stability of the H ion without relativistic considerations, which is in direct opposition to the conclusions of Turbiner [16]. This paper is organized as follows. In Section II a review and present what is known about the
Sinc function analytic properties. In Section III, a novel method for constructing the Coulomb
Page 3 integrals necessary to compute efficiently all matrix elements is described. In Section IV exhibits basis definitions and consistencies of this scheme. Section V discusses results of one-electron systems including a pure Gaussian calculation, a pure
Sinc calculation, and a Gaussian-
Sinc mixed basis. And finally the analysis of the H atom, H and the H in magnetic fields up to 10,000 a.u. is presented. II. THEORY: THE
SINC
BASIS
For a complete treatise of the subject of the
Sinc function see Stenger [17]. An outline follows with useful and important properties of the
Sinc function for general considerations and utility.
IIA. Discrete Representation
The
Sinc function that we consider is normalized and defined as:
Sinc ϑ ( ) = d Sin πϑ ( ) πϑ , (2.1) where ϑ is the dimensionless parameter that has a place of an index or sub index written in the form of x/d , and d is the distance between lattice sites in one-dimension. The Sinc function also is zero at all lattice sites except its own, in other words dSinc n ( ) = δ n . The Sinc function is orthonormal on a regular lattice, δ nm = Sinc x d − n ( ) −∞∞ ∫ Sinc x d − m ( ) dx . (2.2) The origin of the lattice is chosen to be zero here; transformations between different origins will be discussed in Section II.D. These methods offer a discrete representation where the sampling of the band-limited function is the content of the representation, f n = F x ( )
Sinc x d − n ( ) dx −∞∞ ∫ f x ( ) = f n Sinc x d − n ( ) n = −∞∞ ∑ f n = d dk π % F ( k ) e iknd −π / d π / d ∫ . (2.3) The normal convention for the Fourier transform is used. f n is equal to the value at nd of the convolution of the function, F(x), with a band-pass at a minimum wavelength of . A band-limited function is perfectly represented on an appropriate lattice, (cid:0) f x ( ) = F x ( ) , whereby the sampling of the function equals the coefficients up to an overall norm, f n = d F ( nd ) . (2.4) The Sinc function behaves as a discrete delta function on the projection of the function into the band-pass sampled by Shannon’s sampling theorem. Engineers call this a low-pass filter of a function [18]. For utility we name these functions Eikons of the functions, because they are similar to the images of the function being represented. Thus, the coefficients on the lattice are
Page 4 the sampling of the low-pass filter of the function, which is seen clearly in Eq. (2.3). The continuous inner product of Eikons is the same as their sampling. The continuous inner product is
F G = F * ( x ) G ( x ) dx −∞∞ ∫ . (2.5) Where Eikons, F and G , are specified by f and g , F G = dx f n * Sinc ( x d − n ) g m Sinc ( x d − m ) m ∑ n ∑ −∞∞ ∫ = f n * g nn ∑ . (2.6) The convolutions of the Sinc functions are always translations of band-limited functions; thereby an integer translation of
Sinc is a Kronecker delta function on the lattice. These expressions are exact for band-limited functions.
IIB. Scaling Relationship
The
Sinc function is a universal scaling function. The scaling relationship is due to their band pass nature. A wider band pass admits a narrower band-pass; therefore one can express a
Sinc at one scale with a lattice of
Sincs at a finer scale:
Sinc xd − n ( ) = d ' Sinc md ' d − n ( ) m ∑ Sinc x d ' − m ( ) d ' < d , (2.7) The Sinc function is completely expressible in terms of other
Sinc functions at any sub-scale. This property can be very useful in making some calculations tractable. The
Sinc function is a unique universal function. Only band-limited functions have a scaling relationship that allows them to be expressed in terms of themselves at different origins and scales. All band-limited functions can be expressed as
Sinc functions. Therefore, the
Sinc function is the unique universal scaling function of band-limited functions.
IIC. Incompleteness
The
Sinc function is a basis that spans band-limited functions. The inner product with a Sinc is necessarily a projection of the state into the band-limited space; necessarily the
Sinc basis is not complete. Instead, the semi-completeness improves monotonically with the quality of the lattice,
Sinc x − x ' d ( ) = d Sinc x d − n ( ) Sinc x ' d − n ( ) n ∑ . (2.8) This should be contrasted with the standard closure relation for orthogonal basis sets where the left hand side of the equation is a Dirac-delta function. This point has been treated incorrectly in the past. Light [13] claims the Sinc basis is complete. This cannot be true because the
Sinc
Page 5 function only forms a complete basis up to the band edge [17]. This property is the root approximation used by the DVR community [9].
IID. Lattice Derivative:
The lattice derivative is defined by infinitesimal translations of the Eikon. This procedure leads to semi-local derivative matrix, which when multiplied by the Eikon gives the Eikon of the derivative. One could reconstruct the Eikon continuously into space by using the
Sinc as the kernel of the interpolation vector. Translations are sampling of the band-limited function at a different origin, and because of Shannon’s theorem each sampling is equally good. Referring to Eq. (2.3), the Eikon is translated by, a , by the operator, (cid:0) T nm a ( ) = Sinc n − m + a d ( ) , (2.9) where T ( a ) operating on the Eikon f gives, T a ( ) f x ( ) = T nm a ( ) f m Sinc x d − n ( ) m ∑ = f x − a ( ) . (2.10) This translation is equivalent to sampling the interpolation of the Eikon at a different origin. We express infinitesimal translations to derive the lattice derivative: D (1) = lim ε→ T ε ( ) − T −ε ( ) ε . (2.11) The limit can be taken without respect to any particular Eikon. The derivative matrix is then, (cid:0) D nm (1) = d n = m − ( ) n − m n − m , n ≠ m ⎧ ⎨ ⎪ ⎩ ⎪ . (2.12) The expression for D is semi-local, since an isolated single nonzero lattice element is a Sinc function, which is semi-local but zero on the lattice sites. The derivative is necessarily related to the translation operator via: (cid:0) d T x ( ) dx = D (1) T x ( ) . (2.13) This expression shows that the Sinc basis can span the derivative of all of its constitute members. We agree with DVR methods [9] that using the exact expression for the second derivative is highly accurate and preferred for the purposes of calculating the derivatives in electronic structure. This expression for the second derivative can be computed via the multiplication of two derivative matrices summed to infinity. Therefore, the second derivative matrix is given,
Page 6 (cid:0) D nm (2) = d − π n = m − − ( ) n − m n − m ( ) , n ≠ m ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ . (2.14) Additionally, all these derivative matrices obey the derivative chain rule exactly, as they should (cid:0) D ( n + = D (1) D ( n ) . (2.15) IIE. Commutation Relations
The commutation relationship of the derivative and the position operator is correct on the lattice up to a factor describing the Fourier component on the edge of Shannon’s sampling criterion. This represents the failure of the position operator in a band-limited space more than a true failure of Calculus, since the derivative operator maps the state into the same band-limited space, but the Fourier transform of the position operator is not compact in Fourier space. From Calculus the continuous derivative satisfies the product rule. Specifically on a continuous domain, ∂∂ x , x [ ] F x ( ) = F x ( ) . (2.16) Let us now consider the lattice derivative: The sampling of the derivative is equal to the derivative matrix acting on the sampled function. We propose to check this correspondence by applying the commutator in Eq. (2.16). We derive the commutation relation of the derivative and position operator analytically where D (1) is the derivative matrix and X is: X nm = d n δ nm . (2.17) From this the commutation relation follows, D (1) , X ⎡⎣ ⎤⎦ ij = D (1) ik X kjk ∑ − X ik D (1) kj = − ( ) i − k i − k k δ kj − i δ ik − ( ) k − j k − j k ∑ = − ( ) i − j i − j j − i − ( ) i − j i − j = i = j − ( − i − j i ≠ j ⎧⎨⎪⎩⎪ . (2.18) Eq. (2.18) is equal to the following form, [ D , X ] = I − A , (2.19) Where the A matrix is given by, Page 7 A ij = ( − i − j . (2.20) We refer to A as the alternating matrix. The null space of A is the number of basis elements minus one. The alternating vector is its span and has a wavelength of Shannon’s sampling theorem, specifically . The alternating matrix stands as a reminder of the sampling theorem, and can only be suppressed with proper encapsulation of the information on the lattice. Alternatively, one can notice that the alternating vector is not translational invariant, so it is an ambiguous object that touches on all possible representations of the same alternating state under translation of the basis. III. COULOMB INTERACTION AND MATRIX INTEGRALS
An integral equation is presented to calculate Coulomb terms between Gaussian and Gaussians,
Sincs and
Sincs , and Gaussians and
Sincs . The equations are modular and built on units of the basis in one dimension since these functions are separable in three dimensions. The final expression of the matrix and interaction elements is tractable with available computer technologies. The interaction term is defined by the following wave functions, (cid:0) H ( r r ) and (cid:0) G ( r r ) , d v r d ′ v r G ( v r ) H ( ′ v r ) v r − ′ v r ∫∫ . (3.1) The Coulomb potential can be expressed in terms of a Gaussian Fourier integral, = d v r d ′ v r G ( v r ) H ( ′ v r ) d v k (2 π ) d α ∞ ∫ e − v k /(4 α ) − i v k ⋅ ( v r − ′ v r ) α ∫∫∫ = π d α ∞ ∫ d v k (2 π ) e − v k /(4 α ) α % G ( v k ) % H ( − v k ) −∞∞ ∫ ⎡⎣⎢ ⎤⎦⎥ . (3.2) The Fourier transform of a function is defined in the normal convention, % F ( v k ) = d v rF ( v r ) e − i v k ⋅ v r ∫ F ( v r ) = d v k (2 π ) % F ( v k ) e i v k ⋅ v r ∫ . (3.3) Beginning with a six-dimensional integral, one can immediately reduced this to a four-dimensional integral. At all points this expression of the Coulomb integral is exact. We restrict our final calculation to only separable basis elements in three dimensions, F ( v r ) = F x ( x ) F y ( y ) F z ( z ) . (3.4) Eq. (3.2) separates inside the alpha integral into three one-dimensional integrals, Page 8 I i ( α ) = dk (2 π ) e − k /(4 α ) α % G i ( k ) % H i ( − k ) −∞∞ ∫ ⎡⎣⎢ ⎤⎦⎥ . (3.5) The next step is a reduction of this to a one-dimensional “alpha integral” that is the product of one-dimensional Fourier integrals, dr d ′ r G ( r ) H ( ′ r ) r − ′ r ∫∫ = π d α I x ( ∞ ∫ α ) I y ( α ) I z ( α ) . (3.6) Obviously the original interaction element expression is effectively intractable to numerical integration. However, Eq. (3.6) for separable basis elements is numerically tractable and very efficient. It takes about a second on a single CPU for any term to machine precision. Further improvements in speed can be afforded by analytically expressing Eq. (3.5). The Gaussian- Sinc terms are separated inside the “alpha integral”, leaving one-dimensional integrals to be computed that can be of the
Sinc - Sinc , Gaussian-
Sinc , or Gaussian-Gaussian variety. The
Sinc-Sinc term has an alternative but equal useful expression,
Sinc ( x d − n ) Sinc ( x d − m ) = d dk π dk π e − i k ( x − dn ) + k ( x − dm ) ( ) −π / d π / d ∫ −π / d π / d ∫ . (3.7) The Gaussian- Sinc term is, βπ ( ) e −β ( x − x ) Sinc ( x d − n ) = βπ ( ) e −β ( x − x ) d dk π e − ik ( x − dn ) −π / d π / d ∫ . (3.8) The Gaussian-Gaussian term is also Gaussian, β π ( ) β π ( ) e −β ( x − x ) e −β ( x − x ) . (3.9) This completes the list of all necessary basis constructions needed to construct all interaction terms in the Gaussian– Sinc basis. All corresponding Fourier transforms of the functions listed above can be analytically expressed. The Fourier transform of the
Sinc-Sinc term is analytically expressible. A single
Sinc function has a Fourier Transform of a top hat and a double Sinc function has a Fourier transform of a triangle:
Page 9 − ie iknd − e i ( m − n )( kd −π ) + e i ( m − n ) π ( ) π ( m − n ) , 0 < k < π / d , m ≠ nie iknd − e i ( m − n )( kd + π ) + e − i ( m − n ) π ( ) π ( m − n ) , − π / d < k < m ≠ ne iknd π − kd ( ) π , − π / d < k < π / d , m = n ⎧⎨⎪⎪⎪⎪⎩⎪⎪⎪⎪ . (3.10) The expression for the Gaussian- Sinc
Fourier transform is more difficult and presents a larger problem generally. This expression is written in terms of Error functions, d β π ( ) e iknd e −β ( nd − x ) Erf − k + π d − i β ( nd − x )2 β ⎡⎣ ⎤⎦ − Erf − k − π d − i β ( nd − x )2 β ⎡⎣ ⎤⎦ ( ) . (3.11) The ultimate problem with Error functions is that they are not analytically integrable. This presents a problem for further analytic reduction. Analytic expressions are desirable because of the speed and accuracy for which they can be calculated. The analytic expression for the Fourier Transform of a Gaussian-Gaussian term is well known to be a Gaussian. The position of the Gaussian is a beta weighted average and the distance between them has Gaussian suppression, β β β + β ( ) ( ) e − k β + β ( ) −β β ( x − x ) β + β − ik β x + β x β + β . (3.12) This is a trivial point in this context, but may be of more general interest to the Gaussian community. This formula allows us to integrate any Gaussian integral. In fact it is straightforward to use Eq. (3.5) to solve for any angular symmetry type of a Gaussian-Gaussian interaction term generally under any condition. One could also imagine different betas in every direction. At this point we can see why plane-wave basis elements are intractable. The Fourier transform of a plane wave is a Dirac delta, and equation Eq. (3.5) prescribes two Dirac Deltas per integral per dimension. Therefore the result is undefined. An adaptive integral procedure [19] can integrate these integrals at a rate of one element per second. However, this can be substantially speed up for all one-body matrix elements allowing for the very rapid calculation times. We have calculated all matrix elements in this paper with numerical precision of 9 significant digits at a rate of about 1000 per second. IV. THE SCHRÖDINGER WAVE EQUATION IVA. Gaussian-
Sinc
Basis
The Gaussian basis extends and compliments the
Sinc basis allowing for a more accurate treatment of the “cusp” of the wave function without the linear dependences of broad Gaussians. The Gaussian basis elements are placed at the atomic centers to accurately span the atomic core structure. The Gaussian exponents are restricted such that inter-atomic overlap is below a defined threshold, significantly simplifying computation because these Gaussian-Gaussian terms from one atom to another can usually be neglected.
Page 10
IVB. Upper Bounded
Without magnetic fields, we are using the true matrix elements of the basis functions mentioned and therefore remain variationally bound. One can prune the Sinc basis without loss of this principle. We do not need to assume periodicity or other boundary conditions. The construction guarantees that the wave function goes to zero at infinity. The
Sinc function is the proper way to sample the wave function in a digital representation. The
Sinc function is the interpolation kernel for a sampling of a function. The tricky point is the lattice scale, however it is conceivable to integrate even multi-resolution
Sinc functions. Wavelets and multi-resolution run parallel in the
Sinc formalism [ref].
IVC. Definitions Momentum Operator : The momentum operator P is proportional to the derivative operator; therefore we can define it as, P = − i h D (1) . (4.1) Canonical Commutation Relations : Using Eq. (2.19), we evaluate the canonical commutations relationships used in Quantum Mechanics as, [ X , P ] = i h ( I − A ) . (4.2) Schrödinger Wave Equation in a Magnetic Field:
From Eq. (2.14) and Eq. (3.1) we construct the non-relativistic Schrodinger equation in a zero magnetic field for quantum mechanical problems in our Eikon method as, H = P m + V ⎛⎝⎜ ⎞⎠⎟ . (4.3) To include magnetic fields, we begin with the standard construction, P → P + ec A , (4.4) Which give us the Hamiltonian (cid:0) H = m P + ec A ( ) + V . (4.5) For a uniform magnetic field, where (cid:0) A = r × B , we can rewrite as for a magnetic field in the z-direction, H = − h m D (2) − ie h mc ( xD y (1) − yD x (1) ) B + e mc x + y ( ) B + V . (4.6) The mixed basis requires an overlap matrix to solve the generalized Eigen-value problem, which has been solved for real and complex Hamiltonians using the LAPACK computational eigen-solver package [20]. Page 11
V. ONE-ELECTRON SYSTEMS VA. Hydrogen Benchmark Calculation with Pure Gaussians
We have applied our new integral technology to compute the non-interacting eigen-states of Hydrogen using a pure Gaussian basis in order to benchmark the methodology. Previous Gaussian algorithms require recursive and/or many function evaluations to compute matrix element and scale poorly with angular symmetry,
O(L ) computational complexity per matrix element (where there are O(L ) matrix elements), with angular symmetry types [21]. Our integral methodology allow us to solve for any matrix element directly and computationally very efficiently, with O (L) computational complexity per matrix element. Calculating the non-interacting Eigen-state of Hydrogen demonstrates the quality of a solution of core states in this computational scheme. 100 GTOs are used with estimated error on each matrix element below 10 -11 . We used the simple formula for the Gaussian exponents, (cid:0) β i = α β i − ( ) (5.1) Where (cid:0) α = and (cid:0) β = with 25 s-type Gaussians and 3x25 p-type Gaussians, a total of 100 GTO primitives. The 1s Eigen-state has an energy of -0.49999999948, the 2s Eigen-state energy was -0.249999999967 and the 2p Eigen-state energy was -0.249999999942. The total calculation took approximately ten seconds of CPU time. VB. Calculation of Hydrogen Ground State with
Sincs
Next, the new
Sinc methodology is applied to the ground state eigen-energies of Hydrogen. We calculated the Hydrogen system with a radius of 7.5 au from the nucleus in spherically pruned basis of lattice elements with lattice spacing of d=0.5, 1.0, and 2.0 au. The results of this calculation are presented in Table I. Resultant errors found are 0.6% error at d=0.5, 4% error at d=1.0, and 15% error at d=2.0. For reference, we compare this to previous wavelet calculations [22-29]. Wavelets could be considered the state-of-art when it comes to semi-analytical grid basis. As a specific example, consider the results of Tymczak et. al. [30]. Tymczak reports 2% errors on the computation of the same system. We see that Tymczak’s sparsely filled 128 =2,097,152 grid results were comparable to the one grid, yet the one grid is significantly smaller. Additionally, it was stated in this work that there was a sensitivity of the calculated eigen-energies to the position of the nuclear center in respect to the grid, which would cause significant issues with the calculation of atomic forces using these methods. This problem is rectified with the Sinc basis. In Figure 1 we plot a grid of d=0.5 calculation of the probability density centered and translated. There is convergence of the energies between the centered and translation of the one-grid calculations of increasing range of the wave-function. From this we conclude that any variations of the energies are strictly due to edge effects and not any sensitivity of the solution to the position of the singularity. We find it difficult to compare our results to most other calculation because our results do not assume any symmetry. Almost every paper solving hydrogen has assumed separation of variables in spherical or cylindrical coordinates, which is irrelevant to almost all molecular systems. For instance, finite element methods [31] and DVR [32] solved Hydrogen to very high accuracy in that symmetry. Such papers are useful as benchmarks for this work, but cannot be compared equally to a calculation conducted in three dimensions. Three-dimensional calculations are much more flexible to molecular configurations and are therefore a much more useful tool in solving the electronic structure without bias.
Page 12
Additionally, in Figure 2 we plot the error in ground state energy of hydrogen as the grid spacing decreases, as can be seen we obtained exponential convergence with decreasing grid spacing. For comparison, we also show in this figure the convergence of pure Hydrogen in a magnetic field and the H ion. VC. Hydrogen in the Gaussian-
Sinc
Basis
In this section we extend our calculation of the ground state of Hydrogen using a mixed Gaussian-
Sinc basis set. To minimize the linear dependence but still achieve a span of the cusp with the Gaussian basis we have restricted the minimum exponent. Using Eq. (5.1), we varied the alpha and beta values with a set number of Gaussians until the energy was minimized. Our results are reported in Table III for the Hydrogen atom using the
Sinc -Gaussian basis set. As is expected, as we decrease the grid spacing the minimum exponents increases. One problem that this addresses is linear dependence of the basis set. The major contribution to this linear dependence is the Gaussian
Sinc overlap. This dependence becomes significant when the Gaussian exponent and the grid spacing become comparable. However, with careful optimization we can easilly avoid this issue, as is shown in the table.
VD. Results for Hydrogen in a Magnetic Field
In Figures 3 and 4 we plot the cross sections of the ground state wave functions at various magnetic fields strengths. The convergence of the solutions with increasing sampling rate and range of the neighborhood around the atoms are given in Table IV and Table V for 1 au and 10 au magnetic fields respectively. Figure 3 shows the cross sections along all three axes of the wave function at magnetic field of 1 au. The energy of this calculation has been known to be -0.3312 Hartrees [33] and we found -0.33114 Hartrees. At 10 au of magnetic field we plot the best solution in Figure 4. The known calculated energy from reference [33] is 3.2522 Hartrees and our result was 3.25222 Hartrees. Energies remained within tolerance under translation by 0.5 au transverse to the magnetic field.
VE. Results for the H ion in the Gaussian-Sinc basis We calculated the ground state configuration H in a Gaussian- Sinc basis, as shown in Table VI. From the table we see that we agree with the known energies to an accuracy of six significant digits, where the known energy of the system is reported by Scott et. al. [34] to be -0.60263 Hartree at 1.998 au separation. In this reference the author details the technique of searching a hyper-dimensional trial wave function parameter space to find the minimum energy. The advantage of our method is that we do not need to conduct such a complex multidimensional optimization.
VF. Result of the H ion in a Magnetic Field We report on the calculation of the H ion system in a magnetic field of 1 and 10 au. Results shown in Table VII and Table VIII are directly compared to established work on one electrons systems by Turbiner et. al. [16]. Specifically, Our calculation at 45 is of lower energy then previous calculation and is consistent under translation of 0.5 au, which is a validation of the utility of our method. VG. Results for the Trigonal H ion in a Magnetic Field
In Figure 5, the ground state wave function of trigonal H ion for a zero magnetic field and a magnetic field of 1 au. Freezing the positions of the atoms, the electron completely leaves the
Page 13 orbit of the third proton with this magnetic field. This is indicative of the instability of the trigonal H ion in a magnetic field. We also find that the linear H is always significantly more stable then the trigonal H , therefore we will only consider the linear configuration in the next section for which we will explore these systems in more detail.
VH. Energetic’s in Very High Magnetic Fields
We have conducted a large-scale scan from zero to 10,000 au of magnetic field for the H, H and H systems. From this we can conclude H system is unstable to decay into H without relativistic considerations. The results are reported in Table IX. In Figure 6 the calculated energies of the one electron species with respect to the energy of Hydrogen within the same field are plotted. Fitting the data to a functional form, (cid:0) F ( B ) = a + bB + cB d + eB , (5.2) which is also shown in Figure 6. We find that the energies converge to parallel lines with energy difference of approximately 0.2 Hartree for all field predictions up to the relativistic limit. All our calculations agree with reference [16] result to a hundredth of a Hartree for all magnetic fields considered. We disagree with the conclusion of reference [16] that the H ion becomes the more stable species at high magnetic fields and believe that the incorrect conclusion in this reference is do to the scarcity of data points in the calculation. We warn that calculations above 5,000 au are questionable at best and its better to fit trends below that field value. The difficulty of trans-10,000 au field calculations places very little weight on their final numbers of a bound state to H . Whereas a trend in a more trustworthy domain of field with a reliable extrapolation scheme has more weight. No aspect of our analysis biased the result to lead to parallel lines above 10,000 au of field. Even so, there never is a transition to an H ground state. We therefore conclude that H is unstable to decay in this non-relativistic theory. VI. CONCLUSION
We have demonstrated the use of a novel Gaussian-Sinc basis set for the calculation of the electronic structure within a uniform magnetic field in three dimensions. With this method, one can achieve four to six significant digits in the ground state energy for any proton configuration. One can extend this technology to nonzero magnetic fields to show that H is unstable in all magnetic fields without relativistic corrections. Additionally, we are in the process of coding this method to include a many-body description (Hartree-Fock level) that will allow us to study atoms and molecules in ultra-high magnetic fields.
VII. ACKNOWLEDGEMENTS
The authors would like to acknowledge the support given by Welch Foundation (Grant J-1675), the ARO (Grant W911Nf-13-1-0162), the Texas Southern University High Performance Computing Center (http:/hpcc.tsu.edu/; Grant PHY-1126251) and NSF-CREST CRCN project (Grant HRD-1137732). Dr. Jerke would also like to thank University of Houston for hosting him as a Visiting Assistant Professor while beginning to consider information theory to be the foundation of Quantum Mechanics. We would also like to thank G. Reiter for many fruitful discussions.
Page 14
VIII. FIGURES
Fig. . A pure
Sinc half-grid calculation of Hydrogen with the analytic solution superimposed in yellow. Translation of the potential successfully leads to translation of the solution’s wave function.
Page 15
Fig. 2.
The convergence of the ground state energy of systems with and without magnetic fields. a) Hydrogen at zero field; the answer is converging towards the analytically known energy of -0.500000. b) Hydrogen in 1 au magnetic field: the true energy is not known well enough to establish a linear convergence plot, but there is rough convergence to -0.3312 Hartree, Thirumalai et. al. [33]. c) we also include H in zero-field. The convergence approaches the true energy of - Page 16
Fig. 3. The cross sections of the ground state Hydrogen wave functions in a B=1au magnetic field using a Gaussian-
Sinc basis. The zero magnetic field solution is drawn in green. The x-y verses z components are slightly compressed. This calculation was generated with d =0.533 au and a box size of 7.5 au. For reference 1 au = 2.35052 x 10 Tesla.
Page 17
Fig. 4.
The cross sections of the ground state Hydrogen wave function in a B=10 au magnetic field using a Gaussian-
Sinc basis. The zero magnetic field solution is drawn in green. The x-y verses z components are highly compressed. This figure was generated with d=0.300 and a box size of 4.21 au. For reference 1 au = 2.35052 x 10 Tesla
Page 18
Fig. 5. a) Trigonal electron wave function of H using a Gaussian-
Sinc
Basis. The three protons are distributed in an equilateral triangle with length 5.25 au. The energy is -0.34 au. The lattice spacing is 1 au. b) Trigonal electron wave function of H in B=1 au magnetic field. This is the same proton configuration as is shown in figure 6a. The lattice spacing is halved, but box size is unchanged. The third proton receives almost no electron density. The size of the electron wave function decreases with increasing magnetic field. For reference 1 au = 2.35052 x 10 Tesla
Page 19
Fig. 6. The energies of H and H with respect to the energy of Hydrogen in the presence of a magnetic field using the Gaussian-Sinc basis set.. Page 20
IX. TABLES TABLE I:
Hydrogen in the
Sinc basis for decreasing grid spacing with a radius of 7.5 au. The number of basis elements on the three grids of decreasing quality are 14,000, 1700, and 250 respectively.
N d (au) T (Hartree) V (Hartree) E (Hartree) 250 2.000 0.284773 -0.702619 -0.417846 1700
TABLE II:
Results under translations with increasing
Sinc basis set range. All calculations are conducted with a grid of d = 1 au. We find the translational quality is dependent only on the edge effects that reduce with an increasing range of the
Sinc basis, as expected. Range (au) Energy (Hartree) Energy ( ½ d shift) 3.0 -0.397029 -0.416588 5.0 -0.475833 -0.475425 7.0 -0.479425 -0.479045 9.0 -0.479587 -0.479263 11.0 -0.479595 -0.479324 13.0 -0.479597 -0.479364 TABLE III:
The Hydrogen Ground State using a mixed Gaussian-
Sinc basis. We calculate the core part of the wave function using 3, 6 and 9 s-type Gaussians where we use the equation (5.1) and optimize the variables. The deviation against the correct energies is shown to be exponentially convergent in Figure 3. No. of Gaussian Grid Spacing d (au) α (cid:0) β Condition Number Energy 3 2.0 0.53 4.8 105 -0.498528 6 2.0 0.30 3.2 1905 -0.499967 9 2.0 0.25 2.7 8700 -0.499998 3 1.0 1.90 4.8 146 -0.499831 6 1.0 1.15 3.2 2386 -0.499995 9 1.0 1.05 2.7 6995 -0.499999 3 0.5 7.1 4.8 180 -0.499978 6 0.5 5.6 3.3 960 -0.499997 9 0.5 5.0 3.2 910 -0.499999
Page 21
TABLE IV:
The convergence of Hydrogen in a magnetic field of B=1 au. For reference 1 au = 2.35052 x 10 Tesla d (au) Range (au) E (Hartree) 1.00 4.00 -0.32638 0.89 4.50 -0.32965 0.80 5.00 -0.33073 0.73 5.50 -0.33105 0.67 6.00 -0.33113 0.62 6.50 -0.33116 0.57 7.00 -0.33116 0.53 7.50 -0.33116
TABLE V:
Convergence of Hydrogen in a B=10 au magnetic field. For reference 1 au = 2.35052 x 10 Tesla. d (au) Range (au) E(Hartree) 0.56 2.25 3.18700 0.50 2.53 3.24029 0.45 2.81 3.24826 0.41 3.09 3.25164 0.37 3.37 3.25226 0.35 3.66 3.25224 0.32 3.94 3.25223 0.30 4.22 3.25222
TABLE VI:
Convergence of H at equilibrium in zero magnetic field. d (au) Range (au) E (Hartree) 1.00 4.00 -0.598735 0.80 5.00 -0.602098 0.67 6.00 -0.602555 0.57 7.00 -0.602616 Page 22
TABLE VII:
The calculation of H without a magnetic field under all angles and translations. The calculations are completed with the grid centered on the atom and the grid translated by 0.5 au. Each line was independently scanned and shown to be identical within tolerance. B= 0 au
Not translated Translated Reference θ E T R eq θ E T R eq θ E T R eq o -0.6026 2.00 0 o -0.6026 2.00 0 o -0.60263 2.00 45 o -0.6026 2.00 45 o -0.6026 2.00 45 o -0.60263 2.00 90 o -0.6026 2.00 90 -0.6026 2.00 90 o -0.60263 2.00 TABLE VIII:
A calculation of a translated and origin centered H in the presence of magnetic fields. The total binding energy, E T , equilibrium distance R eq at the specified angle θ . Comparison is made with reference [16]. For reference 1 au = 2.35052 x 10 Tesla
B= 1.0 au
Not translated Translated Reference θ E T R eq θ E T R eq θ E T R eq o -0.4750 1.75 0 o -0.4750 1.76 0 o -0.47496 1.75 45 o -0.4623 1.69 45 o -0.4622 1.69 45 o -0.45925 1.67 90 o -0.4506 1.64 90 o -0.4506 1.64 90 o -0.44956 1.63 B=10.0 au
Not translated Translated Reference θ E T R eq θ E T R eq θ E T R eq o o o o o o o o o Page 23
TABLE IX:
Results for Hydrogen, the H ion and the H ion in magnetic fields up to 10,000 au. Results are compared to known calculation from references [16]. Energy (Hartree) Reference Energy (Hartree) B (au) H H H H H H
1 -0.3311 -0.4750 - -0.3312 -0.47496 - 10 3.2522 2.8251 - 3.25225 2.82512 3.3042 100 46.21 44.86 45.68 46.2104 44.8545 45.6806 1000 492.34 488.56 489.57 492.341
Page 24 X. REFFERENCES [1] T. H. Dunning Jr, "
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen" , The Journal of Chemical Physics (2), 1007-1023 (1989). [2] W. J. Hehre, R. Ditchfield, and J. A. Pople, " Self—consistent molecular orbital methods. XII. Further extensions of gaussian—type basis sets for use in molecular orbital studies of organic molecules" , The Journal of Chemical Physics (5), 2257-2261 (1972). [3] G. Kresse and J. Furthmüller, " Software VASP, Vienna (1999)" , Phys. Rev. B (11), 169 (1996). [4] R. Dovesi, R. Orlando, B. Civalleri, C. Roetti, V. R. Saunders, and C. M. Zicovich-Wilson, " CRYSTAL: a computational tool for the ab initio study of the electronic properties of crystals" , Zeitschrift für Kristallographie (5/6/2005), 571-573 (2005). [5] M. Challacombe, E. Schwegler, C. J. Tymczak, A. N. Niklasson, N. Bock, K. Nemeth, V. Weber, C. K. Gan, and G. Hinkelman, FreeON code suite, 2009 http://en.wikipedia.org/wiki/FreeON http://savannah.nongnu.org/projects/freeon/ [6] Æ. Frisch, M. J. Frisch, and G. W. Trucks,
Gaussian 03 user's reference . (Gaussian, Incorporated, 2003). [7] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, "
QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials" , J. Phys: Condensed Matter (1), 395502 (2009). [8] C. Mundy, F. Mohamed, F. Schiffman, G. Tabacchi, H. Forbert, W. Kuo, J. Hutter, M. Krack, M. Iannuzzi, and M. McGrath, (2000). [9] D. T. Colbert and W. H. Miller, " A novel discrete variable representation for quantum mechanical reactive scattering via the S ‐ matrix Kohn method" , The Journal of Chemical Physics , 1982 (1992). [10] C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, " Introducing ONETEP: Linear-scaling density functional simulations on parallel computers" , The Journal of Chemical Physics (8), 084119 (2005). [11] C. E. Shannon and W. Weaver,
The mathematical theory of communication . (University of Illinois Press, Urbana, 1963). [12] R. G. Littlejohn, M. Cargo, T. Carrington Jr, K. A. Mitchell, and B. Poirier, "
A general framework for discrete variable representation basis sets" , The Journal of Chemical Physics , 8691 (2002). [13] J. C. Light, I. P. Hamilton, and J. V. Lill, "
Generalized discrete variable approximation in quantum mechanics@f[sup a]@f[sup )]" , The Journal of Chemical Physics (3), 1400-1409 (1985). [14] J. C. Light and T. Carrington Jr, " Discrete-variable representations and their utilization" , Advances in Chemical Physics , 263-310 (2000). [15] C. Cattani, "
Shannon Wavelets Theory" , Mathematical Problems in Engineering (2008).
Page 25 [16] A. V. Turbiner and J. C. L. Vieyra, "
One-electron molecular systems in a strong magnetic field" , Physics reports (6), 309-396 (2006). [17] F. Stenger,
Numerical methods based on sinc and analytic functions . (Springer, 1993). [18] S. G. Mallat, "
A theory for multiresolution signal decomposition: the wavelet representation" , Pattern Analysis and Machine Intelligence, IEEE Transactions on (7), 674-693 (1989). [19] M. Galassi, B. Gough, G. Jungman, J. Theiler, J. Davies, M. Booth, and F. Rossi, " The GNU Scientific Library Reference Manual, 2007" , 35 (2007). [20] E. Anderson,
LAPACK Users' guide . (Siam, 1999). [21] L. E. McMurchie and E. R. Davidson, "
One-and two-electron integrals over Cartesian Gaussian functions" , Journal of Computational Physics (2), 218-231 (1978). [22] K. Cho, T. A. Arias, J. D. Joannopoulos, and P. K. Lam, " Wavelets in electronic structure calculations" , Physical Review Letters (12), 1808 (1993). [23] J. P. Modisette, P. Nordlander, J. L. Kinsey, and B. R. Johnson, " Wavelet based in eigenvalue problems in quantum mechanics" , Chem. Phys. Lett. (5-6), 485-494 (1996). [24] L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, and R. Schneider, "
Daubechies wavelets as a basis set for density functional pseudopotential calculations" , The Journal of Chemical Physics (1), 014109-014114 (2008). [25] S. Goedecker and O. V. Ivanov, "
Frequency localization properties of the density matrix and its resulting hypersparsity in a wavelet representation" , Physical Review B (11), 7270 (1999). [26] T. A. Arias, " Multiresolution analysis of electronic structure: semicardinal and wavelet bases" , Reviews of Modern Physics (1), 267 (1999). [27] M. Holmstrom, " Solving Hyperbolic PDEs Using Interpolating Wavelets" , SIAM J. Sci. Comput. (2), 405-420 (1999). [28] B. Poirier and A. Salam, " Quantum dynamics calculations using symmetrized, orthogonal Weyl-Heisenberg wavelets with a phase space truncation scheme. III. Representations and calculations" , The Journal of Chemical Physics , 1704 (2004). [29] A. M. N. Niklasson, C. J. Tymczak, and H. Rˆder, "
Multiresolution density-matrix approach to electronic structure calculations" , Physical Review B (15), 155120 (2002). [30] C. J. Tymczak and X. Q. Wang, " Orthonormal Wavelet Bases for Quantum Molecular Dynamics" , Phys Rev Lett , 3654 (1997). [31] L. R. Ram -‐ Mohan, S. Saigal, D. Dossa, and J. Shertzer, "
The finite ‐ element method for energy eigenvalues of quantum mechanical systems" , Computers in Physics (1), 50-59 (1990). [32] B. I. Schneider and N. Nygaard, " Discrete variable representation for singular Hamiltonians" , Physical Review E (5), 056706 (2004). [33] A. Thirumalai and J. S. Heyl, " Hydrogen and helium atoms in strong magnetic fields" , Physical Review A (1), 012514 (2009). [34] T. C. Scott, M. Aubert-Frécon, and J. Grotendorst, " New approach for the electronic energies of the hydrogen molecular ion" , Chemical physics (2), 323-338 (2006).