Application of CUDA technology for calculation of ground states of few-body nuclei by Feynman's continual integrals method*

Abstract

The possibility of application of modern parallel computing solutions to speed up the calculations of ground states of few-body nuclei by Feynman's continual integrals method has been investigated. These calculations may sometimes require large computational time, particularly in the case of systems with many degrees of freedom. This paper presents the results of application of general-purpose computing on graphics processing units (GPGPU). The energy and the square modulus of the wave function of the ground states of several few-body nuclei have been calculated using NVIDIA CUDA technology. The results show that the use of GPGPU significantly increases the speed of calculations.

Keywords

• CUDA
• Feynman's continual integrals method
• few-body nuclei

1. Introduction

2. Theory

(1) * The work was supported by grant 15-07-07673-a of the Russian Foundation for Basic Research (RFBR). is a propagator − the probability amplitude for the particle of mass m to travel from the point q0 to the point q in time t . Here S[q(t)] and Hˆ are the action and the Hamiltonian of the system, respectively, Dq(t) is the integration measure [8, 9]. For the time-independent potential energy the transition to the imaginary (Euclidean) time t  i gives the propagator KE q, ;q0 ,0 with the Euclidean action  1 KE q, ;q0 ,0   DEq() exp  SE q()

    m  dq 2  SE q()   d     V (q) .

0  2  d   Integration over q with the periodic boundary condition q  q0 allows to find the energy E0 of the ground state in the limit    [10]    Hˆ     KE q, ; q,0 dq  Sp exp    exp  En    exp  E  g(E)dE ,     n   0     KE q, ; q, 0 dq  exp   E0 ,    ,   

  E  KE q, ; q,0   n (q) 2 exp   En    E (q) 2 exp    g(E)dE .

n   Econt   Here g(E) is the density of states with the continuous spectrum E  Econt . For the system with a discrete spectrum and finite motion of particles the square modulus of the wave function of the ground state may also be found in the limit    [10] together with the energy E0

KE q, ;q,0  0 (q) 2 exp  E0 ,    .

 

Outside of the classically allowed region the square modulus of the wave function 0 (q) 2 of the ground state with E  Econt may be significantly less than E (q) 2 for the states with the continuous spectrum E  Econt . The ground state term in the formula (6) will not dominate despite the much more rapid decrease of the exponential factors exp     exp   E0 , E  E0 . Therefore, in this case  E      the formula (7) is in general applicable only for the region not far beyond the classically allowed ground state region.

Such situation may occur in the description of bound states of few-particle systems (for example, two protons and a neutron) when the existence of bound states of some of them (e.g., proton plus neutron) is possible.

The contribution of states with the continuum spectrum may be eliminated by introducing infinitely high walls in the potential energy located about the range of the nuclear forces beyond the classically allowed region. Introduction of the boundary condition 0 (q)  0 at these walls will not have a significant effect on the energy E0 and 0 (q) 2 far away from the walls.

Feynman’s continual integral (2) may be represented as the limit of the multiple integral K q, ; q0 ,0  lim 

N N

 1 N  mqk  qk1 2  exp     k1  2

 V qk  C N dq1dq2  dqN 1,

(8) (2) (3) (4) (5) (6) (7) where  m 1 2 qk  q(k ), k  k, k  0, N , qN  q, C    2   .

Here  N 1 -fold integral corresponds to averaging over the “path” of the particle as a broken line in the plane q,  with the vertices qk , k , k  1, N 1. For the approximate calculation of the continual integral (8) the continuous axis  is replaced by the grid   k  ka, k  0, N with the step a and the Euclidean propagator of a free particle K E(0) q, ; q0 ,0 is separated  a N  KE q, ; q0 ,0  K E(0) q, ; q0 ,0 exp  V (qk )  k1  0,N

, The denoted by angle brackets averaging over  N 1 -dimensional vectors Q  q1, , qN1 with the distribution law W q0;q1, , qN1;qN   m N qk  qk1 2  W q0 ; q1, , qN 1; qN   C N exp     k1 2a  may be calculated using the Monte Carlo method [15]. The standard algorithm for simulation of the random vector consists in a sequential choice of the values of its components from the conditional distributions W1 q1  , W2 q2 | q1  , W3 q3 | q1, q2  , …, WN1 qN1 | q1, q2 , , qN2  . Here Wk qk | q1, q2 , , qk1  is the probability density for the values of the quantity qk given the values of quantities q1, q2 , , qk1 . In this case the quantity qk is normally distributed with the mean value Mq , variance Dk and standard deviation k 

