Fully numerical electronic structure calculations on diatomic molecules in weak to strong magnetic fields
Susi Lehtola, Maria Dimitrova, and Dage Sundholm
TL;DR
This paper introduces fully numerical electronic structure calculations for diatomic molecules in magnetic fields, benchmarking against Gaussian basis sets and analyzing errors across various field strengths.
Contribution
It develops and applies a fully numerical method for diatomic molecules in magnetic fields and compares its accuracy to Gaussian basis set approaches.
Findings
Basis set truncation error increases with magnetic field strength.
Fully numerical results serve as benchmarks for Gaussian basis sets.
Reliable Gaussian basis results require high-angular-momentum functions at strong fields.
Abstract
We present fully numerical electronic structure calculations on diatomic molecules exposed to an external magnetic field at the unrestricted Hartree-Fock limit, using a modified version of a recently developed finite element program, HelFEM. We have performed benchmark calculations on a few low-lying states of H2, HeH+, LiH, BeH+, BH, and CH+ as a function of the strength of an external magnetic field parallel to the molecular axis. The employed magnetic fields are in the range of atomic units, where T. We have compared the results of the fully numerical calculations to ones obtained with the LONDON code using a large uncontracted gauge-including Cartesian Gaussian (GICG) basis set with exponents adopted from the Dunning aug-cc-pVTZ basis set. By comparison to the fully numerical results, we find that the basis set truncation error in the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40| Molecule | spin state | basis | configuration | field strength () |
|---|---|---|---|---|
| \ceH2 | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| \ceHeH+ | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| LiH | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| quintet | FEM | |||
| GTO | ||||
| Molecule | spin state | basis | configuration | field strength () |
|---|---|---|---|---|
| \ceBeH+ | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| quintet | FEM | |||
| GTO |
| Molecule | spin state | basis | configuration | field strength () |
|---|---|---|---|---|
| \ceBH | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| quintet | FEM | |||
| GTO | ||||
| \ceCH+ | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| quintet | FEM | |||
| GTO | ||||
| Spin state | triple- | quadruple- | quintuple- | |
|---|---|---|---|---|
| singlet | 1 | 0.57 | 0.16 | 0.03 |
| 2 | 24.16 | 4.00 | 1.77 | |
| triplet | 1 | 1.03 | 0.34 | 0.14 |
| 2 | 10.37 | 2.01 | 1.00 | |
| quintet | 1 | 2.48 | 0.75 | 0.24 |
| 2 | 103.30 | 29.00 | 25.04 |
| Molecule | spin state | basis | configuration | field strength () |
|---|---|---|---|---|
| \ceH2 | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| \ceHeH+ | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| LiH | singlet | FEM | ||
| GTO | ||||
| triplet | FEM | |||
| GTO | ||||
| quintet | FEM | |||
| GTO | ||||
| Spin state | basis | configuration | field strength () |
|---|---|---|---|
| singlet | FEM | ||
| GTO | |||
| triplet | FEM | ||
| GTO | |||
| quintet | FEM | ||
| GTO | |||
| Spin state | basis | configuration | field strength () |
|---|---|---|---|
| singlet | FEM | ||
| GTO | |||
| triplet | FEM | ||
| GTO | |||
| quintet | FEM | ||
| GTO | |||
| Spin state | basis | configuration | field strength () |
|---|---|---|---|
| singlet | FEM | ||
| GTO | |||
| triplet | FEM | ||
| GTO | |||
| quintet | FEM | ||
| GTO | |||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
\newref
eqname = eq. , names = eqs. \newreffigname = Fig. , names = Figs. \newrefsecname = Section , names = Sections
Fully numerical electronic structure calculations on diatomic molecules in weak to strong magnetic fields
\nameSusi Lehtola, Maria Dimitrova, and Dage Sundholm CONTACT: Susi Lehtola, email: [email protected] Dimitrova, email: [email protected] Sundholm, email: [email protected] Department of Chemistry, University of Helsinki, P.O. Box 55 (A. I. Virtasen aukio 1), FI-00014 University of Helsinki, Finland
Abstract
We present fully numerical electronic structure calculations on diatomic molecules exposed to an external magnetic field at the unrestricted Hartree–Fock limit, using a modified version of a recently developed finite element program, HelFEM. We have performed benchmark calculations on a few low-lying states of \ceH2, \ceHeH+, \ceLiH, \ceBeH+, \ceBH, and \ceCH+ as a function of the strength of an external magnetic field parallel to the molecular axis. The employed magnetic fields are in the range of atomic units, where T. We have compared the results of the fully numerical calculations to ones obtained with the London code using a large uncontracted gauge-including Cartesian Gaussian (GICG) basis set with exponents adopted from the Dunning aug-cc-pVTZ basis set. By comparison to the fully numerical results, we find that the basis set truncation error in the gauge-including Gaussian basis set is of the order of 1 kcal/mol at zero field, that the truncation error grows rapidly when the strength of the magnetic field increases, and that the largest basis set truncation error at exceeds 1000 kcal/mol. Studies in larger Gaussian basis sets suggest that reliable results can be obtained in GICG basis sets at fields stronger than , provided that a sufficient coverage of higher-angular-momentum functions is included in the basis set.
**Keywords Magnetic field, finite element, Hartree–Fock, intermediate regime, basis set truncation error.
**
1 Introduction
Because of the intrinsically weak strength of the magnetic interaction, most molecules exposed to even the strongest static magnetic fields available in laboratory conditions can be studied computationally with excellent accuracy using perturbation theory. For example, the most powerful continuous magnetic field obtained under laboratory conditions is produced by the 45-tesla magnet at the National High Magnetic Field Laboratory in Florida,1 for which meV that is small compared to the thermal energy , except at very small temperatures K.
Stronger fields, however, have to be considered explicitly in electronic structure calculations, as the magnetic field ceases to be a small perturbation and becomes of equal or greater importance than the Coulomb interactions between the electrons. The interest in such studies arose when Kemp et al. discovered circularly polarised light coming from a white dwarf, and estimated the magnetic field strength to be about T at the surface.2 Zeeman splittings that are consistent with the lines of light atoms and molecules such as H, He, O, CH, and \ceC2 have also been observed in the atmospheres of magnetic white dwarfs.3; 4; 5; 6; 7; 8 In addition to these applications in astrophysics, the modeling of matter in extreme magnetic fields can be beneficial also for other purposes, as e.g. the behavior of some solid-state systems with a high dielectric constant mimics that of atoms in strong magnetic fields.9; 10
However, the computational modeling of the electronic structure of matter in strong magnetic fields is challenging, in general, since the magnetic interaction confines the electrons in directions perpendicular to the field. Furthermore, since the magnetic field couples explicitly to the spin, also the ground-state configuration depends on the field strength, yielding a rich state diagram with several state crossings. But, in extreme conditions such as in the atmospheres of highly magnetized neutron stars a.k.a. magnetars where the field strength can reach billions of tesla,11 electronic structure calculations become again easier, as the magnetic interaction becomes dominant and Coulomb interactions can be treated perturbatively to good accuracy. This is known as the Landau regime in which the orbitals are elongated along the magnetic field axis due to the strong perpendicular confinement, and the energy barrier to two- and three-dimensional molecular structures is so large that only one-dimensional chains can be formed.12 Several computational studies on the structure of atoms and chains of atoms in the strong-field regime have been reported.13; 14; 15; 16; 17; 18; 19; 20
However, the more interesting and computationally challenging situations appear when the Coulomb and the magnetic-field interactions are of similar magnitude; this is called the intermediate regime,21 where the choice of the basis set is difficult. The traditional linear combination of atomic orbitals (LCAO) approach, in which the atomic orbitals are by convention isotropic, works well for zero-field calculations. When the magnetic interactions become dominant over Coulomb interactions, Landau orbitals are often employed. However, in the presently investigated intermediate regime, neither the traditional LCAO expansion nor Landau orbitals are expected to yield a fast convergence to the basis set limit.
Gaussian-type orbital (GTO) basis sets are the common choice in calculations at zero field, because the electron interaction integrals can be calculated in a GTO basis set in a straightforward manner. Calculations using GTOs yield good results even at finite magnetic field strengths, provided that a gauge factor is explicitly included in the basis set; that is, if gauge-including atomic orbitals (GIAO) a.k.a. London orbitals are used.22; 23; 24 Calculations with gauge-including GTOs have shown that significant changes in molecular structure may be observed in the intermediate regime:25 for instance, a new type of paramagnetic bond can be created by the interaction of two atoms with high spin multiplicity.26 Lange et al. have shown that \ceH2 exists as a triplet-state molecule when the magnetic field is stronger than 175 000 T, with the molecule strongly preferring an orientation perpendicular to the field.26 This paramagnetic bonding mechanism explains why Xu et al. could observe molecular hydrogen in white dwarf atmospheres,27 even though the temperatures therein may reach over K.
However, the use of conventional isotropic gauge-including GTOs might lead to large basis set truncation errors at finite fields, as the electronic structure of atoms deforms continuously from spherical symmetry at zero field to needle-like structures at strong fields, and as the latter shape is not readily expandable in such a basis set. Therefore, finite-field calculations are typically performed only up to about one atomic unit of field strength T in this kind of basis set. Furthermore, uncontracted basis sets are used for a better description of the deformation of the orbitals caused by the magnetic field.28; 29; 30; 31
Although the problems encountered with isotropic GTO basis sets in the description of needle-like orbitals can be circumvented by employing anisotropic GTO basis sets,32; 33 which have been successfully used to study atoms34; 35; 36; 37; 38 and molecules39; 40; 41 in strong magnetic fields, they have some shortcomings as well. It is more difficult to optimize their exponents than those of conventional isotropic basis sets.42; 43; 44 Anisotropic basis sets are also considerably larger than isotropic ones, because the latter consist of GTOs with the same exponents in the three spatial directions, whereas the former employ a set of and exponents parallel and perpendicular to the field, respectively, yielding basis functions. More general approaches, such as using a full matrix of exponents in , are possible as well, which lead to an even larger number of basis functions. Moreover, the use of anisotropic exponents requires dedicated software,32; 33; 34 because the calculation of one- and two-electron integrals becomes more complicated and cannot be performed with standard molecular integrals codes.
Yet, all the problems inherent in either GTO basis set approaches can be circumvented by switching to the use of a fully numerical basis set. Fully numerical approaches typically allow for adaptively refined basis sets, and a systematical approach to the complete basis set limit. Even highly anisotropic problems can be readily treated with such approaches. Furthermore, if a practically complete basis set is used – as is typically done in numerical approaches – the use of gauge-including basis functions is not necessary for maintaining the gauge invariance of the magnetic field. For instance, fully numerical calculations on atoms and \ceH2+ in weak to strong magnetic fields have been reported by Ivanov and Schmelcher,45; 46; 47; 48; 49; 50; 51 whereas calculations of atoms in magnetar-level magnetic fields have been reported by Schimeczek et al.52 Static magnetic properties employing linear response at the self-consistent field level of theory have also been recently reported by Jensen et al.53
Here, we present a fully numerical approach for electronic structure calculations on diatomic molecules with an explicit, finite magnetic field along the molecular axis. The matrix elements of the magnetic field interaction Hamiltonian have been derived and implemented in the recently published, freely available HelFEM program54 for fully numerical calculations on atoms55 and diatomic molecules.56 HelFEM uses a variational method for solving the self-consistent field equations of Hartree–Fock (HF) or Kohn–Sham57 density functional theory (DFT),58 and supports hundreds of exchange-correlation functionals via an interface to the Libxc library.59 The present approach is ideal for studies with magnetic fields strengths belonging to the intermediate regime, as Coulomb and magnetic interactions are treated self-consistently on the same footing in HelFEM.
The partial wave expansion originally proposed by McCullough60; 61 is used in HelFEM for calculations on diatomic molecules, in which the “radial” direction62 is further described with one-dimensional finite element shape functions; see ref. 56 and \secrefnummeth for more details. Analogous grid-based approaches have also been pursued previously; they have been reviewed in ref. 62. The complete basis set (CBS) limit can be reached systematically in a variational fashion with the HelFEM program, as the results converge smoothly and monotonically when the size of the basis set is increased,55; 56 thus enabling electronic structure studies even at strong magnetic fields. The prolate spheroidal coordinate system used in the present work has also been previously used in studies of the \ceH2+ ion in a strong magnetic field parallel to the molecular axis.63; 64
Numerical calculations on diatomic molecules have a long history of being used for establishing basis set limits and LCAO basis set truncation errors in the absence of electromagnetic fields as well as for finite electric fields, as has been reviewed elsewhere.62 In the present work, we study the complete basis set (CBS) limits of the singlet, triplet, and quintet states of \ceH2, \ceHeH+, \ceLiH, \ceBeH+, \ceBH and \ceCH+ in weak to strong magnetic fields by performing fully numerical calculations with HelFEM. Guided by the CBS limit results obtained with HelFEM, we assess the accuracy of commonly-used gauge-including Cartesian Gaussian (CICG) basis sets, as the question of the reliability of GTO basis sets is an important issue for studies at finite magnetic fields. The Gaussian-basis calculations are done with the London program.
The study is undertaken at the unrestricted HF (UHF) level of theory, which is sufficient for the present purpose of establishing the basis set truncation error in commonly employed basis sets. Although both HelFEM and London also support DFT calculations, the basis set requirements of HF and DFT are known to be close to identical. UHF calculations allow for an unbiased comparison of the fully numerical and GTO approaches.
The organization of the present work is the following. The basic theory is briefly presented in \secrefmagfield. The computational approach is presented in \secrefcompmeth: the numerical approach as implemented in HelFEM is described in \secrefnummeth, the gauge-including isotropic Cartesian Gaussian basis set methods as implemented in the London code65 are presented in \secrefgtometh, and the details on the calculations are given in \secrefcompdetails. The results are discussed in \secrefresults, and the main conclusions thereof are summarized in \secrefsummary. Atomic units are employed, if not specified otherwise.
2 Quantum Systems in Finite Magnetic Fields
As is well known from Maxwell’s equations, any magnetic field can be defined in terms of a vector potential as . The presence of an electromagnetic field changes the Hamiltonian for an electron as
[TABLE]
where is the field-free molecular Hamiltonian and is the momentum operator which commutes with for a uniform magnetic field. In the rotational-symmetry preserving symmetric gauge, the uniform magnetic field is generated by the vector potential
[TABLE]
However, the vector potential has a gauge degree of freedom : any potential generates the same magnetic field as . The gauge of can be chosen freely for calculations, with some choices being more favorable for computation than others. Typically, the Coulomb gauge
[TABLE]
is used, as it leads to the least number of terms in the equations. This choice still leaves the freedom to choose e.g. a gauge origin , as replacing with in (LABEL:Auniform) again generates the same field . The inclusion of the gauge origin dependence yields the final expression for the vector potential
[TABLE]
Including this dependence on the gauge origin of in the Hamiltonian, the many-electron Hamiltonian can now be written in terms of the sums of the orbital and spin angular momenta
[TABLE]
where is the total spin, is the position of the :th electron with respect to the global gauge origin , and is the canonical angular momentum.
Now, choosing and the magnetic field to coincide with the axis as , (LABEL:Hamvec) simplifies to
[TABLE]
as then only the projections of the orbital and spin angular momentum operators and , respectively, remain in the orbital-Zeeman () and spin-Zeeman () terms. Next, as the magnetic field direction coincides with the molecular axis, the matrix elements of (LABEL:HamZ) can be made real. The Hamiltonian then has no explicit dependence on the azimuth angle around the molecular axis, and by Noether’s theorem the -momentum is conserved. This leads to a blocking of the molecular orbitals, which have an azimuth-angular dependence of with integer values of , analogously to the field-free cases discussed in refs. 55, 56, and 62.
The terms in (LABEL:Hamvec,_HamZ) that are linear in describe paramagnetic response, while the terms are diamagnetic. The paramagnetic terms can either increase or decrease the energy as the field strength grows, whereas the diamagnetic term is positive, always increasing the energy of all electronic systems. All systems become diamagnetic at strong fields due to the quadratic character of the diamagnetic terms.66 The diamagnetic term leads to a confinement in the directions perpendicular to the magnetic field, which causes atoms and molecules to become elongated in the parallel direction in strong magnetic fields.
Due to the orbital-Zeeman and spin-Zeeman effects, electron configurations change dramatically with increasing field strength. High-spin states with occupied -spin () orbitals are favored in strong magnetic fields, and orbitals with large negative angular momentum quantum numbers are stabilized. For instance, for atoms in strong magnetic fields, the orbital can lie energetically much lower than the orbital of the same electronic shell.28 Here and in the rest of the manuscript, we use the sign convention that the energy of orbitals with and spin are lowered by the external magnetic field due to the orbital-Zeeman and spin-Zeeman effects, respectively.
3 Computational Methods
3.1 Numerical calculations
The numerical approach used in HelFEM follows the traditional variational basis set approach of quantum chemistry; that is, the program solves the Roothaan–Hall/Pople–Nesbet equations67; 68
[TABLE]
where and are the Fock matrix and molecular orbital coefficients for spin , are the corresponding orbital energies, and is the overlap matrix. However, instead of a basis set composed of atomic orbital basis functions with global support, HelFEM employs a basis set of finite element shape functions, which are local; see ref. 55 for details on the presently employed basis functions. The diatomic approach, in which the “radial” functions are expanded in finite-element functions and the “angular” part is expanded in spherical harmonics, has been described in ref. 56. In the present work, HelFEM has been modified such that the magnetic field interaction terms of (LABEL:HamZ) are added into the core Hamiltonian and Fock matrices, as is described below. The implementation and calculations otherwise follow the procedures described in refs. 55 and 56, where the computational approach is described in more detail.
Although we consider only molecular calculations in this work, we also describe the implementation of magnetic fields for atomic calculations in HelFEM. In both cases, the component of the angular momentum is given by
[TABLE]
which is easily recognized by noting that in either case the dependence of the orbitals is of the form ,55; 56 where is the magnetic quantum number.
3.1.1 Atoms
The numerical basis set for atomic calculations is55
[TABLE]
where are the radial shape functions and are spherical harmonics. The matrix elements of (LABEL:HamZ) are then obtained as
[TABLE]
and
[TABLE]
as the parabolic confining term is given by
[TABLE]
in spherical polar coordinates. The only contribution that has to be added to the expressions in ref. 55 is the term, which can be handled via Gaunt coefficients with the expansion
[TABLE]
For details on the evaluation of the integrals in (LABEL:Lz-ij,_x2y2-atdef), see ref. 55.
3.1.2 Diatomic molecules
For diatomic molecules, the numerical basis set is of the form56
[TABLE]
where are transformed prolate spheroidal coordinates. While maintains the same form as (LABEL:Lz-ij), , the matrix element for the parabolic confining term is obtained as
[TABLE]
which can again be evaluated using (LABEL:sin2) and the additional expansion
[TABLE]
For details on the evaluation of the integrals in and (LABEL:x2y2-diat), see ref. 56.
3.2 Gaussian basis set calculations
Gauge-independent electronic structure LCAO calculations on atoms and molecules exposed to an external magnetic field can be performed by modifying the basis set to include an explicit gauge transformation factor22; 23; 24
[TABLE]
resulting in GIAOs, or London orbitals. The complex prefactor in (LABEL:psiGauge) corresponds to a gauge-origin transformation from to ((LABEL:vecPot)), which makes the wave function invariant to first-order changes in the gauge origin.
The London program employs GIAOs formed of Cartesian GTOs.69 A primitive Cartesian GTO centered at is given by
[TABLE]
where are integers, the sum of which corresponds to the angular character of the shell (e.g. for , , , and shells, respectively), is the exponent, is the corresponding normalization factor, and is the square of the distance from the center of the basis function, which typically resides at a nucleus.
3.3 Computational details
We have studied the total energy of a few low-lying states of six diatomic molecules as a function of the strength of the external magnetic field. Since the aim of this study is to assess the quality of calculations employing commonly-used Gaussian basis sets, the unrestricted Hartree–Fock level of theory is sufficient, and fixed internuclear distances () have been employed. The studied molecules include \ceH2 (), \ceHeH+ , \ceLiH (), \ceBeH+ (), \ceBH ), and \ceCH+ () in their singlet, triplet, and quintet spin states; the internuclear distances () were chosen based on the equilibrium bond length of the molecular species for at the HF/def2-SVP level of theory.70 Correlation effects or the fact that the equilibrium bond distances shrink with increasing strength of the magnetic field along the molecular axis have not been considered in this work.
As was already stated in \secrefgtometh, the GTO calculations were performed with the London program developed by Tellgren et al.65; 69 An uncontracted gauge-including Cartesian basis set was used, with isotropic exponents adopted from the triple- correlation consistent basis set augmented with diffuse functions (aug-cc-pVTZ).71; 72 This basis set was chosen for the present work, as it has also been used in previous studies.26 A linear dependency threshold of was used in the canonical basis set orthonormalization procedure. The direct inversion of the iterative subspace (DIIS) method73; 74 was used for accelerating the convergence of the self-consistent field procedure. An orbital rotation convergence threshold of was used in the London calculations. In problematic cases, the lowest state was obtained by restarting the calculations from converged density matrices for the desired configuration for a different field strength.
The fully numerical calculations were carried out using the HelFEM program, which employs finite-element basis functions, as was detailed above in \secrefnummeth. The DIIS method was likewise used in the HelFEM calculations, in combination with the ADIIS algorithm.75 An orbital rotation convergence threshold of was used in the HelFEM calculations.
For simplicity, the same numerical basis set was used for all molecules. The numerical basis set was determined at zero field by the proxy approach described in ref. 56. Taking the largest basis set parameters for each molecule, it was found that energies beyond microhartree accuracy at zero field are obtained for all molecules with a linear element grid containing three 15-node radial Lobatto elements, with and . A closer study showed that \ceBH and \ceCH+ have a ground state for some field strengths. The initial basis set was therefore increased to include orbitals with , which was determined by extrapolating from and .
Although magnetic flux density is traditionally measured in units of tesla, in computational studies the field strength is instead expressed in atomic units . One atomic unit of magnetic flux density is equal to T, where J, C, m/s and the inverse fine structure constant is .
The investigated magnetic field strengths are in the range of in steps of . The GTO results obtained for magnetic fields larger than may not be entirely reliable, because the employed basis set is not able to span the orbitals accurately at strong fields. The extreme case is certainly beyond the scope of the isotropic GTO approach studied in the present work. Due to the prevalence of state crossings in the calculations with magnetic field, we have explicitly tracked several electronic configurations for the studied molecules by fixing the symmetry of the occupied orbitals in the fully numerical calculations, as described in refs. 55 and 56. The ground state energy has then been extracted by identifying the lowest energy at each value of the magnetic field strength.
In order to study the basis set convergence at finite magnetic field, larger numerical basis sets were formed by increasing the number of partial waves by simultaneously in each subchannel, whereas five radial elements were used in all calculations. Although the accuracy of the fully numerical calculations was found to decrease monotonically with increasing field strength, the basis set convergence was found to be exponential with the basis set size. While the numerical calculations are fully converged to the basis set limit at , calculations with the largest basis set consisting of five radial elements and , , are estimated to have millihartree accuracy at , which is sufficient for the purpose of the present study.
4 Results
4.1 General results
Total energies as a function of the magnetic field strength calculated using GTO and numerical basis sets are shown in \figrefabsE. The obtained ground state configurations are given in \tabrefgs-B-a for \ceH2, \ceHeH+, LiH, in \tabrefgs-B-b for \ceBeH+, and in \tabrefgs-B-c for BH and \ceCH+. The GTO basis set truncation errors in \figrefdiffE were obtained by comparing the total energies of the GTO calculations with London with the converged fully numerical values from with the HelFEM code.
Although many states have been explicitly considered in the present work, the data in \figrefabsE, diffE correspond only to the lowest energy at each value of the field. Plots of the energies of the individual configurations are shown in the Supporting Information, in addition to detailed data on how the energy ordering of the , , and orbitals changes with magnetic field.
The curves in \figrefdiffE display oscillations, which can tentatively be attributed to the GTO basis set. As the same basis set is used for all field strengths, the basis set truncation error can be smaller for some values, for which the GTO expansion is better able to describe the electronic wave function than for other values of . This basis set artifact can result in the observed oscillatory behavior.
The basis set errors of the GTO calculations at zero field are of the order of 1 kcal/mol (1 eV = 23.06054801(14) kcal/mol) with the employed gauge-including fully uncontracted Cartesian triple- basis set, but the accuracy decays rapidly with increasing field strength. The poor performance of GTO basis sets for absolute energies is well-known, whereas relative energies are usually more accurate due to cancellation of errors. However, in the presence of strong magnetic fields, relative energies such as spin-state splittings may also become inaccurate, since the accuracy of the GTO calculations varies strongly with the spin state, as will be seen below.
4.2 Singlet states
The difference of the GTO and the numerical energies is 0.5–2.3 kcal/mol at for all studied singlet states. The GTO truncation errors approach 1000 kcal/mol at strong magnetic fields for all but the lightest molecules, \ceH2 and \ceHeH+.
The ground-state orbital occupations of the London and HelFEM calculations agree well for the singlet states. The problematic cases turned out to be \ceLiH and \ceBeH+. \ceH2 and \ceHeH+ have the electron configuration in the investigated magnetic field strength of . In the GTO calculations, the energy of the singlet state of \ceLiH has a kink at , because the electron configuration changes from to a state where the highest occupied and spin-orbitals have neither nor symmetry. In contrast, the HelFEM calculations employ axial symmetry that prevents symmetry breaking. In the fully numerical calculations on \ceLiH, the orbital configuration changes from to at . This configuration was not obtained in the London calculations.
Symmetry breaking is a well-known artifact of the Hartree–Fock method that originates from the neglect of electron correlation. A full wave function description includes all configurations which may be spatially dissimilar, the superposition of which, however, has the same symmetry as the system itself. As Hartree–Fock considers only one of these configurations, which does not reflect the full symmetry of the system, it is sometimes advantageous to further break the symmetries that are fulfilled by the exact wave function. Because the character of the symmetry-broken GTO calculation does not match that of the numerical calculations, the data corresponding to symmetry-broken configurations have been excluded from the comparison.
Calculations on the singlet state of \ceBeH+ yield the same orbital occupation with both methods. The crossing between the and the states occurs at – significantly earlier than for \ceLiH. The GTO calculations yielded a symmetry-broken state at . By tracing the symmetry-broken solution down to weaker fields, it turned out to be lower in energy than the state with the configuration all the way down to , where it crosses the state. As the symmetry-broken state is again a Hartree–Fock artifact, we did not consider it further and excluded it from the graphs.
The singlet state of \ceBH has a configuration at weak fields. Both the London and HelFEM calculations predict that a state crossing occurs at , where a orbital becomes the highest occupied orbital.
The \ceCH+ cation has a weak-field ground state with occupied orbitals of only symmetry. For field strengths , the ground-state configuration becomes , as for \ceBH. Symmetry analysis of the London calculations showed that the orbital becomes more stable than the orbital at fields stronger than .
4.3 Triplet states
The total energy of the triplet state of the molecules calculated with the two approaches agree well up to about . When the field strength is increased further, the basis set truncation error increases again and approaches 1000 kcal/mol at .
The \ceH2 molecule did not behave well in the GTO calculations, probably due to lack of basis functions needed for an accurate description of the orbital. The basis set truncation error is 11.5 kcal/mol at . For the other molecules, the basis set truncation errors lie in the range of 0.2 – 3.3 kcal/mol for . The monotonically increasing basis set truncation errors reach 30 kcal/mol at . At strong magnetic fields, the errors reach hundreds of kcal/mol for \ceH2 and \ceHeH+, and 1300 kcal/mol for \ceBH.
State crossings occur for many of the triplet states. In the HelFEM calculations, the triplet state of \ceH2 exhibits crossings between the and configurations at . As was mentioned above, the state is not reproduced accurately by the GTO calculation, which can tentatively be attributed to the insufficient number of functions in the used basis set for hydrogen. This also leads to a different state-crossing point compared to the numerical calculation.
The calculations show that there is a state crossing at for the triplet state of \ceBeH+ and at for the triplet state of \ceLiH, where the configuration becomes the ground state.
The triplet state of \ceBH was found to be challenging. The London calculations yielded an energy of at , which is lower than the energy of obtained in the fully numerical calculation. Field-free calculations with a development version of Q-Chem 5.1.177 yielded the same energy as London, whereas the total energy obtained with the Erkale78; 79 and Psi480 codes for a symmetry-abiding wave function is higher than the one obtained in the fully numerical calculation. Despite the significant energy difference, the Psi4 and Q-Chem calculations converged onto true local minima. The different energies that were obtained with the London, Erkale, Psi4 and Q-Chem programs suggest that the calculations on the triplet state of \ceBH converged to different solutions, even though the orbital configurations supposedly match. However, once the magnetic field is switched on, the problem with multiple solutions vanishes. The same solutions are obtained with the London and HelFEM codes for .
The electron configuration of the triplet state of \ceBH is in weak fields. The configuration of the electrons changes from to at according to the London calculations and at in the HelFEM calculations. The state with the configuration crosses the state at in the fully numerical calculations; this state crossing was not obtained in the GTO calculations.
The triplet state of \ceCH+ behaves in the same way as the triplet state of \ceBH, except at . The weak field configuration is . In the GTO calculations, a state crossing occurs at , where the state with the becomes the ground state. This state crossing occurs already at in the fully numerical calculations.
4.4 Quintet states
The quintet states are well described by the GTO calculations in the zero-field case, while the basis set truncation errors increase significantly at finite field. Large oscillations in the basis set truncation error are observed for \ceBH and \ceCH+.
In the GTO calculations, the ground-state electron configuration of the quintet state of \ceLiH is at fields weaker than ; this configuration is an excited state in the fully numerical calculations. The configuration is the ground state of \ceLiH in the range of in the GTO calculations, whereas a state crossing from the configuration to occurs already at in the fully numerical calculations. This state containing a orbital was not found in the GTO calculations. Instead, a second state crossing is seen in the GTO calculations at , where the configuration becomes the ground state again. The second state crossing is unexpected and is probably another artifact due to the incomplete GTO basis set.
The quintet state of \ceBeH+ has a configuration in the whole range of the investigated magnetic field strengths in the GTO calculations. The ground state in the fully numerical calculations, however, changes to at , similarly to the case of LiH.
Many ground-state electron configurations were obtained for \ceBH in the GTO calculations, whereas the numerical calculations yielded only two ground-state configurations. For , the highest occupied orbital is a destabilized orbital, whose orbital energy increases with increasing strength of the magnetic field. The obtained electron configuration at weak magnetic fields is in both the HelFEM and London calculations. For magnetic fields stronger than , the HelFEM calculations yielded a ground-state configuration of . In the GTO calculations, the ground-state configuration in the range of is , whose orbital is stabilized when increasing the magnetic field strength. For , the obtained configuration is . The two latter states are excited states in the fully numerical calculations. The sharp peak at in the curve for the quintet state in \figrefdiffE appears to originate from variations in the basis set truncation error with respect to the strength of the magnetic field.
For , a ground-state configuration is obtained in both the London and HelFEM calculations, whereas at stronger fields the ground-state configuration of the GTO calculations is again . Due to the orbital-Zeeman term in the Hamiltonian, it is unexpected that a high-angular-momentum orbital with becomes less stable than a orbital with , when increasing the strength of the magnetic field. Thus, the second state crossing in the GTO calculations is most likely an artifact due to the employed GTO basis set.
The ground-state configuration of the quintet state of \ceCH+ is at weak fields. The first state crossing to occurs at in the HelFEM calculations and at in the London calculations. The small peak at in the curve for the quintet state in \figrefdiffE appears to be caused by variations in the basis set quality with respect to the strength of the external magnetic field. A second state crossing is obtained in the GTO calculations at , where the configuration becomes the lowest one. This is again most likely an artifact due to the employed GTO basis set as in the \ceBH case.
4.5 Basis set study
Basis set truncation errors were studied by performing calculations on \ceBeH+ as the representative molecule due to its rich state diagram. Preliminary calculations showed that higher-angular-momentum functions become more important in stronger magnetic fields, and that deviations from the numerical reference energies became much smaller when using larger basis sets.
The influence of higher-angular momentum functions was assessed by performing calculations using uncontracted aug-cc-pVQZ and aug-cc-pV5Z basis sets without or functions, in order to make the calculations more tractable. This yields triple-, quadruple- and quintuple- basis sets with the compositions of , , and for H, and , , and for Be, respectively. The results of the calculations in these basis sets are shown in \tabrefbste. The basis-set truncation errors are larger for stronger magnetic fields, as expected.
The calculations show that the triple- basis set describes the singlet and triplet states of \ceBeH+ better than its quintet state. The basis is not large enough for the quintet state, as it fails to reproduce the occupied orbital, which is, however, found with the quadruple- and quintuple- basis sets. This manifests as the huge error – over 100 kcal/mol – at . But, the energy difference from the numerical reference value for the same configuration is only 19.09 kcal/mol.
At strong magnetic fields, the ground state might have occupied orbitals with larger quantum numbers. However, states with occupied orbitals that have values larger than 2 ( orbitals) have not been investigated in the present work. The benchmark study suggests that rather accurate energies can be obtained for molecules in magnetic fields stronger than , provided that a GTO basis set providing sufficient coverage for the occupied orbitals in its range of exponents is used. For instance, a proper description of the orbital in the quintet state of \ceBeH+ requires a large set of functions, as well as higher-angular-momentum functions to describe the polarization of the orbital.
5 Summary and Conclusion
We have performed fully numerical Hartree–Fock (HF) calculations on diatomic molecules exposed to an external magnetic field along the molecular axis with the HelFEM program. The main aim of the study was to assess the basis set truncation errors in unrestricted HF energies obtained using large gauge-including basis sets of uncontracted Gaussian-type orbitals (GTO). Although calculations at the density functional theory (DFT) level can also be performed with HelFEM, for simplicity we chose the HF level for the present study, since HF and DFT calculations have similar basis set requirements.
We performed calculations on a few low-lying states of \ceH2, \ceHeH+, \ceLiH, \ceBeH+, \ceBH, and \ceCH+ as a function of the magnetic field strength in the range of , using the recently published HelFEM program55; 56; 54 that has been extended to calculations on atoms and diatomic molecules in finite magnetic fields in the present work. Due to state crossings, several electronic states were studied as a function of the magnetic field strength. As the aim was to study basis set effects, fixed internuclear distances were employed, even though bond lengths are known to decrease significantly when increasing the strength of the magnetic field.
In the present work, we compared energies calculated using the numerical approach implemented in the HelFEM code with values calculated with the London program, employing London- a.k.a. gauge-including atomic orbitals (GIAO) formed from commonly-used primitive Gaussian basis sets. As the HelFEM calculations employ local basis sets close to the complete basis set limit, gauge factors such as those used in London are not necessary in the fully numerical approach. The energies obtained in the fully numerical UHF calculations can be used in other studies as reference data.
The calculations show that the Cartesian forms of the uncontracted aug-cc-pVTZ basis sets71; 72 have basis set truncation errors of the order of 1 kcal/mol in the absence of the external magnetic field, whereas the truncation errors of the gauge-including basis become significantly larger when the strength of the magnetic field is increased. The largest basis set truncation error at is more than 1000 kcal/mol. The calculations also show oscillations in the energy as a function of the magnetic field, which can be possibly explained by the basis set truncation error being smaller for some magnetic field strengths, for which the exponents happen to be closer to optimal values than at other field strengths.
Calculations on \ceBeH+ with larger basis sets showed that the basis set truncation errors can be made significantly smaller by employing basis sets with more high-angular-momentum functions. The problem of determining optimal atomic-orbital basis sets for finite fields could be investigated in future work.
The comparison of the energies calculated using the two approaches revealed lots of technical and formal difficulties. The ground-state electron configurations change when increasing the strength of the external magnetic field. When symmetry restrictions are not employed, optimization algorithms may not always be able to find the lowest-lying configuration, which leads to discontinuities in the energy as a function of the magnetic field strength. The discontinuities can, however, be resolved by tracking the various states by reading in converged orbitals for the wanted configuration at a different value of .
We also discovered that in some cases the lowest state corresponds to a symmetry-broken solution. Symmetry-broken solutions were a problem only in the GTO calculations, as symmetry restrictions are not supported in the London program, whereas in the fully numerical HelFEM program the symmetry of the occupied orbitals were explicitly enforced. HelFEM also supports calculations of states with broken symmetry. However, locating symmetry-broken solutions in the extended basis set is not trivial.55; 56
Acknowledgment
We thank Trygve Helgaker and Erik Tellgren for a copy of the London code, and Stella Stopkowicz and Florian Hampe for their tool for determining orbital symmetry. This work was supported by The Academy of Finland through projects 311149 and 314821 as well as by the Magnus Ehrnrooth Foundation. The authors acknowledge CSC – the Finnish IT Center for Science as well as the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) for computer time.
Supporting Information
The total energy as a function of the magnetic field strength for all states obtained in the calculations with the finite-element (FEM) calculations with HelFEM, and the gauge-including Gaussian-type orbital (GTO) calculations with London are shown in \figrangerefHH-1HH-3 for \ceH2, in \figrangerefHeH+-1HeH+-3 for \ceHeH+, in \figrangerefLiH-1LiH-5 for \ceLiH, in \figrangerefBeH+-1BeH+-5 for \ceBeH+, in \figrangerefBH-1BH-5 for \ceBH, and in \figrangerefCH+-1CH+-5 for \ceCH+. The resulting ground-state configurations in the FEM and GTO calculations have been identified in \tabrefgs-B-light for \ceH2, \ceHeH+, and \ceLiH, and in \tabrefgs-B-BeH+,gs-B-BH,gs-B-CH+ for \ceBeH+, \ceBH and \ceCH+, respectively.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 11 J.R. Miller, IEEE Trans. Appl. Supercond. 13 , 1385–1390 (2003).
- 22 J.C. Kemp, J.B. Swedlund, J.D. Landstreet and J.R.P. Angel, Astrophys. J. 161 , L 77–L 79 (1970).
- 33 J.R.P. Angel, Astrophys. J. 216 , 1–17 (1977).
- 44 G.D. Schmidt, P. Bergeron and B. Fegley, Astrophys. J. 443 , 274 (1995).
- 55 S. Jordan, P. Schmelcher, W. Becken and W. Schweizer, Astron. Astrophys. 36 , L 33–L 36 (1998).
- 66 D.T. Wickramasinghe and L. Ferrario, Publ. Astron. Soc. Pacific 112 , 873–924 (2000).
- 77 D.T. Wickramasinghe, G. Schmidt, L. Ferrario and S. Vennes, Mon. Not. R. Astron. Soc. 332 , 29–33 (2002).
- 88 J. Liebert, H.C. Harris, C.C. Dahn, G.D. Schmidt, S.J. Kleinman, A. Nitta, J. Krzesiski, D. Eisenstein, J.A. Smith, P. Szkody, S. Hawley, S.F. Anderson, J. Brinkmann, M.J. Collinge, X. Fan, P.B. Hall, G.R. Knapp, D.Q. Lamb, B. Margon, D.P. Schneider and N. Silvestri, Astron. J. 126 , 2521–2528 (2003).