An efficient solution for Dirac equation in 3D lattice space with the conjugate gradient method
aa r X i v : . [ nu c l - t h ] J u l An efficient solution for Dirac equation in 3D lattice space withthe conjugate gradient method
B. Li ( 李 博 ), Z. X. Ren ( 任 政 学 ), and P. W. Zhao ( 赵 鹏 巍 ) ∗ State Key Laboratory of Nuclear Physics and Technology,School of Physics, Peking University, Beijing 100871, China
Abstract
An efficient method, preconditioned conjugate gradient method with a filtering function (PCG-F), is proposed for solving iteratively the Dirac equation in 3D lattice space for nuclear systems.The filtering function is adopted to avoid the variational collapsed problem and a momentum-dependent preconditioner is introduced to promote the efficiency of the iteration. The PCG-Fmethod is demonstrated in solving the Dirac equation with given spherical and deformed Woods-Saxon potentials. The solutions given by the inverse Hamiltonian method in 3D lattice space andthe shooting method in radial coordinate space are reproduced with a high accuracy. In comparisonwith the existing inverse Hamiltonian method, the present PCG-F method is much faster in theconvergence of the iteration, in particular for deformed potentials. It may also provide a promisingway to solve the relativistic Hartree-Bogoliubov equation iteratively in the future. ∗ [email protected] . INTRODUCTION During the past decades, new experimental facilities with radioactive beams have ex-tended our knowledge of nuclear physics from the stable to the unstable nuclei far from thestability line. The density functional theory (DFT) has been proved to be an importantmicroscopic approach for self-consistent description of nuclei [1]. Starting from a universalenergy density functional, DFT can provide a satisfactory description for nuclei all over thenuclide chart. The covariant density functional theory (CDFT), which includes the Lorentzsymmetry, has attracted a lot of attention in nuclear physics [2–6]. Its starting point isa standard effective Lagrangian density, where nucleons can be coupled with either finite-range meson fields [7, 8] or zero-range point-coupling interactions [9, 10]. The CDFT bringsmany advantages to describe the nuclear systems, such as the natural inclusion of spin-orbitinteractions [11] and the self-consistent treatment of the time-odd fields [12, 13].An essential ingredient of DFT is to solve the so-called Kohn-Sham equation for nucleons.In many cases, it has been solved by an expansion of a finite set of basis functions, suchas, the eigenfunctions of a harmonic oscillator. This method has been used successfullyfor many investigations in the literature, but it has limitations: (a) The convergence ofthe number of basis functions depends on the parameters of the basis, and this requires acareful optimization of these parameters, for example, in describing nuclear states with largedeformations; (b) For heavy nuclei, the required number of basis functions becomes large,and the construction and diagonalization of the single-particle Hamiltonian matrix witha large dimension in each step of the iteration leads to a surge of computational costs; (c)There are specific difficulties in describing nuclei with a large space distribution, for instance,in the case of halo nuclei [4, 14, 15]. Therefore, the methods developed in coordinate spaceare preferred to avoid these limitations.The CDFT has been successfully developed in spherical coordinate space by solving thespherical relativistic Hartree-Bogoliubov (RHB) equation, where the conventional shootingmethod works quite well [16]. It has been applied to investigate the halo and giant halophenomena in nuclei [17, 18]. Recently, this framework has also been used to explore thelimits of the nuclear landscape [19]. For deformed nuclei, however, the shooting methodbecomes rather complicated due to the difficulty in solving the coupled channel differentialequations [20]. Therefore, the Dirac Woods-Saxon (DWS) basis was proposed [21] and the2orresponding DWS basis expansion method has been used in the CDFT for deformed haloin nuclei [22–25]. The same basis has also been employed in the development of sphericaland axially deformed relativistic Hartree-Fock theories [26–28], where nonlocal potentialsare involved in. Despite these achievements, the CDFT calculations with further inclusionsof triaxial and/or octupole deformations become very sophisticated in the DWS expansionmethod. Therefore, to develop the CDFT in three-dimensional (3D) lattice space withoutany symmetry limitation is highly desired.The nonrelativistic DFT in 3D lattice space has been realized for a long time with theiterative methods including the imaginary time method [29, 30] and the damped-gradientmethod [31–33]. The basic idea of these iterative methods is searching for the descendingdirection of energy and following it iteratively until a local energy minimum is reached. Forcovariant DFT, however, due to the existence of Dirac sea in Dirac equation, the relativisticground state within the Fermi sea is a saddle point rather than a minimum. A directapplication of the iterative methods for Dirac equation usually encounters the so-calledvariational collapse problem [34, 35]. Therefore, the development of the CDFT on a 3Dlattice becomes a longstanding challenge. Recently, the inverse Hamiltonian method (IHM)has been used to avoid the variational collapse problem [36], and to solve the CDFT in 3Dlattice space [37]. Later on, the Fourier spectral method has been used to solve the fermiondoubling problem [38], another challenge in the numerical implementation of the IHM inlattice space. This new framework is then successfully applied to study the nuclear linear-chain [39] and toroidal structures [40]. Very Recently, the time-dependent CDFT has alsobeen developed on a 3D lattice and applied to investigate the microscopic dynamics of thelinear-chain cluster states [41].Although the IHM has been successful to solve the Dirac equations in 3D lattice space, itis still very time-consuming for heavy nuclei, due to the numerical complexity for calculatingthe inverse of the Dirac Hamiltonian. Moreover, it is practically difficult to apply the IHM toHartree-Fock-Bogoliubov (HFB) calculations in 3D lattice space due to the slow convergencein calculating the inverse of a HFB Hamiltonian [42], although it is feasible in the sphericalcase [43]. Therefore, it is desirable to develop an efficient method to solve the Dirac equationin 3D lattice space.In this paper, inspired by the successful application of the conjugate gradient methodwith a filtering step to solve the Dirac equation for electron systems [44], a preconditioned3onjugate gradient method with a filtering function (PCG-F) is developed to solve the Diracequation for nuclear systems on a 3D lattice. This new method avoids the inverse of Hamil-tonian and, thus, provides an efficient way to solve the nuclear Dirac equation in 3D latticespace. Moreover, it also paves a new way to solve the RHB equation with the powerfulgradient method [45]. The efficiency and accuracy of the newly proposed PCG-F methodare demonstrated in comparison with the inverse Hamiltonian and shooting methods.
II. THEORETICAL FRAMEWORKA. Conjugate gradient method for eigenstate problems
The conjugate gradient method was proposed to solve the system of linear equationsiteratively [46]. Later on, it was extended to solve the eigenstate problem [47], Aφ = λφ, (1)where A is a real n × n symmetric matrix and φ is the eigenstate with the eigenvalue λ . Inthe conjugate gradient method, the eigenstate with the smallest eigenvalue is obtained byminimizing the Rayleigh quotient,min X ∈ R n λ ( X ) = min X ∈ R n ( X, AX )( X, X ) . (2)Here, the trial solution X is updated iteratively starting from an normalized initial guess X (0) . In the i -th iteration, the search direction P ( i ) for updating X ( i ) is determined by P ( i ) = R ( i ) + β ( i − P ( i − , β ( i − = − ( R ( i ) , AP ( i − )( P ( i − , AP ( i − ) (3)with the residual R ( i ) = AX ( i ) − ( X ( i ) , AX ( i ) ) X ( i ) and P (0) = R (0) . As a result, the updated X ( i +1) is provided by X ( i +1) = G a X ( i ) + G b P ( i ) , (4)where the coefficients G a and G b are chosen to minimize λ ( i +1) = ( X ( i +1) , AX ( i +1) ) underthe normalization condition ( X ( i +1) , X ( i +1) ) = 1. B. Locally optimal block preconditioned conjugate gradient method
For a long time, the conjugate gradient method suffers from a poor convergence in theiteration process of finding the eigenstate. Therefore, the preconditioning technique has4een introduced, and this provides the preconditioned conjugate gradient (PCG) meth-ods [48]. Compared to other types of PCG methods, the so-called locally optimal blockPCG method [49], where the local optimization of a three-term recurrence is adopted, hasbeen shown to be effective for evaluating a relatively large number of eigenvalues and eigen-states.In the locally optimal block PCG method, the n lowest eigenstates of a Hamiltonian ˆ h are solved iteratively starting from a sets of normalized guess solutions X (0) k ( k = 1 , , ..., n ).The trial wavefunction X k is updated iteratively with X ( i +1) k = n X l =1 (cid:16) G akl X ( i ) l + G bkl R ( i ) l + G ckl P ( i ) l (cid:17) . (5)Here, R ( i ) l is the residual R ( i ) l = ˆ hX ( i ) l − h X ( i ) l | ˆ h | X ( i ) l i X ( i ) l , (6)and P ( i ) l is the previous search direction P ( i ) l = X ( i ) l − n X l ′ =1 h X ( i − l ′ | X ( i ) l i X ( i − l ′ , P (0) = 0 . (7)To accelerate the convergence of the evolution, the preconditioning technique is usuallyadopted for R ( i ) W ( i ) l = ˆ T − l R ( i ) l , (8)where ˆ T l is the preconditioner. As a result, the wavefunction X k is updated with X ( i +1) k = n X l =1 (cid:16) G akl X ( i ) l + G bkl W ( i ) l + G ckl P ( i ) l (cid:17) , (9)where the coefficient matrices G a , G b and G c are chosen to minimize n P k =1 h X ( i +1) k | ˆ h | X ( i +1) k i under the orthonormalization condition h X ( i +1) k | X ( i +1) l i = δ kl . C. Filtering and preconditioning operators for Dirac equation
The main task for CDFT is to solve the Dirac equation with the Hamiltonian,ˆ h = α · ˆ p + β ( m N + S ) + V − m N , (10)where α and β are the Dirac matrices, m N is the mass of nucleon, and S and V are thescalar and vector potentials, respectively. For the sake of convenience, here the Hamiltonian5s shifted down by a nucleon mass m N . Since the spectrum of the Dirac Hamiltonian ˆ h contains negative- and positive-energy states, a direct application of the PCG method wouldsuffer from the variational collapse problem.To avoid the variational collapse problem, a filtering operator can be used to suppressthe components of negative-energy states in the wavefunctions during the iteration [44]. Inthe present work, the filtering operator is taken as, F (ˆ h ) = 1 D (ˆ h − C ) , (11)where C and D are two parameters to be optimized in the practical calculations. Thefiltering operator is implemented in the PCG method by replacing X ( i ) l and W ( i ) l in Eq. (9)with X ( i ) l → F (ˆ h ) X ( i ) l , W ( i ) l → [ F (ˆ h )] N F W ( i ) l . (12)Note that here the filtering operation on W ( i ) l is carried out by N F times. This is differentfrom the solution of Dirac equation for electron systems, where the states in the Fermi- andDirac sea are well separated due to the negligible spin-orbit splittings. For nuclear systems,however, due to the large spin-orbit interactions, the energy gap between the positive- andnegative-energy states is only two or three times of the potentials S and V . The optimizedvalue of N F will be discussed below in Sec. III.The multiple filtering operations on W ( i ) l could lead to a poor convergence behaviorbecause the components of high positive-energy states can be substantially enlarged. Moti-vated by the fact that the high-energy states are usually dominated by the kinetic energy,here the preconditioner ˆ T in Eq. (8) is chosen as the following momentum-dependent form,ˆ T l = [ ˆ p + g l m N ] , (13)where g l is a optimized factor and will be discussed below in Sec. III. This preconditioneroperation could effectively damp the components of high-energy states, and provides anefficient convergence. Note that for electron systems, such a preconditioner operation is notmandatory [44]. 6 II. NUMERICAL DETAILS
In the present work, the Dirac equation for nucleons is solved in 3D lattice space by thePCG-F method. The large scalar and vector potentials are taken as a Woods-Saxon form, V ( r ) + S ( r ) = V r − R D ( θ, ϕ )) /a ] ,V ( r ) − S ( r ) = − λV r − R ls D ( θ, ϕ )) /a ls ] , (14)where D ( θ, ϕ ) brings in the quadrupole ( β , γ ) and octupole deformations β , D ( θ, ϕ ) =1 + β cos γY ( θ, ϕ ) + 1 √ β sin γ [ Y ( θ, ϕ )+ Y − ( θ, ϕ )] + β Y ( θ, ϕ ) . (15)The parameters for the Woods-Saxon potentials are listed in Table I, which correspond tothe neutron potentials of Ca [50].
TABLE I. The parameters in the Woods-Saxon potentials [see Eq. (14)]. V [MeV] R [fm] a [fm] λ R ls [fm] a ls [fm]-65.796 4.482 0.615 11.118 4.159 0.648 In the present calculations, the coordinate space along the x , y , and z axes is respectivelydiscretized by 28 grids with the mesh size d = 1 fm. The initial guess of the single-particlewavefunctions are taken as the spherical harmonic oscillator wavefunctions for both upperand lower components. To avoid the fermion doubling problem, the spatial derivatives areperformed in the momentum space with the help of the fast Fourier transformation [38].For the filtering operation, we define C = − m N in F (ˆ h ), and this is in analogy toRef. [44], where C = − m e is used for electron systems. It should be noted that the valueof D is in principle irrelative, because it is anyhow absorbed in the coefficient matrices G a , G b , G c in Eq. (9) via the energy minimization and orthonormalization condition. In thepresent work, we define D = ( V + S ) min + 2 m N with ( V + S ) min being the minimum of thepotential V + S in the coordinate space. By this definition, we have F (ˆ h ) ≈ F (ˆ h ) ≈ . F ( E )is shown as a function of energy E . It can be seen that the F ( E ) values are very small7around 10 − ) in the region of negative-energy spectrum, in comparison with those in thepositive-energy spectrum. This could suppress the negative-energy components during theiteration, while the suppression is found to be not sufficient to avoid the variational collapsein the practical calculations. Therefore, as in Eq. (12), the filtering operation on W ( i ) l isperformed by N F times.In Figs. 1(b) and 1(c), the square and fourth power of the filtering function F ( E ) areshown respectively. The suppression of the filtering function on negative-energy states ispromoted obviously with the increasing power. In particular for F ( E ), the suppression onthe negative-energy states reaches to an oder of 10 − . It is found that such a suppressionworks quite well to overcome the variational collapse problem and, thus, we adopt N F = 4in the present work.In Fig. 1, one can also see that the filtering function at high energies becomes larger withthe increasing power. This could lead to a slow convergence of the iteration, because thecomponents of bound states are relatively reduced by the filtering function. Therefore, asin Eq. (13), the preconditioner ˆ T l is introduced, and the factor g l is taken as g l = 0 . λ l − ( V + S ) min ( V + S ) min + 0 . , (16)with the single-particle energy λ l = h X l | ˆ h | X l i .To demonstrate the effects of the preconditioner, the PCG-F method is applied for a Diracequation with a spherical Woods-Saxon potential in 3D lattice space by either considering thepreconditioner ˆ T l as in Eq. (13) or setting ˆ T l as a unit operator, i.e., without preconditoner.The evolution of the maximum energy dispersion h ˆ h i − h ˆ h i for the bound single-particlestates is shown in Fig. 2. For the PCG-F method, it takes only 15 iterations to reducethe maximum energy dispersion to 10 − MeV , while without the preconditioner, it takesmore than 1200 iterations to reach the same level. As a comparison, the IHM calculationsrequire more than 30 iterations to reach this accuracy. Moreover, one can see that theenergy dispersions can drop to around 10 − MeV after 18 iterations for the PCG-F method(after 37 iterations for the IHM), but they are finally fluctuated around 10 − MeV for thecalculations without the preconditioner. Therefore, one can conclude that the preconditionergreatly improves the convergent accuracy and the speed of the iteration.In the previous work with the IHM [38], the convergence of a wavefunction is regardedto be reached if the corresponding energy dispersion is smaller than 10 − MeV . The same8 .00.51.01.50.00.51.01.5 -2500 -2000 -1500 -10000.000.020.04 -2500 -2000 -1500 -1000 -500 0 5000.00.51.01.5 -2500 -2000 -1500 -10000.04.0x10 -4 -4 (a) F(E)
Positive-energy states Negative-energy states F n ( E ) (b) F (E) E [MeV] (c) F (E) FIG. 1. (Color online) The filtering function F n ( E ) with n = 1 (top), n = 2 (middle), and n = 4(bottom) as a function of energy E . The solid and dashed lines respectively correspond to thepositive- and negative-energy states, which are separated by the filtering function values (dottedlines). The insets present a partial enlargement for the negative-energy region. criterion is adopted in the present work, while W ( i ) l or P ( i ) l are removed from the summationin Eq. (9) if the corresponding energy dispersion h X l | ˆ h | X l i − h X l | ˆ h | X l i is smaller than10 − MeV or 10 − MeV , respectively. IV. RESULTS AND DISCUSSIONA. Spherical potential
We first assume the potentials S and V in Eq. (14) are spherical, and solve the corre-sponding Dirac equation with the PCG-F method. The results of other methods, includingthe IHM in 3D lattice space and the shooting method in radial coordinate space, are usedfor comparison. The numerical details used in the IHM is the same as those in the PCG-9
10 20 30 40 50 6010 -11 -9 -7 -5 -3 -1 [ Æ h æ - Æ h æ ] m a x [ M e V ] Iteration number N iter.
PCG-F CG-F (N iter. /30) IHM
FIG. 2. (Color online) Evolution of the maximum energy dispersion for the bound single-particlestates in the spherical Woods-Saxon potential as a function of the iteration number. The solidand dashed lines respectively represent the results with and without the preconditioner, and theabscissa of the latter is scaled by a factor of 30. The results of the inverse Hamiltonian method(dotted line) are also shown for comparison.
F method. For the shooting method, the radial box size R = 20 fm and the mesh size dr = 0 .
01 fm are adopted, and the obtained results can be regarded as exact solutionsthanks to the high accuracy. S i ng l e - pa r t i c l e E ne r g y [ M e V ] Iteration numbers
Shooting IHM
FIG. 3. (Color online) Evolution of the single-particle energies in a spherical Woods-Saxon potentialobtained by the PCG-F method as a function of the iteration numbers. For comparison, the single-particle energies obtained by the inverse Hamiltonian and the shooting methods are shown in theright side together with the spherical quantum numbers.
In Fig. 3, the evolution of single-particle energies obtained by the PCG-F method is shownas a function of the iteration numbers. There are 40 bound states obtained in the PCG-F10ethod, and they are grouped in energy according to the degeneracy due to the sphericalsymmetry. One can see that the single-particle energies given by the PCG-F method arein very good agreement with those give by the IHM and the shooting method. Althoughthe evolution in Fig. 3 is shown up to 20 iterations, to check the numerical stability of thePCG-F methods, the calculation was carried out up to the 3000th iteration. It is foundthat the obtained single-particle energies of the bound states are quite stable after the 15thiteration, and the negative-energy states are eliminated perfectly.
FIG. 4. (Color online) Absolute differences between the single-particle energies in a sphericalWoods-Saxon potential obtained with the PCG-F method and those with other methods includingthe inverse Hamiltonian (a) and shooting (b) methods. The spherical quantum numbers are listedin panel (b). For the degenerate states, only the maximum and minimum deviations are shownwith the connecting bands.
For a more precise comparison, Fig. 4 shows the absolute differences between the single-particle energies obtained with the PCG-F method and those with the inverse Hamiltonianand shooting methods. In Fig. 4(a), one can see that the absolute energy deviations betweenthe PCG-F method and IHM are extremely small for all states, i.e., less than 10 − MeV.This demonstrates that the 3D lattice calculations with these two methods are in accuracyat the same level.In Fig. 4(b), the absolute energy deviations between the PCG-F and shooting methodsare found to be in the range of 10 − ∼ − MeV. This also shows the high accuracy of the11CG-F method, and it can be further improved by reducing the mesh size and/or enlargingthe 3D box size. Moreover, in contrast to the shooting solutions, the spherical degeneracyof the single-particle levels are slightly broken in the 3D lattice calculations. This is becausethe discretized 3D lattice space is not exactly spherical, while the spherical symmetry isexactly fulfilled for the shooting method on radial coordinates. -7 -6 -5 -4 -3 -2 -1 D en s i t y [f m - ] Shooting PCG-F CaNeutron (a) r [fm] (b)
FIG. 5. (Color online) The total density of the lowest 28 levels in the spherical Woods-Saxonpotential as a function of the radial coordinate r in normal (top) and logarithmic scales (bottom).The circles and lines represent the results of the PCG-F and shooting methods, respectively. Apart from the single-particle energies, it is necessary to examine the accuracy of thewavefunctions. This is shown in Fig. 5, where the total density of the lowest 28 levels, i.e.,the neutron density for Ca, are shown as a function of the radial coordinate r in normaland logarithmic scales, respectively. The density obtained with the PCG-F method agreeswith that with the shooting method very well, even at very large r values. In the region of r > B. Deformed potential
Deformation can be introduced to the potentials S and V through the parameters( β, γ, β ) in Eq. (15). We investigate three cases with the PCG-F method, i.e., ( β, γ, β ) =120 . , ◦ , . , ◦ , . , ◦ , . -11 -8 -5 -2 -11 -8 -5 -2 PCG-F (b, g, b )= (0.0, 0(cid:176), 0) (0.3, 0(cid:176), 0) (0.3, 30(cid:176), 0) (0.3, 30(cid:176), 0.1) (a) IHM [ Æ h æ - Æ h æ ] m a x [ M e V ] Iteration number N iter. (b)
FIG. 6. (Color online) The maximum energy dispersions of the bound single-particle states invarious Woods-Saxon potentials obtained by the PCG-F (top) and inverse Hamiltonian methods(bottom) as a function of the iteration number.
Firstly, the convergent behavior of the PCG-F and inverse Hamiltonian methods is exam-ined for spherical and deformed potentials in Fig. 6, where the maximum energy dispersionsof the bound single-particle states are shown as a function of the iteration number. Fordeformed potentials, the convergence for both the PCG-F and inverse Hamiltonian meth-ods becomes slower. This is because, on the one hand, there are more bound states indeformed potentials than in spherical ones, and, on the other hand, the initial guess of thewavefunctions is usually spherical for simplicity.Nevertheless, for all potentials, the iteration of the PCG-F method is more efficient thanthe inverse Hamiltonian method. In particular for deformed potentials, for instance, it takesless than 30 iterations for the PCG-F method to reduce the maximum energy dispersion toaround 10 − MeV , while for the inverse Hamiltonian method, one needs more than 100iterations. This feature should be helpful to save the computational time for the futurestudies on many phenomena with the 3D lattice CDFT, such as exotic deformations [13, 51,52], super-heavy nuclei [53–55], fission [56–58], fusion dynamics [41, 59, 60], etc.13 -14 -13 -12 -11 -14 -13 -12 -11 -50 -40 -30 -20 -10 010 -15 -14 -13 -12 -11 (b, g, b )=(0.3, 0(cid:176), 0) (a) (b, g, b )=(0.3, 30(cid:176), 0) (b) | E L O BP C G - F - E I H M | [ M e V ] Single-particle energy [MeV] (b, g, b )=(0.3, 30(cid:176), 0.1) (c) FIG. 7. (Color online) Absolute differences between the single-particle energies in deformed Woods-Saxon potentials obtained with the PCG-F method and those with the inverse Hamiltonian method.Panels (a), (b), and (c) present the results in axial, triaxial, and triaxially octupole potentials,respectively.
In Fig. 7, the absolute differences between the single-particle energies obtained with thePCG-F method and those with the inverse Hamiltonian method are depicted. Similar to thespherical case [see Fig. 4(a)], the results of the two methods agree with each other in a veryhigh accuracy, i.e., less than 10 − MeV. In particular, the magnitudes of the deviations arenot affected by the shape of the potential. This demonstrates that the 3D lattice calculationsrealized by the PCG-F and the inverse Hamiltonian methods are essentially equivalent.However, considering the high efficiency of the PCG-F method, to develop the CDFT in 3Dlattice space with this method would be very beneficial in the future.Moreover, it is worthwhile to mention the perspective of the present PCG-F method onsolving the RHB equation, where pairing correlations are taken into account with the Bogoli-ubov transformation. Instead of diagonalizing huge matrices in the basis expansion method,the nonrelativistic HFB equation can be solved by the powerful gradient method, which is14uite robust and easily deals with multiple constraints [45]. However, a direct application ofthe gradient method for the RHB equation is inhibited by the quasiparticle states in Diracsea. In this sense, the present PCG-F method, with an appropriate filtering function, seemsa promising way to solve the relativistic Hartree-Bogoliubov equation iteratively.
V. SUMMARY
In summary, an efficient method, PCG-F, has been proposed for solving nuclear Diracequation in 3D lattice space, where a filtering function is adopted to avoid the variationalcollapsed problem and a momentum-dependent preconditioner is introduced to promote theefficiency of the iteration. The method has been demonstrated in solving the Dirac equationwith spherical and deformed Woods-Saxon potentials. In the spherical case, the PCG-F method reproduces the single-particle energies and densities obtained by the shootingmethod in radial coordinate space with a high accuracy. In both spherical and deformedcases, the single-particle energies obtained with the PCG-F and inverse Hamiltonian methodsagree with each other very precisely, but the PCG-F method is much faster to achievethe convergence of the iteration, in particular for deformed potentials. Considering thehigh efficiency of the PCG-F method, to develop the CDFT in 3D lattice space with thismethod would be very beneficial in the future. Moreover, the present PCG-F method seemsa promising way to solve the relativistic Hartree-Bogoliubov equation iteratively. Worksfollowing these directions are in progress.
ACKNOWLEDGMENTS
This work was partly supported by the National Key R&D Program of China (ContractsNo. 2018YFA0404400 and 2017YFE0116700), the National Natural Science Foundationof China (Grants No. 11621131001, 11875075, 11935003, and 11975031), the State KeyLaboratory of Nuclear Physics and Technology, Peking University (No. NPT2020ZZ01),and the China Postdoctoral Science Foundation under Grant No. 2020M670013. [1] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. , 121 (2003).
2] J. Meng, ed., Relativistic Density Functional for Nuclear Structure, International Review ofNuclear Physics, Vol. 10 (World Scientific, Singapore, 2016).[3] P. Ring, Prog. Part. Nucl. Phys. , 193 (1996).[4] J. Meng, H. Toki, S. G. Zhou, S. Q. Zhang, W. H. Long, and L. S. Geng,Prog. Part. Nucl. Phys. , 470 (2006).[5] D. Vretenar, A. V. Afanasjev, G. A. Lalazissis, and P. Ring, Phys. Rep. , 101 (2005).[6] T. Nikˇsi´c, D. Vretenar, and P. Ring, Prog. Part. Nucl. Phys. , 519 (2011).[7] W. Long, J. Meng, N. Van Giai, and S.-G. Zhou, Phys. Rev. C , 034319 (2004).[8] G. A. Lalazissis, T. Nikˇsi´c, D. Vretenar, and P. Ring, Phys. Rev. C , 024312 (2005).[9] T. Nikˇsi´c, D. Vretenar, and P. Ring, Phys. Rev. C , 034318 (2008).[10] P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, Phys. Rev. C , 054319 (2010).[11] M. M. Sharma, G. A. Lalazissis, and P. Ring, Phys. Lett. B , 9 (1993).[12] J. Meng, J. Peng, S.-Q. Zhang, and P.-W. Zhao, Front. Phys. , 55 (2013).[13] P. Zhao and Z. Li, Int. J. Mod. Phys. E , 1830007 (2018).[14] I. Tanihata, H. Hamagaki, O. Hashimoto, Y. Shida, N. Yoshikawa, K. Sugimoto, O. Yamakawa,T. Kobayashi, and N. Takahashi, Phys. Rev. Lett. , 2676 (1985).[15] J. Meng and S.-G. Zhou, J. Phys. G , 093101 (2015).[16] J. Meng, Nucl. Phys. A , 3 (1998).[17] J. Meng and P. Ring, Phys. Rev. Lett. , 3963 (1996).[18] J. Meng and P. Ring, Phys. Rev. Lett. , 460 (1998).[19] X. W. Xia, Y. Lim, P. W. Zhao, H. Z. Liang, X. Y. Qu, Y. Chen, H. Liu, L. F. Zhang, S. Q.Zhang, Y. Kim, and J. Meng, At. Data Nucl. Data Tables , 1 (2018).[20] C. E. Price and G. E. Walker, Phys. Rev. C , 354 (1987).[21] S.-G. Zhou, J. Meng, and P. Ring, Phys. Rev. C , 034323 (2003).[22] S.-G. Zhou, J. Meng, P. Ring, and E.-G. Zhao, Phys. Rev. C , 011301 (2010).[23] L. Li, J. Meng, P. Ring, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C , 024312 (2012).[24] X.-X. Sun, J. Zhao, and S.-G. Zhou, Phys. Lett. B , 530 (2018).[25] K. Y. Zhang, D. Y. Wang, and S. Q. Zhang, Phys. Rev. C , 034312 (2019).[26] W. H. Long, P. Ring, N. V. Giai, and J. Meng, Phys. Rev. C , 024308 (2010).[27] W. H. Long, P. Ring, J. Meng, N. Van Giai, and C. A. Bertulani,Phys. Rev. C , 031302(R) (2010).
28] J. Geng, J. Xiang, B. Y. Sun, and W. H. Long, Phys. Rev. C , 064302 (2020).[29] K. T. R. Davies, H. Flocard, S. Krieger, and M. S. Weiss, Nucl. Phys. A , 111 (1980).[30] P. Bonche, H. Flocard, and P. Heenen, Comput. Phys. Commun. , 49 (2005).[31] P.-G. Reinhard and R. Cusson, Nucl. Phys. A , 418 (1982).[32] J. A. Maruhn, P. G. Reinhard, P. D. Stevenson, and A. S. Umar,Comput. Phys. Commun. , 2195 (2014).[33] Y. Y. Wang and Z. X. Ren, Sci. China-Phys. Mech. Astron. , 082012 (2018).[34] Y. Zhang, H. Z. Liang, and J. Meng, Chin. Phys. C , 113 (2009).[35] Y. Zhang, H. Liang, and J. Meng, Int. J. Mod. Phys. E , 55 (2010).[36] K. Hagino and Y. Tanimura, Phys. Rev. C , 057301 (2010).[37] Y. Tanimura, K. Hagino, and H. Z. Liang, Prog. Theor. Exp. Phys. , 073D01 (2015).[38] Z. X. Ren, S. Q. Zhang, and J. Meng, Phys. Rev. C , 024313 (2017).[39] Z. X. Ren, S. Q. Zhang, P. W. Zhao, N. Itagaki, J. A. Maruhn, and J. Meng,Sci. China-Phys. Mech. Astron. , 112062 (2019).[40] Z. X. Ren, P. W. Zhao, S. Q. Zhang, and J. Meng, Nucl. Phys. A , 121696 (2020).[41] Z. X. Ren, P. W. Zhao, and J. Meng, Phys. Lett. B , 135194 (2020).[42] Y. Tanimura, Covariant Density Functional Calculations for Atomic Nuclei in the 3-dimensional Coordinate Space,Ph.D. thesis, Department of Physics, Tohoku University (2014).[43] Y. Tanimura, K. Hagino, and P. Ring, Phys. Rev. C , 017301 (2013).[44] L. Lin, S. Shao, and W. E, J. Comput. Phys. , 205 (2013).[45] P. Ring and P. Schuck, The nuclear many-body problem (Springer Science & Business Media,2004).[46] M. R. Hestenes and E. Stiefel, J. Res. Natl. Bur. Stan. , 409 (1952).[47] W. W. Bradbury and R. Fletcher, Numer. Math. , 259 (1966).[48] A. V. Knyazev, Electron. Trans. Numer. Anal. , 104 (1998).[49] A. V. Knyazev, SIAM J. Sci. Comput. , 517 (2001).[50] W. Koepf and P. Ring, Z. Phys. A , 81 (1991).[51] P. W. Zhao, N. Itagaki, and J. Meng, Phys. Rev. Lett. , 022501 (2015).[52] J. Zhao, B.-N. Lu, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C , 014320 (2017).[53] S. E. Agbemava, A. V. Afanasjev, T. Nakatsukasa, and P. Ring,Phys. Rev. C , 054310 (2015).
54] Z. Shi, A. V. Afanasjev, Z. P. Li, and J. Meng, Phys. Rev. C , 064316 (2019).[55] X. Meng, B. Lu, and S. Zhou, Sci. China-Phys. Mech. Astron. , 212011 (2020).[56] B.-N. Lu, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C , 011301(R) (2012).[57] B.-N. Lu, J. Zhao, E.-G. Zhao, and S.-G. Zhou, Phys. Rev. C , 014323 (2014).[58] S.-G. Zhou, Phys. Scr. , 063008 (2016).[59] A. Umar and V. Oberacker, Nucl. Phys. A , 238 (2015).[60] L. Guo, C. Shen, C. Yu, and Z. Wu, Phys. Rev. C , 064609 (2018)., 064609 (2018).