Domain Decomposition Method for the N-body Time-Independent and Time-Dependent Schrödinger Equation

Domain Decomposition Method for the N-body Time-Independent and Time-Dependent Schrödinger Equation

E. Lorin School of Mathematics and Statistics, Carleton University, Ottawa, Canada, K1S 5B6 Centre de Recherches Mathématiques, Université de Montréal, Montréal, Canada, H3T 1J4

This paper is devoted to the derivation of a pleasingly parallel Galerkin method for the time-independent -body Schrödinger equation, and its time-dependent version modeling molecules subject to an external electric field BAN (); gauge (); CCT (). In this goal, we develop a Schwarz Waveform Relaxation (SWR) Domain Decomposition Method (DDM) for the -body Schrödinger equation. In order to optimize the efficiency and accuracy of the overall algorithm, i) we use mollifiers to regularize the singular potentials and to approximate the Schrödinger Hamiltonian, ii) we select appropriate orbitals, and iii) we carefully derive and approximate the SWR transmission conditions. Some low dimensional numerical experiments are presented to illustrate the methodology.

-body Schrödinger equation, domain decomposition method, mollifiers, parallel computing.

1 Introduction

This paper is devoted to the derivation of a pleasingly parallel real-space algorithm for solving the -body Schrödinger equation. It is well-known that the numerical computation to the solution to this equation faces the curse of dimensionality, as it requires computations in a -dimensional space. Even for 2 electrons (), smart parallel numerical algorithms must then be developed in order to tackle this problem. At the discrete level, the dimension of the Hamiltonian is basically dependent on the chosen basis functions: the less information is contained in the basis functions, the larger the dimension of the approximate Hamiltonian. The “worst-case scenario”, would be a finite difference approximation. However, elaborated alternatives exist such as the Full Configuration Interaction (FCI) ostlund (), which requires highly non-trivial basis functions, and which then allows for the construction of relatively compact discrete Hamiltonians. The time-independent and dependent Schrödinger equations are computed using algorithms (linear system and eigenvalue solvers) basically requiring many high-dimensional matrix-vector products. The parallelization of such algorithms is a major research field, and numerous parallel libraries exist among which, we can for instance cite lapack, arpack, sparselib, petsc, iml,.... However, the performance of the parallel implementation is often far from ideal, and other approaches should be explored. In this goal, we propose a Schwarz domain decomposition method for solving the -body Schrödinger equation. The principle consists of solving a large number of Schrödinger equations on “small” spatial domains. The interest is double. First, the parallel implementation is expected to be highly efficient, as local Schrödinger equations are solved independently and their corresponding solution is connected to the others only through the so-called transmission conditions DoleanBook (). Secondly, we can benefit from a scaling effect, as the computational complexity of real-space numerical (time-independent or dependent) Schrödinger equation solvers, is usually polynomial in time. We then solve in parallel, several small linear systems associated to local Schrödinger equations, rather than a large one. The price to pay though, is the need for computing several times (Schwarz iterations) the same equations each Schwarz iteration, but with different boundary conditions. More specifically, the chosen DDM is a Schwarz Waveform Relaxation method halpern3 (); GanderHalpernNataf (); halpern2 (); jsc (), which is an elaborated fixed point algorithm. Notice that SWR could as well be used as a preconditioning technique, but we would a priori then not benefit from the scaling effect, as using the latter approach we will still have to solve a huge linear system, although of course very preconditioned. This method is characterized by a choice of transmission conditions (or boundary conditions) on each subdomain, derived from the solution to local wave equations. This popular method is for instance analyzed for a 2-subdomain Schrödinger equation in halpern2 (); jsc (); lorin-TBS (); lorin-TBS2 (); AML ().
A SWR-Galerkin method for solving the -body Schrödinger equation is developed in this paper, using 2 different types of basis functions. The first basis is composed by Gaussian functions. The second basis which is used, is constituted by local Slater’s determinants in the FCI formalism. The corresponding SWR method requires i) the computation of 1-electron orbitals, that is eigenfunctions of local -body Schrödinger Hamiltonians, ii) from which we construct local Slater’s determinants. These Slater determinants can be used as local basis functions. The SWR methods developed in this paper are first applied to solve the stationary Schrödinger equation using the Normalized Gradient Flow (NFG) method bao (). The NGF is a minimization method, which consists in solving a normalized Schrödinger equation in imaginary time (this is why, it is also referred in the literature as the imaginary time method), or equivalently a normalized heat equation, with variable integration times. The SWR method is next applied to the time-dependent Schrödinger equation, modeling a molecule subject to an external electric field. For sufficiently intense fields, the -body wavefunction is expected to be delocalized, requiring then large -dimensional computational domains and then justifying the use of DDM. The purpose of this paper is not to show some high-dimensional simulations but rather to precisely describe aefficient general strategy for addressing the N-body problem. Some numerical results are however proposed in one dimension () for 2 electrons () to illustrate the proposed approach.