Dk [11] where The next trajectory in the simulation is calculated by the formula

Mqk  C1qk1  C2qN ,

Dk  C2 a k  C2 a m , m1 2 , C2   N  k 11 , C1   N  k C2 .

qk =Mqk  kk , k  1, N 1 , where k is a normally distributed random variable with zero mean and unity variance. Sample onedimensional random trajectories for low and large numbers of time steps are shown in Figs. 1a and 1b, respectively.

For large values of  random trajectories may reach the region where the probability density for the states with continuum spectrum is substantially larger than the probability density for the ground state, which may lead to a deviation from the asymptotic behavior (7) and the growth of the error. Therefore, the formula (7) is only applicable for the not very large values of  .

(9) ( 10 ) ( 11 ) ( 12 ) ( 13 ) ( 14 ) ( 15 ) ( 16 ) ( 17 )

For convenience of calculations in the scale of nuclear forces the expressions (6), ( 10 ) − ( 15 ) are represented using dimensionless variables q  q x0 , V  V (q) E0 , m  m m0 ,    t0 , a  a t0 , where x0  1 fm, E0  1 MeV, m0 is the nucleon mass, t0  m0 x02 1.57 1023 sec,

, Dk  k2 , qk =Mqk  kk , k  x0 C2 a m1 2 ,

 N  qk  qk 1 2  W  q0 ; q1, , qN 1; qN   C N exp   ,

 k 1 2a  ln KE q, ; q, 0 

ln 0 (q) 2  E0,    . probability density 0  r1,, rn 

Formulas (2)−( 16 ) are naturally generalized to a larger number of degrees of freedom and few particles including identical ones. The nucleon identity requires symmetrization of trajectories [9], therefore the configurations symmetric with respect to the positions of two neutrons and/or protons will be considered below. The nuclei 3H, 3He and 4He contain only two identical fermions (protons and/or neutrons with opposite spins) and the calculation of their ground states may be carried out without taking into account the Pauli principle by selecting pairs of identical fermions.

It should be noted that the calculation of multiple integrals required to find the multidimensional 2

by Feynman’s continual integrals method continues to be a challenging task. However, the analysis of the properties of 0  r1,, rn  2 allows to choose analytical approximations of 0 r1,, rn  and the application of the formula (7) in a single point in the multidimensional space allows to find the approximate value of the energy of the ground state.

To reduce the multiplicity of integrals in the formula ( 10 ) the calculation should be performed in the center of mass system using the Jacobi coordinates [4, 9].

For a system of two particles (2H nucleus) where r1 and r2 are the radius vectors of a proton and a neutron, respectively.

For a system of three particles, two of which are identical (2 neutrons or 2 protons in 3H and 3He nuclei, respectively) ( 18 ) ( 19 ) (20) (21) (22)

1 R  r2  r1, r  r3  r  r2  .

2 1 In the case of 3H nucleus r3 is the radius vector of a proton, r1 and r2 are the radius vectors of neutrons. In the case of 3He nucleus r3 is the radius vector of a neutron, r1 and r2 are the radius vectors of protons.

For a system of four particles consisting of two pairs of identical particles (2 protons and 2 neutrons in 4He nucleus) 1 1

R1  r2  r1, R2  r4  r3, r  2 r3  r4   2 r1  r2  , where r1 and r2 are the radius vectors of protons, r3 and r4 are the radius vectors of neutrons.

In the calculation of the propagator K q, ; q0 ,0 for the nuclei 2H, 3H, 3He, 4He neutron-proton Vn p (r) , neutron-neutron Vnn (r) and proton-proton Vp p (r) two-body strong interaction potentials have been used. The dependence of the nucleon-nucleon interaction with a hard core on the distance r was approximated by a combination of Gaussian type exponentials similar to the M3Y potential [16, 17] 3 Vnn (r)  Vp p (r)  uk exp  r 2 bk2  ,

