Parallel Implementation of Numerical Solution of Few-Body Problem Using Feynman’s Continual Integrals

. Modern parallel computing algorithm has been applied to the solution of the few-body problem. The approach is based on Feynman’s continual integrals method implemented in C ++ programming language using NVIDIA CUDA technology. A wide range of 3-body and 4-body bound systems has been considered including nuclei described as consisting of protons and neutrons (e.g., 3 , 4 He) and nuclei described as consisting of clusters and nucleons (e.g., 6 He). The correctness of the results was checked by the comparison with the exactly solvable 4-body oscillatory system and experimental data.


Introduction
The wave functions of the ground states of 3,4,6 He nuclei were calculated in Refs.[1,2] by Feynman's continual integrals method in Euclidean time [3,4].The present work is devoted to studying hydrogen, helium, lithium, and beryllium nuclides using the same approach.The few-body nuclei 2,3 H, 3,4 He were considered as consisting of protons and neutrons, whereas the nuclei 6 He, 6 Li, 9 Be were considered as α-cluster nuclei.The algorithm allowing us to perform calculations directly on GPU was developed and implemented in C++ programming language.The energy and the square modulus of the wave function of the ground states for several few-body nuclei have been calculated using NVIDIA CUDA technology.The results show that the use of GPU is very effective for these calculations.

Theory and computing
The energy E 0 and the square modulus of the wave function |Ψ 0 | 2 of the ground state of a system of few particles with coordinates q may be calculated by Feynman's continual integrals method using the propagator K E (q, τ; q, 0) in Euclidian time τ [4] e-mail: anaumenko@jinr.rue-mail: samarin@jinr.ru Here g(E) is the density of states with the continuous spectrum E ≥ E cont .For the system with a discrete spectrum and finite motion of particles the square modulus of the wave function of the ground state may be found in the limit τ → ∞ together with the energy E 0 The theoretical approach is described in detail in Ref. [2].The calculation of K E (q, τ; q, 0) for the fixed τ was performed by parallel calculation of exponentials for every random trajectory q k = f (q, k∆τ), where ∆τ is time step, N = τ/∆τ.The same effective pairwise nucleon-nucleon, nucleon-cluster and cluster-cluster interaction potentials V (r) were used for all the studied nuclei.The Monte Carlo algorithm for numerical calculations was developed and implemented in C++ programming language using NVIDIA CUDA technology.Calculations were performed on the NVIDIA Tesla K40 accelerator installed within the heterogeneous cluster [5] of the Laboratory of Information Technologies, Joint Institute for Nuclear Research, Dubna.The code was compiled with NVIDIA CUDA version 7.5 for architecture version 3.5.
The code implementing Feynman's continual integrals method was initially written for CPU.The calculation time of the ground state energy for 3 He with statistics 10 6 using Intel Core i5 3470 (double precision, 1 thread) and NVIDIA Tesla K40 (single precision) was ∼18377 sec and ∼47 sec, respectively, which corresponds to the impressive performance gain of ∼400 times [2].The typical calculation time of the square modulus of the wave function for the ground state of 3 He with statistics 10 6 on the mesh of 60•60•12 points using NVIDIA Tesla K40 was ∼10 hours, which yields an estimation of the calculation time using Intel Core i5 3470 (double precision, 1 thread) ∼170 days.This performance gain allows us to increase the statistics and the accuracy of calculations, reduce the mesh step, and may even enable calculations impossible before.
To check the correctness of the calculation of the propagator the comparison with the exactly solvable 4-body oscillatory system has been performed.For four particles with masses m 1 = m 2 = m 3 = m 4 = m interacting with each other by oscillator potentials the kinetic energy T and the potential energy V may be written using the normalized Jacoby coordinates x, y, z and the energy E 0 of the ground state is given by Assuming ω = 1, = 1, we obtain E 0 = −U 0 + 9.For example, in the cases of U 0 = 0, E 0 = 9 and U 0 = 15, E 0 = −15+9 = −6, whereas Feynman's continual integrals method with statistics N = 7•10 7 yields E 0 = 9.05 ± 0.1 and E 0 = −5.98 ± 0.02, respectively, which is close to the exact result.

Results
In the calculation of the propagator for the nuclei 2,3 H, 3,4 He nucleon-nucleon two-body effective strong interaction potentials V i− j (r) (i, j = n, p) with repulsive cores, similar to the M3Y potential [6], have been used (Figure 1a) Figure 1.a) Two-body interaction potentials V i− j (r) for the following pairs of particles: neutron-proton with parallel (solid line) and antiparallel spins (dotted line), neutron-neutron (dashed line), α-cluster-neutron (dashdotted line) and α-cluster-α-cluster (dash-dot-dotted line).b) The dependence of the logarithm of the dimensionless propagator KE on the dimensionless Euclidean time τ for 2 H (filled triangles), 3 H (circles), 3 He (dots), 4 He (empty triangles), 6 He (filled squares), 6 Li (empty squares), 9 Be (stars); lines are the results of linear fitting of the data; b 0 is a coefficient [2].
The calculated logarithm of the propagator for nuclei 2,3 H and 3,4 He is shown in Figure 1b.The theoretical binding energies E B = −E 0 obtained using formula (2) are listed in Table 1 together with the experimental values taken from the NRV knowledge base [7].It is clear that the theoretical values are close enough to the experimental ones.The observed difference between the calculated binding energies of 3 H and 3 He is also in agreement with the experimental values.
Momentum distributions after breakup in nuclear reactions show that 6 He nucleus consists of an α-cluster core and a two-neutron cluster (e.g., [8]).Similarly, the 6 Li nucleus consists of an αcluster core and a deuteron cluster.The 9 Be nucleus consists of two α-clusters and a neutron.There are several different approaches to the approximate analysis and solution of the three-body problem (e.g., [9]).Feynman's continual integrals method provides a simpler possibility of calculating the energy and the probability density for the ground state of the few-body system, because it does not require expansion of the wave function in a system of functions.The calculated dependence of logarithm of the propagator on the Euclidean time for nuclei 6 He (α + n + n), 6 Li (α + n + p) and 9 Be (α + n + α) is shown in Figure 1b.The α-cluster-nucleon and α-cluster-α-cluster strong interaction potentials V i− j (r) (i, j = n, p, α) were used in the form of the combination of Woods-Saxon potentials (Figure 1a) where s = 2, 3.The values of parameters are given in Ref. [10].The obtained theoretical energies of separation into cluster(s) and nucleon(s) E S = −E 0 are listed in Table 1 together with the experimental values taken from the NRV knowledge base [7].It can be seen that the theoretical values are close enough to the experimental ones.
The results show that even simple central two-body potentials provide good agreement with experimental data and thus may be used for calculations of the ground states.Moreover, we would like to emphasize the fact that good agreement with experimental data was obtained using the same unified set of potentials with the same parameters for all the studied nuclei, which makes the approach universal and self-consistent.

Table 1 .
Comparison of theoretical and experimental binding energies for the ground states of the studied nuclei