1.1 -body Time-Independent Schrödinger Equation

The stationary -particle Schrödinger equation, under the Born-Oppenheimer approximation, reads in dimensions CAM15-10 (); ostlund (), as follows


with , where is the spatial coordinates of the th electron, and its spin. In (1), the wavefunction is an eigenstate associated to the eigenvalue The Schrödinger Hamiltonian , for electrons and fixed nuclei (Born-Oppenheimer approximation) is given by

where denotes the position of the th nucleus, and its charge. In order to ensure the antisymmetry of the wavefunction , due to Pauli’s exclusion principle

for any , one can consider the traditional FCI approach ostlund () based on Slater’s determinants. Assume that is a set of orthonormal spatial orbitals in , and that , denote the spin coordinates, we can then define orbitals as follows: , , for . Next, -spin orthogonal orbitals are defined among Slater’s (antisymmetric) determinants.

Slater’s determinants , can then be used as basis functions in order to compute the eigenfunctions of . However, it is possible to construct more simple local basis functions . In particular, we will also use in this paper, Gaussian basis functions, in Section 2. Once, the basis functions are selected, any wavefunction is then expanded as . From now on and for the sake of the presentation, we will omit the spin variable and only consider the spatial coordinates.

1.2 -body Time-Dependent Schrödinger Equation

The -body time-dependent Schrödinger equation (TDSE), for in length gauge (LG) gauge () reads


where denotes a given external electric field, under the dipole approximation (wavelength of the electric field much larger than the spatial scale of the -body system). That is, for

which can be rewritten

with , and where

and . If the basis functions are Slater’s determinants, their orthogonality implies that is the identity matrix (when considering Neumann boundary conditions). When we are interested in the interaction of -electrons with an intense external field , it is necessary to include a very large number of 1-electron orbitals (, ), in particular “non-local” ones, as ionization is also expected PBC (); AB12 (). As a consequence, actual computations in all is, in principle, necessary.

1.3 Normalized Gradient Flow (NGF) method

In order to solve the time-independent Schrödinger equation, the method which is proposed is the imaginary time method, also referred in the Mathematics literature as a Normalized Gradient Flow (NGF) method. Let us rewrite the time-dependent Schrödinger in a compact form in :


The NGF method for computing the ground state of the Schrödinger Hamiltonian, consists of solving the time-dependent Schrödinger equation in imaginary time, and normalizing the solution at each time iteration. The converged state minimizes the energy functional

defined by


More specifically the ground state is constructed by solving for and ,


In the above system of equation, are discrete times, is an initial data for the time marching algorithm discretizing the projected gradient method and pointwise , see bao () for and .
However by construction, the computed state is not antisymmetric. As a consequence a constraint must be added in order to ensure that Pauli’s exclusion principle is well satisfied. We denote by an antisymmetrization operator. For instance, if , and we have

We assume that is symmetric, that is for any in

In the general situation, is antisymmetrized using odd permutations. We notice that satisfies the heat equation


as , are linear operators, and as . We have now to show that the following algorithm is energy decreasing for any and :


