{'fraction_non_alphanumeric': 0.045647507751013595, 'fraction_numerical': 0.023467684235630813, 'mean_word_length': 4.438651102464332, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
HeTDSE: GPU based program to solve full-dimensional-time-dependent Schr\"{o}dinger equation for two-electron Helium subjected to strong laser fields

We present a GPU based openACC fortran program named HeTDSE, which provides an efficient way to investigate the non-perturbative electronic dynamics of helium subjected to a strong laser pulse by solving full-dimensional two-electron time dependent Schr\"{o}dinger equation (TDSE). OpenACC is a directive-based programming model for accelerators without modifying the underlying CPU code itself. HeTDSE uses B-Spline basis sets expansion method to construct the radial part of the wavefunction, and the spherical harmonic functions is used to express for the angular part. Adams algorithm is employed for the time propagation. Our example shows HeTDSE running on an NVIDIA Kepler K20 GPU can outperform the one on an Intel E5-2640 single CPU core by a factor of 147.


I. INTRODUCTION
The rapid development in laser technologies opens a way for scientists to probe and even control the fundamental dynamics of electron correlations [1][2][3][4][5][6][7][8]. As the simplest multi-electron atom, Helium is an idea starting point for exploring electron correlation dynamics in multi-electron systerms [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. However, due to its six degrees of freedom, the response of Helium to strong fields is considerably more complicated than that of single-electron atoms, which poses great theoretical and computational challenge. To overcome this difficulty, a set of conventional CPU parallel computing techniques have been developed to numerically solve the time-dependent Schrödinger equation (TDSE) of helium subjected to a laser pulse: Peng et al. investigated the electron correlation effects in two-photondouble-ionization (TPDI) of helium by finite-elementdiscrete-variable-representation (FEDVR) method [24]; Parker et al. used the finite-difference method to calculate the above threshold ionization (ATI) process [25]; Piraux et al. investigated the electron correlation effect by using the Gauss-Sturman function [14]. All of these simulations are very likely a numerical virtual experiment on the servers (about 200 to 1000 CPU cores are used in [24,25]). Thus, more efficient algorithm is needed to further promoting the numerical simulations of Helium (and even more complex multi-electron systems beyond Helium) in strong laser fields investigations. One of this solution is Graphic Processing Unit(GPU) programming, which is the most powerful high performance computing tool so far and widely applied both in science and engineering numerically studies [27][28][29][30][31][32][33][34].
GPU contains hundreds of computing cores and is originally designed for the highly parallel process of graphic * zhaoxi719@ksu.edu rendering [27,28]. Comparing with CPU, the computing performance of GPU can be increased by tens of times with proper optimizing [27][28][29][30][31][32][33][34]. To do this, a so-called "Compute Unified Device Architecture (CUDA)" GPU programming model is present [28]. However, porting of legacy CPU-based codes with CUDA often necessitates explicit compute and data management, thus requiring significant structural changes to existing applications [29]. Therefore, we make the choice to use Ope-nACC, which gives an alternative model as a GPU programming scheme [35]. OpenACC is a set of directivebased extensions to C, C++ and Fortran that allow programmers to annotate regions of code and data for offloading from a CPU host to an attached GPU, without requiring modification to the underlying CPU code itself. Programmers simply insert OpenACC directives before specific code sections to engage the GPU to accelerate the code. This approach enables the compiler to target and optimize parallelism automatically. More example programs and detailed description of openACC can be found in openACC official website [35].
In this work, we present a GPU based openACC fortran program HeTDSE, which is an efficient tool to investigate the non-perturbative electronic dynamics of helium subjected to a strong laser pulse by solving full-dimensional two-electron TDSE. It goes beyond the single-electron-approximation (SEA) approach and includes the response to the field of all two electrons. To build Helium wavefunction, B-Spline basis sets, which were widely used in computational atomic and molecular physics [36][37][38][39][40][41][42][43][44][45][46][47], are used to construct the radial part of the wavefunction, while the spherical harmonic functions is used to express for the angular part. The reason we use B-Spline basis sets to expand Helium radial wavefunction is that, B-Spline function has great advantages of describing both bound and continue states with small number basis sets [36,37,39]. Adams algorithm is employed for the time propagation [48]. Another advantage of using B-Spline basis sets and Adams method is that it is easy to parallelize the code and will get an excellent paralleling scaling with openACC.
The paper is organized as follows: In section 2, we present the theoretical background on HeTDSE. In Section 3 we exhibit an overview of the package structure, the input, output files and the code parallelizing. In Section 4 we show several test applications of HeTDSE. The parallel efficiency is given in Section 5. We present our conclusions in Section 6. Atomic units are used throughout, unless otherwise stated.

II. MATHEMATICAL SETUP AND ALGORITHM
A. B-Spline function B-splines are functions designed to generalize polynomials for the purpose of approximating arbitrary functions, we use B-Spline basis sets to construct the Helium wavefunction in HeTDSE. Thus, we begin this Section with a brief description of B-Spline function, very details of B-Spline function can be found in ref [39].
A B-Spline function is defined by the order k and a set of the breakpoints {t i }, The radial coordinates of Helium wavefunction is discrete in the interval [r min , r max ] by the breakpoints {t j }. Three types of breakpoints sequence {t j } are used in HeTDSE [39]. The first one is linear sequence The sequence is divided into a equal length. Second, exponential sequence, This sequence tends to a linear sequence as γ → 0 while all points exponentially accumulate close to r min as γ → ∞.
Third, Linear-parabolic sequence. A useful sequence adapted to a good description of both the bound and the continuum states associate a linear spreading at large distances with a quadratic sequence close to the origin [39].
Other types sequences can be found in [39] and its references.
The derivative of a B-spline of order k can be expressed as a linear combination of B-Splines of the same order The B-Spline function order k = 7 is used throughout this work, so we don't write out k in B-Spline functions for simplification.

B. Time independent Schrödinger equation
The Helium eigenstate ϕ n (r 1 , r 2 ) and its corresponding eigenvalue E n is solution of time independent Schrödinger equation (TISE) of Helium: where H 0 (r 1 , r 2 ) is the laser-free Hamiltonian where −∇ 2 i 2 − 2/r i with i = 1, 2 denotes the kinetic energy and potential energy of electron i, and 1 |r1−r2| is the electron-electron correlation term. The time-independent Helium wavefunction is expanded as where N is the number of the B-Spline basis sets for each time independent wavefunction ϕ n (r 1 , r 2 ), l 1 (l 2 ) denotes the angular momenta for electron 1(2), L is the total orbital angular momentum, M is its z-component, S is the total spin, {c i } is the expansion coefficient and each i corresponds to a set of {n 1 , l 1 , n 2 , l 2 }. Coupled spherical harmonic functions is used to express for the angular part of the time-independent wavefunction: where m 1 (m 2 ) is the z-component of l 1 (l 2 ) and Y lm (r) is usual spherical harmonic functions.
The wavefunction ϕ n (r 1 , r 2 ), eigenvalue E n as well as the expansion coefficient {c i } in (8) can be got by directly diagonalizing eq. 6. For this purpose, there are a set of matrix integrals would be done: the kinetic energy integral matrix element: the potential energy integral matrix element: the overlap integral: and the coupling term integral: K ij , P ij and O ij are straightforward to discretize and calculated. However, the calculation of C ij is different. To calculate C ij we expand the electron-electron correlation term a truncated multipole series: Thus, each of the terms are handled in a similar manner to the one-electron operators [37]. In HeTDSE, all the integrals are carried by Gauss-Lagrange integration method, which has been widely used in other works [26].

C. Time dependent Schrödinger equation
We solve the Helium TDSE within dipole approximation and length gauge. The full-dimension-TDSE of helium is where the total Hamiltonian of atom-laser system is The interaction Hamiltonian of the helium and laser pulse is where the vector potential of the laser pulse is hereε is the laser polarization direction, ω and φ are the frequency and the carrier envelope phase, respectively. f (t) is the temporal envelope.
The total time dependent wavefunction Ψ(r 1 , r 2 , t) can be expanded in terms of the field-free atomic eigenfunctions Here N total is the total number of the wavefunction ϕ n (r 1 , r 2 ). Substituting equation 19 into 15, the TDSE is converted into a set of coupled partial differential equations, which can be solved by the Adams method, details of this algorithm can be found in [48]. The energy difference E mn = E m − E n and the transition dipole element between ϕ n (r 1 , r 2 )| and |ϕ m (r 1 , r 2 ) , can be calculated from the solution of eq. 6.

D. Absorbing boundary
The absorbing layer is defined by using the radius R = r 2 1 + r 2 1 . Let r mask be the radius where the absorbing layer starts. An absorbing layer between r mask and r max is used to smoothly bring down the wavefunction and to prevent the unphysical reflection from the boundary. The absorbing function has the following form: Although we use an absorbing to avoid the nonphysical reflecting, the simulation box still need to be set big enough so that the physical system is not perturbed by the absorbing boundaries. Thus, there are two interpretation of the ionization yield. First, we directly calculate the single (double) ionization yield by summing all the possibilities of the wavefunction with an eigenvalue larger than −2.0(0.0): Alternative, the ionization probability is calculated as

III. DESCRIPTION OF THE PACKAGE
HeTDSE code package contains 9 fortran files and 4 input files. The fortran driver programs, functions, subroutines, input and output files are all introduced briefly in this section.

A. fortran program files
These fortran codes should be run one by one: Firstly, running eigen-equation.f to solve Eq. 6 to get wavefunction ϕ n (r 1 , r 2 ), eigenvalue E n as well as the expansion coefficient {c i }; Then, with the output files of eigen-equation.f, running dipole.f90 to get the transition dipoles element ϕ n (r 1 , r 2 )| (r 1 Y 10 (r 1 ) + r 2 Y 10 (r 2 ))| ϕ m (r 1 , r 2 ) ; Next, running matrix.f90 to prepare input files for tdse.f90; Finally, running tdse.f90 to solve the Eq. (20) to get the time-dependent wavefunction Ψ(r 1 , r 2 , t).
There are other five fortran programs order.f, rsg.f, wig.f, angl16.f90 and SUBROUTINE.f90 in HeT-DSE. These five programs are the "support codes", we DO NOT suggest the users to modulate them.

B. lower-level functions and subroutines
The lower-level functions and subroutines in this program are: PREQUAN : Get the index of the one electron functions B n (r) Y lm from 1 to n × (l max + 1).
gauleg: Calculate the Gauss-Lagrange integration. SingleInteg2012 : Calculate the integration P ij , K ij or O ij .
DmultiInteg2012 : Calculate the electron-electron integration C ij .
HAMILTON2012 : Construct the Hamiltonian. RSG: Diagonalize the matrix, get the energy level and wavefunctions.
ANG: Calculate the angle part of transition dipole element.
abr : Set the absorbing boundary.

C. Input files
There are totally four input files in HeTDSE, eigenequation.input, dipole.input, matrix.input and tdse.input, which contain input parameters used by eigneequation.f, dipole.f90, matrix.f90 and tdse.f90, respectively. In this subsection, we will present how to set these parameters in these input files one by one.
eigen-equation.input. Line 1: Set total angular momentum L in Eq. 8. Line 2: Set total spin in Eq. 8. Line 3: Set max angular momentum for each electron l max .
Line 4: Set the order of the B-Spline function k. Line 5: Set number of B-SPLINE function breakpoints, that is to say, n in Eq. 2, 3 and 4.
Line 6: Set total number of the basis sets for the He wavefunction N in Eq. 8.
Line 7: i 0 in Eq. (4). Line 8: Set the simulation box size in radial direction, r max .
dipole.input Line 1: Set total angular momentum L of ϕ n (r 1 , r 2 )| and |ϕ m (r 1 , r 2 ) in Eq. 8, respectively. Line 2: Set total spin M of ϕ n (r 1 , r 2 )| and |ϕ m (r 1 , r 2 ) in Eq. 8, respectively. Output files of dipole.f90: There is one output file after running dipole.f90: the transition dipole moment elements between a pair of states with neighbouring total angular momentums L and L + 1.
Output files of matrix.f90: There are three output files after running matrix.f90 : 1. eigenval.out: This file stores all the eigenvalues in the order L = 0, 1, 2, 3 4. c.dat: This file stores the solution a n (t) (n = 1, ..., N total ) of the coupled partial differential equations Eq. 20 at each time step. In principle, if we get a n (t), all the physical information can be retrieved.

E. OpenACC parallelizing implementation
In this subsection, we will explain the details of Ope-nACC implementation in HeTDSE. In HeTDSE package, more than 99% computation time would be paid to calculate transition dipole moment d mn and time propagation. Thus, we will focus on accelerating the two calculations with openACC. In dipole.f90 programs, four do-loops are needed to get all the transition dipole moment elements d mn (see Eq. 21 for mathematical expression), and each loop would run 2000-10000 times. The code and corresponding openACC accelerated implementation is shown below: although the computing scale is large, the code structure itself is really simple (nothing but a sum calculation) and openACC can achieve a high performance parallelizing scaling. All four do-loops are parallelized directly by inserting the OpenACC directive "!$ACC KERNELS", then the data transfer between the host and the GPU memory is automatically executed. The calculation in the area (line 2 to line 12) of directive is executed and accelerated on GPU. Optimization for time propagation is similar, here is the Adams time propagation code at each time step and corresponding openACC accelerated implementation:

! $ACC END PARALLEL LOOP
As we known, the data transfer between the host and GPU memory affects the computational time. To further minimize the cost of data transfer, we use the Ope-nACC "DATA COPY" before time propagating starts, and "END DATA" is used to release the GPU memory at the end of time propagating:

! $ACC END DATA
The datacopy (namelist) is the directive that copies the data from host to GPU memory, then data on GPU memory is used without the data transfer back and forth between the host and GPU every time step. Clearly, by inserting dotacopy, the calculational time is much saved.
The advantages of using B-Spline basis set and Adams method in HeTDSE are emphasized again at the end of this section: First, it is convenient to implement ope-nACC. Second, the parallelizing scaling has a high performance, which will be shown in section 6.

IV. SAMPLE RESULTS
In order to verify the accuracy of our program, we compare our results with previous literatures. In the calculations, the radius of the cavity is r 1max = r 2max = 70 and it is described by 30 B-Spline functions of order 7. We use l i=1,2 = 0, 1, 2, 3, 4 and L = 0, 1, 2, 3, 4 in below simulation examples. Linear-parabolic breakpoints sequence is chosen. Total number of basis sets N total = 11000 are used during the time dependent simulations.

A. ground and bound states calculation of Helium
In table 1, we show the eigenvalues of few low bound states for different total angular momenta. Table 1 shows that, for all the calculated levels, at least the accuracy up to two digits after the decimal point has been obtained. The density distributions of different atomic levels in coordinate space are shown in Fig. 1, where the results agree well with that in [23,43,44].

B. excited states dynamics
Now we turn to the second example: excited states dynamics. Here we focus on the carrier-envelope phase (CEP) effect on band-band state transition induced by a laser pulse. The CEP is a crucial parameter in describing the characteristics of a laser pulse, we can control the dynamic process of matter-laser interaction by measuring or adjusting the CEP [49][50][51][52][53][54]. Especially, the CEP effect on the bound-bound transition of an atom has been investigated both in theoretically and experimentally [49][50][51][52][53][54].
Here we try to reproduce the result from [52]. In [52], the authors used Hylleraas coordinates to reconstruct the wavefunction of Helium, and they introduced a parameter M to quantify the CEP effect: where P (φ max ) and P (φ min ) are, respectively, the maximum and minimum populations for a given excited state.
A large value of M corresponding to a strong CEP effect. In this simulation, the laser parameters are the same as [52]. We use HeTDSE get the value M for 1 D state after the laser ends, which is shown in Fig 2. Our result matches well with that in Fig.1(b) from [52].

C. excitation and ionization yields
Next, we calculate the excitation and ionization yields of helium in a strong laser pulse. Our basis sets covers the energy range located beyond the double-ionization threshold. The initial state is the ground state of Helium 1S 2 . The laser pulse has a duration of 3.8 fs and the peak intensity of 2.97 × 10 14 W/cm 2 , which is the same as those in [38,44]. Eq. 23. The present results in Fig. 3 are accordant with the data from Hasbani [38] and Scrinzi [44].

D. electrons wavepacket dynamics
We then calculate the electrons wavepacket dynamics both in coordinate space after the laser pulse. The frequency and FWHM of the laser pulse are ω = 1.0 and 2 optical cycles, respectively, with the intensity of 1.0 × 10 13 W/cm 2 . The continue state is collected as we get the wave function of helium at the end of the laser pulse. Using this wave function, the density distribution of the electrons in coordinate space can be obtained by ρ(r 1 , r 2 , t) = Ψ(r 1 , r 2 , t)Ψ * (r 1 , r 2 , t)r 2 1 r 2 2 dΩ 1 dΩ 2 .
(26) The density distribution of continuous states in coordinate space is shown in Fig. 4 when the laser is end. There is an evident single ionization characteristics in the figures, where the wave packet moves outside along the coordinate axis with time.

V. PARALLEL SCALING
To test the parallel efficiency of HeTDSE, we compare the serial CPU program (runs at Intel xeon E5-2640 CPU with 2.5GHz clock speed and 15MB L3 cache) and parallel GPU implementations (runs at NVIDIA K20 GPU with 2493 cores). The speedup factor for four simulation cases shown in Table II. The larger the basis number, the larger computation cost needs. All the simulations are done with PGI fortran compiler, the laser is 3.8 fs (which contains about 4000 time steps) and the simulation box is R max = 70. A speed up of 147 is achieved if 4300 basis sets are used. It indicates that as the simulation system size increases, this improvement becomes more and more pronounced.

VI. CONCLUSION
In this work, we present a program which solves the full-dimension-TDSE of helium using OpenACC+GPU simulation acceleration environment. We introduce how to convert the full-dimension-TDSE into coupled partial differential equations. These partial differential equations are solved by Adams method. Our program has two advantages: first, the codes are easily parallelled by adding few detectives and have a speed up of 147 on GPU, HeTDSE dose not have to use a super computer or a computer cluster, even a desktop computer with an openACC-enable GPU can run HeTDSE efficiently; second, we can transplant our program to other accelerators without rewriting the codes. By comparing with the literatures of the excited state dynamics and ionization yield of helium, the accuracy of our program has been verified. Our codes can be used to investigate the nonperturbative electronic dynamics of helium subjected to a strong laser pulse. Besides, for the programming for accelerators such as CUDA is difficult, we hope that HeT-DSE to be an example to help more researchers to handle the GPU calculation more easily by using OpenACC.