k1 Vn p (r)  Vnn (r) . (23) (24) (25) (26) The values of the parameters u1  500 MeV, u2  102 MeV, u3  2 MeV, b1  0.59 fm, b2  1.40 fm, b3  2.94 fm and   1.2 MeV were determined from the condition of the absence of bound states of two identical nucleons as well as the approximate equality of the energy E0 found from (5) and (7) to the experimental values of the binding energies for the nuclei 2H, 3H, 3He, 4He (e.g., [18]). The plots of the potentials (25) and (26) are shown in Fig. 2.

3. Implementation

For numerical calculations the Monte Carlo method was used. The algorithm was developed and implemented in C++ programming language using NVIDIA CUDA technology.

The principal scheme of the calculation of the ground state energy for the one-dimensional case is shown in Fig. 3. The calculation of the propagator ( 18 ) is performed using L sequential launches of the kernel. Each kernel launch simulates n random trajectories in the space evolving from the Euclidean time   0 to  j  jt0 / a , where j  1, L (see Fig. 1). All trajectories start at the same point q0 in the space and in the moment  j return back to the same point q0 according to the probability distribution described above. The choice of the initial point q0 is arbitrary. In the case of the multidimensional space q0 must be replaced with the set of coordinates in the multidimensional space. All threads in a given kernel launch finish at approximately the same time, which makes the scheme quite effective in spite of the possible delays associated with the kernel launch overhead. Besides, the typical number of kernel launches L required for the calculation of the ground state energy usually does not exceed 100.

Starting from the certain time lin  Llint0 / a , where 1  Llin  L , the obtained values of the logarithm of the propagator b01 ln KE (21) tend to lie on the straight line, the slope of which gives the value of the ground state energy. The time lin is then used in the calculation of the square modulus of the wave function.

The principal scheme of the calculation of the square modulus of the wave function for the onedimensional case is shown in Fig. 4. Similarly, the calculation is performed using M sequential launches of the kernel. Each kernel launch simulates n random trajectories in the space from the Euclidean time   0 to the time lin determined in the calculation of the ground state energy. All trajectories start at the same point qi in the space and in the moment lin return back to the same point qi according to the probability distribution described above. Here i  1, M , where M is the total number of points in the space in which the square modulus of the wave function must be calculated. In the case of the multidimensional space qi must be replaced with the set of coordinates in the multidimensional space. One of the benefits of the approach is that the calculation may be easily resumed at a later time. For example, initially the square modulus of the wave function may be calculated with a large space step to obtain the general features of the probability distribution, and later new intermediate points are calculated and combined with those calculated previously. This may be very useful because the calculation of the square modulus of the wave function is generally much more time-consuming since it requires calculation in many points in the multidimensional space.

An important feature of the algorithm allowing to effectively use graphic processors is low consumption of memory during the calculation because it is not necessary to prepare a grid of values and store it in the memory.

To obtain normally distributed random numbers the cuRAND random number generator was used. According to the recommendations of the cuRAND developers each experiment was assigned a unique seed. Within the experiment, each thread of computation was assigned a unique sequence number. All threads between kernel launches were given the same seed, and the sequence numbers were assigned in a monotonically increasing way.

4. Results and discussion

Calculations were performed on the NVIDIA Tesla K40s accelerator installed within the heterogeneous cluster [19] of the Laboratory of Information Technologies, Joint Institute for Nuclear Research, Dubna. The code was compiled with CUDA version 7.5 for architecture version 3.5. Calculations were performed with single precision. The Euclidean time step a  0.01 was used. Additionally, NVIDIA GeForce 9800 GT accelerator was used for debugging and testing purposes.

The dependence of logarithm of the propagator b01 ln KE on the Euclidean time  is shown in Fig. 5 for nuclei 2H (a), 3H (b), 3He (c) and 4He (d). Different symbols correspond to different statistics n: empty circles (105), filled circles (106, 5·106, 107).

The behavior of the curves may be easily explained if we note that in all these cases only the energy of the ground state is negative and therefore only the first term in (4) increases with the increase of  , whereas the energies of the excited states are positive and hence the other terms in (4) decrease with the increase of  .

The results of linear fitting of the straight parts of the curves are shown in Fig. 5eh. According to the formula (21) the angular coefficient of the linear regression equals the binding energy. The obtained theoretical binding energies are listed in Tab. 1 together with the experimental values taken from [18]. It is clear that the theoretical values are close enough to the experimental ones, though obtaining good agreement was not the goal. As can be seen from Fig. 2, the difference between neutronneutron Vnn (r) and proton-proton Vp p (r) potentials is very small. Nevertheless, the difference between the calculated binding energies of 3H and 3He is observed in agreement with the experimental values.