We notice first that:

Then following Theorem 2.1 in bao (), the energy defined in (4) satisfies

We then conclude that

Proposition 1.0

Assuming that is a symmetric potential, the algorithm (8) is convergent to an antisymmetry state of minimal energy.

the NGF algorithm will converge to the minimum energy antisymmetric state.

Remark 1.0

A more straightforward approach is simply to notice that if the initial data is antisymmetric (), then the solution to the heat equation will be antisymmetric as long as the potential is symmetric. This is a simple consequence of the uniqueness of the Cauchy problem associated to (6). Then, as mentioned in Theorem 2.2 from bao (),


where is defined as

1.4 Organization of the paper

This paper is organized as follows. In Section 2, we present the construction of Gaussian local basis functions. We then propose in Section 3, a methodology to construct local Slater’s determinants which can be used as local basis functions. Some properties of local Slater’s determinants, as well as the efficient construction of local Hamiltonians is discussed in this section as well as in A. Section 4 is devoted to the derivation and implementation of the Schwarz Waveform Relaxation algorithm for solving the -body Schrödinger equation. Some mathematical properties of the SWR will be recalled in this section, and their computational complexity will be discussed in B. Sections 5 and 6 are devoted to some numerical experiments for solving the time-independent and time-dependent -body Schrödinger equations, in one dimension. More specifically, the experiments are performed using local Gaussian basis functions in Section 5, and local Slater basis functions in Section 6. We finally conclude in Section 7.

2 Local Gaussian basis functions

The domain decomposition method for solving the -body Schrödinger equation which is proposed in this paper is based on a Galerkin method. The choice of the local basis functions is of crucial matter in order to make the computation as efficient as possible. Before considering complex basis functions in Section 3, we study the methodology with simple basis functions. A natural choice is to use Gaussian functions.
In order to simplify the notations, we will consider here, the case , . The extension of the following ideas is straightforward for arbitrary and and is shortly discussed at the end of this subsection. We denote by an infinite sequence of open intervals, such that: and , for , and , for any and in . Naturally we have . We denote by the set of one-dimensional Gaussian functions, defined by


where is a subdomain-dependent (-index) positive number for Electron (), and is a sequence of Gaussian centers. When the ’s are subdomain and particle independent, we will use the notation . Now, we can construct local basis functions for any . From any localized orbitals , , with in (basis function indices) and in (subdomain indices), we define by:

In term of support, we have

If and is taken subdomain independent, the local basis functions are actually identical in all the subdomains, which is quite convenient from a computational point of view, as we only need to construct once for all, a unique free-particle Hamiltonian. The weakness of this approach is that naturally, as the local basis functions do not contain any particular information, a large number should be used. In Fig. 1, we present in a given subdomain, the local Gaussian basis functions. The construction to Gaussian basis functions for particles in dimensions is naturally straightforward by considering the tensor products of local Gaussian functions: . The analysis of convergence of the Galerkin method applied to the Schrödinger equation, and using Gaussian basis functions was presented in GalerSchro ().

Figure 1: Gaussian basis functions in one subdomain for and .

In order to directly construct antisymmetric basis functions (at least locally) it is possible to construct (spinless) Slater-like Gaussian basis functions ostlund (), from any localized Gaussian functions , , with in and in ,


In the next section, we consider more elaborated antisymmetric basis functions using the traditional Slater’s determinants computed from 1-electron orbitals.

3 Local orbitals and Local Slater’s Determinants as basis functions

This section is devoted to the construction of Local Orbitals (LO’s) and Local Slater’s Determinants (LSD’s).
As we have done in the previous section, we will detail the case and , that is a two-body problem in one-dimension. This is motivated by the fact that the extension to the general case (arbitrarily and case) can be deduced from CAM15-09 (); CAM15-10 () and does not present any fundamental difficulty, but would complexify the notations. The material presented here will be used for the Schwarz Waveform Relaxation (SWR) Domain Decomposition Method (SWR-DDM) presented in Section 4. The local orbitals and Slater’s determinants will allow for the construction of local Hamiltonians and local wavefunctions, from which we will reconstruct a global wavefunction. The basic idea is to construct local, in the sense subdomain dependent, Slater’s determinants from local 1-electron orbitals. This procedure can be applied to any subdomain, or only in some of the subdomains, typically those containing the nuclei.

