Domain Decomposition Method for the Nbody TimeIndependent and TimeDependent Schrödinger Equation
Abstract
This paper is devoted to the derivation of a pleasingly parallel Galerkin method for the timeindependent body Schrödinger equation, and its timedependent 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.
keywords:
body Schrödinger equation, domain decomposition method, mollifiers, parallel computing.Contents
 1 Introduction
 2 Local Gaussian basis functions
 3 Local orbitals and Local Slater’s Determinants as basis functions
 4 Schwarz waveform relaxation domain decomposition method for the Schrödinger equation
 5 Numerical experiments: Gaussian local basis functions
 6 Numerical experiments: Local Slater’s Determinants
 7 Conclusion
 A Antisymmetric wavefunction reconstruction
 B Computational complexity analysis of the SWRDDM solver for body equation
1 Introduction
This paper is devoted to the derivation of a pleasingly parallel realspace algorithm for solving the body Schrödinger equation. It is wellknown 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 “worstcase scenario”, would be a finite difference approximation. However, elaborated alternatives exist such as the Full Configuration Interaction (FCI) ostlund (), which requires highly nontrivial basis functions, and which then allows for the construction of relatively compact discrete Hamiltonians. The timeindependent and dependent Schrödinger equations are computed using algorithms (linear system and eigenvalue solvers) basically requiring many highdimensional matrixvector 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 socalled transmission conditions DoleanBook (). Secondly, we can benefit from a scaling effect, as the computational complexity of realspace numerical (timeindependent 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 2subdomain Schrödinger equation in halpern2 (); jsc (); lorinTBS (); lorinTBS2 (); AML ().
A SWRGalerkin 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 1electron 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 timedependent 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 highdimensional simulations but rather to precisely describe aefficient general strategy for addressing the Nbody problem. Some numerical results are however proposed in one dimension () for 2 electrons () to illustrate the proposed approach.
1.1 body TimeIndependent Schrödinger Equation
The stationary particle Schrödinger equation, under the BornOppenheimer approximation, reads in dimensions CAM1510 (); ostlund (), as follows
(1) 
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 (BornOppenheimer 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 TimeDependent Schrödinger Equation
The body timedependent Schrödinger equation (TDSE), for in length gauge (LG) gauge () reads
(2) 
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 1electron orbitals (, ), in particular “nonlocal” 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 timeindependent 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 timedependent Schrödinger in a compact form in :
(3) 
The NGF method for computing the ground state of the Schrödinger Hamiltonian, consists of solving the timedependent Schrödinger equation in imaginary time, and normalizing the solution at each time iteration. The converged state minimizes the energy functional
defined by
(4) 
More specifically the ground state is constructed by solving for and ,
(5) 
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
(6) 
as , are linear operators, and as . We have now to show that the following algorithm is energy decreasing for any and :
(7) 
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 (),
(8) 
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 timeindependent and timedependent 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 onedimensional Gaussian functions, defined by
(9) 
where is a subdomaindependent (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 freeparticle 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 ().
In order to directly construct antisymmetric basis functions (at least locally) it is possible to construct (spinless) Slaterlike Gaussian basis functions ostlund (), from any localized Gaussian functions , , with in and in ,
(10) 
In the next section, we consider more elaborated antisymmetric basis functions using the traditional Slater’s determinants computed from 1electron 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 twobody problem in onedimension. This is motivated by the fact that the extension to the general case (arbitrarily and case) can be deduced from CAM1509 (); CAM1510 () 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 (SWRDDM) 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 1electron 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 1electron orbitals, which will allow for the construction of the compact support localized orbitals, LO’s, denoted by ^{1}^{1}1Top index refers to subdomain , and bottom index to full orbital index. Typically should satisfy, for any and
(11) 
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:
(12) 
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:
(13) 
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 1electron orbitals persubdomain. The key points are i) the number of nuclei, ii) their location, and iii) in the timedependent case, the intensity of the external electric field. Notice that for any subdomain , we select 1electron 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 onedimensional 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 CAM1509 (). Rather than postprocessing the full domain 1electron 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 timedependent Schrödinger equation for intense fieldparticle interaction. In that case, mollifiers will allow for an arbitrarily accurate smoothing of the nucleus singularities. Notice that in 1d, the singularity treatment is different than in 3d. 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 CAM1443 (), and such that:
(14) where , which also satisfies
(15) 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 CAM1443 (). In 1d, 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:
(16) where refers to the order of the mollifier and the scaling factors are explicitly defined in CAM1443 (). For instance, for , we represent in Fig. 2 for a unique domain (Left) and (Right). In particular it is proven in CAM1443 (), 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 CAM1509 ()
(17) 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 onedimensional oneelectron 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 fieldparticle 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.
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.
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 CAM1510 ().
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 twodimendional subdomain , the local basis functions could be also denoted by , where denotes the basis function index, and the onedimensional 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 lorinTBS (); lorinTBS2 (). 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 timedependent or timeindependent 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 (nonlocal) classical, quantum and relativistic wave equations.
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.
Twosubdomain overlapping zones. Assume first that is a given function. The Schwarz Waveform Relaxation algorithm (SWR) with and , reads in LG and for subdomains, as
(18) 
where is a transmission operator defined at .
Multisubdomain 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 :
Notice also that a special treatment of the TC at the crosspoints 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 Galerkinmethod framework, managing such smooth regions is then straightforward. Another simple approach is presented in stcyr (), allows to avoid the discretization of the righthandside 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 lorinTBS (), halpern2 (). The usual types of transmission conditions (TC) are now shortly recalled. The most simple approach is a Dirichletbased 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 Robinbased 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:
(19) 
In the numerical simulations, will be taken interfaceindependent, and will simply be denoted by . This method will be referred as the RobinSWR method. Along this paper, the RobinSWR we will be used, as it is known for allowing for a good compromise between convergence rate and computational complexity, see halpern2 (). RobinSWR can be seen as an approximation of Optimal SWR (OSWR) which are based on transparent or high order absorbing TC’s and reads, at