The comparison of the square modulus of the wave function for 2H calculated on GPU using NVIDIA CUDA technology within Feynman’s continual integrals method and the square modulus of the wave function calculated on CPU within the shell model is shown in Fig. 6a. The same potentials (25), (26) were used. Good agreement between the curves confirms that the code based on Feynman’s continual integrals method using CUDA technology provides correct results.

It should be mentioned that the wave function cannot be measured directly, though the charge radii and charge distributions obtained from experiments may provide some information on its behavior. To compare the results of calculations with the experimental charge radii and charge distributions the wave function must be integrated.

The probability density distribution for the three-body configurations of 3He 0  R; r  2 (p + p + n) with   0 , 45 , 90 is shown in logarithmic scale in Fig. 7a,b,c, respectively, together with the potential energy surface (linear scale, lines). The vectors in Jacobi coordinates are shown in Fig. 7d.

The theoretical charge distribution for 3He obtained by integration of the wave function is compared with experimental data taken from [18] in Fig 6b. As can be seen, the agreement is very good. The obtained theoretical charge radius Rc2h 1/2  1.94 fm is also very close to the experimental value 1.9664  0.0023 fm. 2 The probability density distribution 0  R1; r ; R2   0  R1x ,0,0;0,0, rz ;0, R2 y  R1x ,0 2 for the symmetric tetrahedral configuration of four nucleons in the nucleus 4He

R1  r  R , R1  R , R1   R1x , 0, 0, r  0, 0, rz , R2  0, R2 y  R1x , 0 2 2 (27) is shown in logarithmic scale in Fig. 7e together with the potential energy surface (linear scale, lines). The vectors in Jacobi coordinates are shown in Fig. 7f.

Note also that the presence of the repulsive core in the nucleon-nucleon interaction reduces the probability of finding nucleons in the center of mass of the system for the considered symmetric configurations. This should lead to a smoother increase in the concentration of nucleons and the density of electric charge when approaching the center of the nucleus.

The analysis of the properties of 0  r1,, rn  2 allows to choose analytical approximations for it, e.g., as the product of the Gaussian type exponentials. The obtained approximations may be used in dynamic calculations.

The code implementing Feynman’s continual integrals method was initially written for CPU. The comparison of the calculation time of the ground state energy for 3He using Intel Core i5 3470 and NVIDIA Tesla K40s with different statistics is shown in Tab. 2. Even taking into account that the code for CPU used only 1 thread and a different random number generator, the time difference is impressive. This fact allows to increase the statistics and the accuracy of calculations in the case of using CUDA technology.

The comparison of the calculation time of the square modulus of the wave function for the ground state of 3He using Intel Core i5 3470 and NVIDIA Tesla K40s with the statistics 106 and the number of points in the space 60·60·12 is shown in Tab. 3. The value ~ 177 days for CPU is an estimation based on the performance gain in the calculation of the ground state energy. It is evident that beside the performance gain the use of CUDA technology may allow to reduce the space step in the calculation of the wave functions, as well as greatly simplify the process of debugging and testing, and in certain cases it may even enable calculations impossible before.

5. Conclusion

In this work an attempt is made to use modern parallel computing solutions to speed up the calculations of ground states of few-body nuclei by Feynman’s continual integrals method. The algorithm allowing to perform calculations directly on GPU was developed and implemented in C++ programming language. The method was applied to the nuclei consisting of nucleons, but it may also be applied to the calculation of cluster nuclei. The energy and the square modulus of the wave function of the ground states of several few-body nuclei have been calculated by Feynman’s continual integrals method using NVIDIA CUDA technology. The comparison with the square modulus of the wave function for 2H calculated on CPU within the shell model was performed to confirm the correctness of the calculations. The obtained values of the theoretical binding energies are close enough to the experimental values. The theoretical charge radius and charge distribution for 3He nucleus are also in good agreement with the experimental data. The results show that the use of GPGPU significantly increases the speed of calculations. This allows to increase the statistics and the accuracy of calculations as well as reduce the space step in calculations of wave functions. It also greatly simplifies the process of debugging and testing. In certain cases the use of CUDA enables calculations impossible before.

References

Penionzhkevich Yu.E. Reactions Involving Loosely Bound Cluster Nuclei: Heavy Ions and New Technologies // Phys. Atom. Nucl. 2011 . Vol. 74 . P. 1615 − 1622 .