We denote by () the coordinate of the th particle. The Schrödinger Hamiltonian reads, for fixed nuclei

where denotes the position of the th nucleus and its charge. Antisymmetry of the wavefunction reads

3.1 Local FCI procedure

We denote by the set of 1-electron orbitals, which will allow for the construction of the compact support localized orbitals, LO’s, denoted by 111Top index refers to subdomain , and bottom index to full orbital index. Typically should satisfy, for any and


By construction, we will assume that . In order to solve the stationary Schrödinger equation, we choose the FCI model for a -electron problem. The latter is based on (spinless) Slater’s Determinant basis functions (SD’s) which are defined as follows. From any localized orbitals , , with in and in , Slater’s determinants as follows:


Notice that in practice the number of determinants to compute can be reduced. For instance, for only indices should be considered. In the following, we will denote by , the support of any function with respect to its -variables. As a consequence:


In other words, the support of any is compact and is strictly included in the union of with the subdomains having an edge in common with . By construct, is naturally antisymmetric.

3.2 Local orbital construction

The domain decomposition introduced above, allows for an adaptive selection of 1-electron orbitals per-subdomain. The key points are i) the number of nuclei, ii) their location, and iii) in the time-dependent case, the intensity of the external electric field. Notice that for any subdomain , we select 1-electron localized orbitals, . Then, from two sets of localized orbitals, , , we can construct LSD’s (12). From a practical point of view, we consider a finite number , of one-dimensional subdomains partially covering : . Notice that this will force us to impose absorbing conditions at the global computational domain boundary ABC (); MOLPHYS (). Then, for each subdomain , we will select LSD’s , among determinants. Notice however that the procedure which is presented below, may only be relevant for subdomains containing at least one nucleus. In the other subdomains, local Gaussian basis functions could be considered. The stationary wavefunction , solution to the Schrödinger equation, will then be searched in each , in the form

where are the unknown coefficients. We now detail the procedure to construct the localized orbitals under the condition (11), for . We consider as a generic example the case of the -molecule, corresponding to .

The approach which is proposed is based on ideas presented in CAM15-09 (). Rather than post-processing the full domain 1-electron orbitals, we directly construct the smooth localized orbitals with compact support, and with orthogonality properties. This is possible thanks to the use of i) infinite potentials at the subdomain boundary, and ii) of mollifiers when a subdomain contains a nucleus singularity. We proceed as follows. We consider the two following situations, for a given subdomain , with .

  • contains a nucleus singularity. Only a few subdomains belong to this first category, in particular when we are interested in the time-dependent Schrödinger equation for intense field-particle interaction. In that case, mollifiers will allow for an arbitrarily accurate smoothing of the nucleus singularities. Notice that in 1-d, the singularity treatment is different than in 3-d. Indeed in the latter case, we benefit from the fact that a Coulomb potential, up to a multiplicative constant is a fundamental solution to Poisson’s equation. This property allows for an accurate and efficient treatment of the localized orbitals. The Coulomb potential is then approximated by a smooth function , thanks to mollifiers as defined in CAM14-43 (), and such that:


    where , which also satisfies


    As a consequence, a smooth approximation of the Coulomb potential using (15), can be constructed with in . Notice that this property is also fundamental for efficiently computing the -dimensional integrals in order to construct the global discrete Hamiltonian CAM14-43 (). In 1-d, the fundamental solution of the Poisson equation is and the latter property does not occur anymore. Instead, we directly computed using (14) with defined by:


    where refers to the order of the mollifier and the scaling factors are explicitly defined in CAM14-43 (). For instance, for , we represent in Fig. 2 for a unique domain (Left) and (Right). In particular it is proven in CAM14-43 (), that for any smooth function

    for some positive sequence .

    Figure 2: (Left) Mollifiers (), (Right) and for .

    Once is computed, we introduce a barrier potential as in CAM15-09 ()


    where i) the smooth function is equal to for and for , ii) , and iii) and are imposed. The support of the localized orbitals is then , where denotes the coordinates of the center of the subdomain . We typically choose to ensure that the localized orbitals are not null at ’s boundary. A contrario, taking too large will lead to a loss of computational efficiency due to a large localized orbital support. In , we then solve the following one-dimensional one-electron eigenvalue problem

    where (resp. ) is the orbital (resp. subdomain) index and where . Notice that the choice of the localized orbitals is motivated by physical considerations. When we are interested in field-particle interaction, for subdomains containing the nuclei, we will select the localized orbitals corresponding to the lower energy states, as they will be predominant in the overall wavefunction in the vicinity of the nucleus singularities.

  • does not contain any nucleus singularity. In that case, the regularization of the Coulomb potential through mollifiers is naturally useless. The localized orbitals are then directly obtained by solving

    Similarly to the previous case (subdomain containing the nuclei), the selected localized orbitals will strongly depend on the relative position of the nuclei . Alternatively, for those subdomains, we can use local Gaussian basis functions as described in Section 2.

Once the localized orbitals are computed, we can construct the discrete Schrödinger Hamiltonian. The approach which is proposed below will benefit from i) the fact that the localized orbitals are selected accordingly to the position of the nuclei, ii) the compact support of the LO’s and iii) their orthogonality property (more specifically their extension by to all ). This last point necessitates some precisions. First, we notice that by construction for any , the supports of and of with are disjoint, so that these LO’s are trivially orthogonal. By construction, the LO’s of any are also orthogonal to each other. For or , the orthogonality of the LO’s and of is not, a priori, ensured. However, by construction for any and

Then, as the LO’s (smoothly) vanish at the boundary of their support, for small enough, we expect that to be small. For and (which is relatively very large) and () 5 subdomains, we represent for subdomain (resp. ) (resp. ), for in Figs. 3.

Figure 3: First 4 eigenstates LO’s (left), (right) for .

3.3 Augmented local bases

By construction the local Slater’s determinant basis functions are null at the boundary of the subdomains. This can constitute an issue if the overlapping zone between two subdomains is too narrow, as in such zones the basis functions are basically null or very small, see Fig. 4 (Left). In order to fix this issue, a simple solution consists of adding Gaussian basis functions all around the subdomains Fig. 4 (Right). It will then ensure that in any overlapping zone the local wavefunctions could be properly transmitted from one subdomain to another thanks to the transmission conditions.

Figure 4: (Left) Local basis functions overlapping issue. (Right) Additional Gaussian basis functions ensuring a proper transmission.

Notice that in A, we present a general strategy to efficiently compute the integrals involved in the construction of the local Hamiltonians, using the formalism proposed in CAM15-10 ().

3.4 Important remarks about subdomain and local basis function indices

In order to lighten the presentation, some compact notations will be used along the paper.

  • Functions will systematically refer to basis functions for a unique domain problem, that is without DDM.

  • In a given subdomain (-index), the local basis functions (-index) can be also denoted by . For a given two-dimendional subdomain , the local basis functions could be also denoted by , where denotes the basis function index, and the one-dimensional subdomain indices. This notation was already used in Subsection 3.2 to define local Slater’s determinants.

  • For and , , and , will also be denoted by . In this case, the subdomains will be reindexed as .

In general, for the local basis functions the top index always refers to the basis function index, and the bottom one to the subdomain index.

4 Schwarz waveform relaxation domain decomposition method for the Schrödinger equation