Skobelev N.K. , Penionzhkevich Yu .E., Voskoboinik E.I. et al. Fusion and Transfer Cross Sections of 3He Induced Reaction on Pt and Au in Energy Range 10−24 .5 MeV // Phys. Part. Nucl.

Lett . 2014 . Vol. 11 . P. 208 − 215 .

Wu Y. , Ishikawa S. , Sasakawa T. Three-Nucleon Bound States: Detailed Calculations of 3H and 3He // Few-Body Systems . 1993 . Vol. 15 . P. 145 − 188 .

Dzhibuti R.I. , Shitikova K.V. Metod gipersfericheskikh funktsiy v atomnoy i yadernoy fizike [Method of Hyperspherical Functions in Atomic and Nuclear Physics]. Moscow, Energoatomizdat, 1993 . 269 P.

Kievsky A. , Marcucci L.E. , Rosati S . et al. High-Precision Calculation of the Triton Ground State Within the Hyperspherical- Harmonics Method // Few-Body Systems . 1997 . Vol. 22 . P. 1- 10 .

Viviani M. , Kievsky A. , Rosati S. Calculation of the ɑ -Particle Ground State // Few-Body Systems . 1995 . Vol. 18 . P. 25 − 39 .

Voronchev V.T. , Krasnopolsky V. М ., Kukulin V.I. A Variational Study of the Ground and Excited States of Light Nuclei in a Three-body Model on the Complete Basis . I. General Formalism // J. Phys. G. 1982 . Vol. 8 . P. 649 − 666 .

Feynman R.P. , Hibbs A.R. Quantum Mechanics and Path Integrals . New York, McGraw-Hill , 1965 . 382 P.

Blokhintsev D.I. Osnovy kvantovoy mekhaniki [Principles of Quantum Mechanics] . Moscow, Nauka, 1976 . 608 P.

10. Shuryak E.V. Stochastic Trajectory Generation by Computer // Sov. Phys. Usp. 1984 . Vol. 27 . P. 448 − 453 .

11. Samarin V.V. , Naumenko M.A. Izuchenie osnovnykh sostoyaniy nuklidov 3,4,6He metodom feynmanovskikh kontinual'nykh integralov [Study of Ground States of 3,4,6He Nuclides by Feynman's Continual Integrals Method] . Izv. RAN, Ser. Fiz. [Bull. Russ. Acad. Sci. Phys.] . 2016 . Vol. 80 , No. 3. P. 315 − 322 .

12. NVIDIA CUDA. URL: http://developer.nvidia.com/cuda-zone/ (accessed: 09 . 02 . 2016 ).

13. Sanders J. , Kandrot E. CUDA by Example: An Introduction to General-Purpose GPU Programming . New York, Addison-Wesley, 2011 . 290 P.

14. Perepyelkin E.E. , Sadovnikov B.I. , Inozemtseva N.G. Vychisleniya na graficheskikh protsessorakh (GPU) v zadachakh matematicheskoy i teoreticheskoy fiziki [Computing on Graphics Processors (GPU) in Mathematical and Theoretical Physics]. Moscow, LENAND, 2014 . 176 P.

15. Ermakov S.M. Metod Monte-Karlo v vychislitel'noy matematike: vvodnyy kurs [Monte Carlo Method in Computational Mathematics . Introductory Course]. St. Petersburg, Nevskiy Dialekt , 2009 . 192 P.

16. Satcher G.R. , Love W.G. Folding Model Potentials from Realistic Interaction for Heavy-Ion Scattering // Phys. Rep . 1979 . Vol. 55 , No. 3. P. 185 − 254 .

17. Alvarez M.A.G. , Chamon L.C. , Pereira D . et al. Experimental Determination of the Ion-Ion Potential in the N=50 Target Region: A Tool to Probe Ground-State Nuclear Densities // Nucl . Phys. A. 1999 . Vol. 656 , No. 2. P. 187 − 208 .

18. NRV Web Knowledge Base on Low-Energy Nuclear Physics . URL: http://nrv.jinr.ru/ (accessed: 09 . 02 . 2016 ).

19. Heterogeneous Cluster of LIT, JINR . URL: http://hybrilit.jinr.ru/ (accessed: 09 . 02 . 2016 ).