We first decompose , in overlapping hypercubes where is an integer parameter, and apply a Schwarz waveform relaxation algorithm lorin-TBS (); lorin-TBS2 (). We present two different approaches. The first one leads to an a posteriori antisymmetric wavefunction, and second one ensures a priori Pauli’s exclusion principle (see A). In the following, we denote i) the artificial interfaces by , for any , and ii) by the overlapping regions , for any . For convenience, we also denote the Cartesian product , such that for all . We now denote by the solution to the -body TDSE in , at time and Schwarz iteration . For any , , we construct a basis of local basis functions (Gaussian functions or Slater’s determinants) in , denoted by in order to compute . Basically, we will solve local time-dependent or time-independent local Schrödinger equations and reconstruct a global solution to the global Schrödinger equation. We then never compute the global solution from a global discrete Hamiltonian, but rather by computing local wavefunctions (one per subdomain) using discrete local Hamiltonians, Figs. 5. SWR algorithms are in particular, studied in halpern3 (); GanderHalpernNataf () and allow for a consistent decoupling on smaller subproblems of high dimensional (non-local) classical, quantum and relativistic wave equations.

Figure 5: (Left) Domain decomposition: Overlapping subdomains are represented in red. Blue subdomains do not overlap. (Right) Domain decomposition with overlapping region on with in

4.1 Schwarz Waveform Relaxation algorithm for the TDSE

We detail the DDM algorithm first for 2 subdomains , with , then for zones where more than subdomains overlap.

Two-subdomain overlapping zones. Assume first that is a given function. The Schwarz Waveform Relaxation algorithm (SWR) with and , reads in LG and for subdomains, as


where is a transmission operator defined at .

Multi-subdomain overlapping zones. The proposed decomposition requires a special treatment in zones, generically denoted see Fig. 6 (Left), where more than 1 subdomain overlap with . We denote by the set of indices of the subdomains, distinct from , sharing the zone with . Notice that for interior subdomains (that is excluding the subdomains of the external layer) . The approach which is proposed is actually an averaging process. Let us generically denote by the interface of involved in the transmission conditions. The condition we impose at , thanks to the operator , is defined by:

In order to clarify the process, let us detail the case and , with a total of subdomains. At a given interface of , located at the righttop of a given subdomain with , we assume that there are subdomains involved in the transmission condition, namely , see Fig. 6. Then we impose at :

Figure 6: (Left) Domain decomposition with and : treatment of transmission conditions on interfaces belonging to more than subdomains. The black segment belonging to . (Right) Domain decomposition with and : smooth subdomain boundary

Notice also that a special treatment of the TC at the cross-points can improve the convergence deterioration at these locations. We do not explore this issue, but we refer to the related literature CP1 (); CP2 (). A simple way to circumvent this difficulty consists in regularizing the corners of the subdomains, as shown in Fig. 6 (Right). In a Galerkin-method framework, managing such smooth regions is then straightforward. Another simple approach is presented in stcyr (), allows to avoid the discretization of the right-hand-side of the transmission conditions.

Selection of the transmission conditions. From the convergence and computational complexity points of view, the selection of the transmission operator , is of crucial matter lorin-TBS (), halpern2 (). The usual types of transmission conditions (TC) are now shortly recalled. The most simple approach is a Dirichlet-based TC’s, where is an identity operator. That is, we impose:

In the literature, this is referred as the Classical Schwarz Waveform Relaxation (CSWR) algorithm. CSWR is in general convergent, and trivial to implement, but unfortunately exhibits usually very slow convergence rate and possibly stability issues at the discrete level halpern2 (); halpern3 (); jsc (). The CSWR method also necessitates an overlapping between the subdomains. As a consequence, more appropriate TC’s should be derived, such as Robin-based TC’s. Denoting , the outward normal vector to , and for in imaginary time (heat equation halpern3 ()), and (real time, Schrödinger equation halpern2 ()), the Robin TC’s read:


In the numerical simulations, will be taken interface-independent, and will simply be denoted by . This method will be referred as the Robin-SWR method. Along this paper, the Robin-SWR we will be used, as it is known for allowing for a good compromise between convergence rate and computational complexity, see halpern2 (). Robin-SWR can be seen as an approximation of Optimal SWR (OSWR) which are based on transparent or high order absorbing TC’s and reads, at