Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis (III) hfbtho (v3.00): a new version of the program

TL;DR

The paper introduces version 3.00 of the HFBTHO code, enhancing nuclear structure calculations with new features like Gogny force implementation, fission fragment analysis, and MPI support for large-scale studies.

Contribution

This update integrates comprehensive nuclear interaction models and advanced computational features into the HFBTHO code for improved nuclear physics simulations.

Findings

Implementation of full Gogny force in HFB calculations

Calculation of nuclear fission fragment properties

Enhanced computational capabilities with MPI interface

Abstract

We describe the new version 3.00 of the code HFBTHO that solves the nuclear Hartree-Fock (HF) or Hartree-Fock-Bogolyubov (HFB) problem by using the cylindrical transformed deformed harmonic oscillator basis. In the new version, we have implemented the following features: (i) the full Gogny force in both particle-hole and particle-particle channels, (ii) the calculation of the nuclear collective inertia at the perturbative cranking approximation, (iii) the calculation of fission fragment charge, mass and deformations based on the determination of the neck (iv) the regularization of zero-range pairing forces (v) the calculation of localization functions (vi)MPI interface for large-scale mass table calculations.

Figures

Tables

Table 1. Table 1: Benchmark of the DFT solvers hfbtho and hfodd for a Hartree-Fock-Bogolyubov calculation in 120 Sn with the D1S Gogny functional. We choose on purpose a deformed basis with N max = 20 subscript 𝑁 max 20 N_{\rm max}=20 and β = 0.4 𝛽 0.4 \beta=0.4 . The spherical-equivalent of the oscillator length is b 0 = 2.1339661000 subscript 𝑏 0 2.1339661000 b_{0}=2.1339661000 fm ( b ⟂ = 2.0035166 subscript 𝑏 perpendicular-to 2.0035166 b_{\perp}=2.0035166 fm and b z = 2.4208988 subscript 𝑏 𝑧 2.4208988 b_{z}=2.4208988 fm). The deformation of the initial solution is set to the same value as the basis deformation, that is β 0 = 0.4 subscript 𝛽 0 0.4 \beta_{0}=0.4 .

	hfbtho	hfodd
$E_{\mathrm{tot}}$ (MeV)	-1385.166990	-1385.166981
$E_{\mathrm{kin}}^{\left(n\right)}$ (MeV)	1416.363321	1416.363301
$E_{\mathrm{kin}}^{\left(p\right)}$ (MeV)	888.527464	888.527452
$E_{\mathrm{Vol}}$ (MeV)	4004.862420	4004.862313
$E_{\mathrm{SO}}$ (MeV)	11-72.659299	11-72.659315
$E_{\mathrm{Gogny}}^{\left(\mathrm{dir}\right)}$ (MeV)	-7274.909836	-7274.909691
$E_{\mathrm{Gogny}}^{\left(\mathrm{exc}\right)}$ (MeV)	-326.636358	-326.636390
$E_{\mathrm{Gogny}}^{\left(\mathrm{pair}\right)}$ (MeV)	11-20.714701	11-20.714651
$r_{\mathrm{rms}}^{\left(n\right)}$ (fm)	-1114.626410	-1114.626410
$r_{\mathrm{rms}}^{\left(n\right)}$ (fm)	-1114.438686	-1114.438686

Table 2. Table 2: Benchmark of the DFT solvers hfbtho and hfodd for a constrained Hartree-Fock-Bogolyubov calculation in 152 Dy with the D1S Gogny functional. We choose a deformed basis with N max = 14 subscript 𝑁 max 14 N_{\rm max}=14 , β = 0.3 𝛽 0.3 \beta=0.3 . The spherical-equivalent of the oscillator length is b 0 = 2.12100948454 subscript 𝑏 0 2.12100948454 b_{0}=2.12100948454 fm ( b ⟂ = 2.0230038 subscript 𝑏 perpendicular-to 2.0230038 b_{\perp}=2.0230038 fm and b z = 2.3314948 subscript 𝑏 𝑧 2.3314948 b_{z}=2.3314948 fm). The deformation of the initial solution is set to the same value as the basis deformation, that is β 0 = 0.4 subscript 𝛽 0 0.4 \beta_{0}=0.4 . We impose a constraint of ⟨ Q ^ 20 ⟩ = 20 delimited-⟨⟩ subscript ^ 𝑄 20 20 \langle\hat{Q}_{20}\rangle=20 b.

	hfbtho	hfodd
$E_{\mathrm{tot}}$ (MeV)	-1816.287980	-1816.287978
$E_{\mathrm{kin}}^{\left(n\right)}$ (MeV)	1742.872954	1742.873034
$E_{\mathrm{kin}}^{\left(p\right)}$ (MeV)	1199.286139	1199.286169
$E_{\mathrm{Vol}}$ (MeV)	5224.217198	5224.217576
$E_{\mathrm{SO}}$ (MeV)	11-90.243411	11-90.243458
$E_{\mathrm{Gogny}}^{\left(\mathrm{dir}\right)}$ (MeV)	-9463.410240	-9463.410712
$E_{\mathrm{Gogny}}^{\left(\mathrm{exc}\right)}$ (MeV)	-411.011819	-411.011793
$E_{\mathrm{Gogny}}^{\left(\mathrm{pair}\right)}$ (MeV)	11-17.998801	11-17.998793
$r_{\mathrm{rms}}^{\left(n\right)}$ (fm)	-1115.103611	-1115.103611
$r_{\mathrm{rms}}^{\left(n\right)}$ (fm)	-1114.959125	-1114.959125

Table 3. Table 3: Benchmark of the collective inertia mass tensor and zero-point energy corrections at the perturbative cranking approximation; see text for details.

		hfbtho	hfodd
GCM	M22 (MeV-1b-2)	0.2016.10-2	0.2016.10-2
	M32 (MeV-1b-5/2)	-0.5787.10-3	-0.5787.10-3
	M33 (MeV-1b-3)	0.3878.10-2	0.3878.10-2
	${ϵ}_{\text{ZPE}}$ (MeV)	2.1785	2.1784
ATDHFB	M22 (MeV-1b-2)	0.2822.10-2	0.2826.10-2
	M32 (MeV-1b-5/2)	-0.8667.10-3	-0.8667.10-3
	M33 (MeV-1b-3)	0.4848.10-2	0.4849.10-2
	${ϵ}_{\text{ZPE}}$ (MeV)	1.6819	1.6819

Table 4. Table 4: Benchmark of the characteristics of the fission fragments at the point ⟨ q 20 ⟩ = 320 delimited-⟨⟩ subscript 𝑞 20 320 \langle q_{20}\rangle=320 b and ⟨ q 30 ⟩ = 25 delimited-⟨⟩ subscript 𝑞 30 25 \langle q_{30}\rangle=25 b 3/2 . We list the dimensionless coordinate of the neck ξ N = z N / b z subscript 𝜉 𝑁 subscript 𝑧 𝑁 subscript 𝑏 𝑧 \xi_{N}=z_{N}/b_{z} . z CM1 subscript 𝑧 CM1 z_{\mathrm{CM}1} (resp. z CM2 subscript 𝑧 CM2 z_{\mathrm{CM}2} ) give the position of the center of mass of fragment 1 (resp., fragment 2) and is given in fermis; see text for additional details.

	hfbtho	hfodd
${\xi }_N$	-110.6418	-110.6418
qN	-118.8442	-118.8442
zCM1 (fm)	11-6.3014	11-6.3013
zCM2 (fm)	-119.1059	-119.1059
Z1	56.0055	56.0053
Z2	37.9945	37.9947
A1	141.8508	141.8480
A2	98.1492	98.1397
q20,1 (b)	29.2889	29.2793
q20,2 (b)	15.2483	15.2234
q30,1 (b${}_{}{}^{3/2})$	-111.6485	-111.6521
q30,2 (b${}_{}{}^{3/2})$	11-0.3751	11-0.3832

Equations

$$\hat{V}({\mathbf{r}}_1,{\mathbf{r}}_2)$$

$$+t_0\left(1+x_0{\hat{P}}_{\sigma }\right){\rho }^{\alpha }\left(\frac{{\mathbf{r}}_1+{\mathbf{r}}_2}{2}\right)\delta ({\mathbf{r}}_1-{\mathbf{r}}_2)$$

$$+i{W}_{LS}({\mathbf{\sigma }}_1+{\mathbf{\sigma }}_2)\cdot ({\overleftarrow{\nabla }}_1-{\overleftarrow{\nabla }}_2)\times \delta ({\mathbf{r}}_1-{\mathbf{r}}_2)({\overrightarrow{\nabla }}_1-{\overrightarrow{\nabla }}_2),$$

$${\Gamma }_{{\mathbf{n}}_1{\mathbf{n}}_3}=\sum _{{\mathbf{n}}_2{\mathbf{n}}_4}⟨{\mathbf{n}}_1{\mathbf{n}}_2|\hat{V}\hat{\mathcal{A}}|{\mathbf{n}}_3{\mathbf{n}}_4⟩{\rho }_{{\mathbf{n}}_4{\mathbf{n}}_2},{\Delta }_{{\mathbf{n}}_1{\mathbf{n}}_2}=\frac{1}{2}\sum _{{\mathbf{n}}_3{\mathbf{n}}_4}⟨{\mathbf{n}}_1{\mathbf{n}}_2|\hat{V}\hat{\mathcal{A}}|{\mathbf{n}}_3{\mathbf{n}}_4⟩{\kappa }_{{\mathbf{n}}_3{\mathbf{n}}_4}.$$

$${\tilde{\rho }}_{{\mathbf{n}}_3{\mathbf{n}}_4}=-2{\sigma }_{{\mathbf{n}}_4}{\kappa }_{{\mathbf{n}}_3{\bar{\mathbf{n}}}_4}^{\ast }.$$

$${\tilde{h}}_{{\mathbf{n}}_1{\mathbf{n}}_2}=-2{\sigma }_{{\mathbf{n}}_2}{\Delta }_{{\mathbf{n}}_1{\bar{\mathbf{n}}}_2}^{\ast }.$$

$${\tilde{h}}_{{\mathbf{n}}_1{\mathbf{n}}_2}=\frac{1}{2}\sum _{{\mathbf{n}}_3{\mathbf{n}}_4}⟨{\mathbf{n}}_1{\bar{\mathbf{n}}}_2|\hat{V}\hat{\mathcal{A}}|{\mathbf{n}}_3{\bar{\mathbf{n}}}_4⟩{\sigma }_{{\mathbf{n}}_2}{\sigma }_{{\mathbf{n}}_4}{\tilde{\rho }}_{{\mathbf{n}}_4{\mathbf{n}}_3}.$$

$$V_{{\mathbf{n}}_1{\mathbf{n}}_2{\mathbf{n}}_3{\mathbf{n}}_4}$$

$$=\int d^3{\mathbf{r}}_1\int d^3{\mathbf{r}}_2{\varphi }_{{\mathbf{n}}_1}^{\ast }({\mathbf{r}}_1){\varphi }_{{\mathbf{n}}_2}^{\ast }({\mathbf{r}}_2)e^{-({\mathbf{r}}_1-{\mathbf{r}}_2{)}^2/{\mu }^2}{\varphi }_{{\mathbf{n}}_3}({\mathbf{r}}_1){\varphi }_{{\mathbf{n}}_4}({\mathbf{r}}_2),$$

$$e^{-({\mathbf{r}}_1-{\mathbf{r}}_2)/{\mu }^2}=e^{-({\mathbf{\rho }}_1-{\mathbf{\rho }}_2)/{\mu }^2}{e}^{-(z_1-z_2)/{\mu }^2},{\varphi }_{\mathbf{n}}(\mathbf{r})={\varphi }_{n_r}^{\Lambda }(\rho ,\phi ){\varphi }_{n_z}(z).$$

$$V_{{\mathbf{n}}_1{\mathbf{n}}_2{\mathbf{n}}_3{\mathbf{n}}_4}=V_{n_r^{(1)}{\Lambda }^{(1)},n_r^{(2)}{\Lambda }^{(2)},n_r^{(3)}{\Lambda }^{(3)},n_r^{(4)}{\Lambda }^{(4)}}{V}_{n_z^{(1)}{n}_z^{(2)}{n}_z^{(3)}{n}_z^{(4)}}.$$

$$V_{n_z^{(1)}{n}_z^{(2)}{n}_z^{(3)}{n}_z^{(4)}}=\frac{\mu }{\sqrt{2{\pi }^3}{b}_z}\sum _{n=|n_z^{(2)}-n_z^{(4)}|,2}^{n_z^{(2)}+n_z^{(4)}}{T}_{n_z^{(2)}{n}_z^{(4)}}^n{\bar{F}}_{n_z^{(1)}{n}_z^{(3)}}^n,$$

$$T_{n_z^{(2)}{n}_z^{(4)}}^n\equiv \frac{\sqrt{n_z^{(2)}!n_z^{(4)}!n!}}{\left(\frac{n_z^{(2)}-n_z^{(4)}+n}{2}\right)!\left(\frac{n_z^{(4)}-n_z^{(2)}+n}{2}\right)!\left(\frac{n_z^{(2)}+n_z^{(4)}-n}{2}\right)!},$$

$$\bar{F}_{n_{z}^{(1)}n_{z}^{(3)}}^{n}\equiv\frac{\Gamma(\xi-n_{z}^{(1)})\Gamma(\xi-n_{z}^{(3)})\Gamma(\xi-n)}{\left(1+\frac{\mu^{2}}{2b_{z}^{2}}\right)^{\xi}\sqrt{n!n_{z}^{(1)}!n_{z}^{(3)}!}}{}_{2}F_{1}\Big{(}-n_{z}^{(1)},-n_{z}^{(3)};-\xi+n+1;-\frac{\mu^{2}}{2b_{z}^{2}}\Big{)},$$

$${\varphi }_{n_r}^{\Lambda }(\rho ,\phi )=\sum _{n_y=0}^{2n_r+|\Lambda |}{C}_{n_x{n}_y}^{n_r\Lambda }{\varphi }_{n_x}(x){\varphi }_{n_y}(y),$$

$$C_{n_x{n}_y}^{n_r\Lambda }\equiv \frac{2^{-n_r-|\Lambda |/2}(-1{)}^{n_r}\sqrt{n_r!(n_r+|\Lambda |)!}}{\sqrt{n_x!n_y!}{i}^{n_y}}\sum _{q=0}^{\min \left(n_y,n_r+\frac{|\Lambda |-\Lambda }{2}\right)}\left(\genfrac{}{}{0ex}{}{n_x}{n-q+\frac{|\Lambda |-\Lambda }{2}}\right)\left(\genfrac{}{}{0ex}{}{n_y}{q}\right)(-1{)}^{n_y-q}.$$

$$V_{n_r^{(1)}{\Lambda }^{(1)},n_r^{(2)}{\Lambda }^{(2)},n_r^{(3)}{\Lambda }^{(3)},n_r^{(4)}{\Lambda }^{(4)}}=$$

$$\times V_{n_x^{(1)}{n}_x^{(2)}{n}_x^{(3)}{n}_x^{(4)}}{V}_{n_y^{(1)}{n}_y^{(2)}{n}_y^{(3)}{n}_y^{(4)}}.$$

$${\hat{V}}_i\hat{\mathcal{A}}=V_i(\mathbf{r})\left(W_i+B_i{\hat{P}}_{\sigma }-H_i{\hat{P}}_{\tau }-M_i{\hat{P}}_{\sigma }{\hat{P}}_{\tau }\right)\left(1-{\hat{P}}_{\sigma }{\hat{P}}_{\tau }{\hat{P}}_r\right)\equiv {\hat{V}}_i^{\mathrm{ID}}-{\hat{V}}_i^{\mathrm{IE}}{\hat{P}}_{\tau },$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}=\sum _{{\mathbf{n}}_2{\mathbf{n}}_4}$$

$$-$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau }$$

$$\displaystyle\equiv\sum_{\bm{n}_{2}\bm{n}_{4}}\Big{[}\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau}|\bm{n}_{3}\bm{n}_{4}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau}+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}|\bm{n}_{3}\bm{n}_{4}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}}\Big{]}.$$

$${\hat{V}}_{i,\tau \tau }$$

$${\hat{V}}_{i,\tau \tau }^{\mathrm{SD}}$$

$${\hat{V}}_{i,\tau \tau }^{\mathrm{SE}}$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau }=\sum _{{\mathbf{n}}_2{\mathbf{n}}_4}$$

$$\displaystyle+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}\rangle\big{[}\langle\sigma_{\bm{n}_{1}}|\sigma\rangle\langle\sigma|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma\sigma}+\langle\sigma_{\bm{n}_{1}}|\sigma\rangle\langle\sigma^{\prime}|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma\sigma^{\prime}}$$

$$\displaystyle\phantom{+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}}+\langle\sigma_{\bm{n}_{1}}|\sigma^{\prime}\rangle\langle\sigma|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma^{\prime}\sigma}+\langle\sigma_{\bm{n}_{1}}|\sigma^{\prime}\rangle\langle\sigma^{\prime}|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma^{\prime}\sigma^{\prime}}\big{]}$$

$$\displaystyle+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}\rangle\big{[}\langle\sigma_{\bm{n}_{1}}|\sigma\rangle\langle\sigma|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma\sigma}+\langle\sigma_{\bm{n}_{1}}|\sigma\rangle\langle\sigma^{\prime}|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma\sigma^{\prime}}$$

$$\displaystyle\phantom{+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}}+\langle\sigma_{\bm{n}_{1}}|\sigma^{\prime}\rangle\langle\sigma|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma^{\prime}\sigma}+\langle\sigma_{\bm{n}_{1}}|\sigma^{\prime}\rangle\langle\sigma^{\prime}|\sigma_{\bm{n}_{3}}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma^{\prime}\sigma^{\prime}}\big{]}\Big{\}}.$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau \sigma \sigma }$$

$$\displaystyle\quad\quad+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}+\hat{V}_{i,\tau\tau^{\prime}}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma\sigma}+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}|\bm{n}_{3}\bm{n}_{4}\rangle\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma^{\prime}\sigma^{\prime}}\Big{]}\vskip 6.0pt plus 2.0pt minus 2.0pt$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau \sigma {\sigma }^{\prime }}$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau \sigma \sigma }$$

$$\displaystyle\quad\quad+\big{(}\langle\bm{n}_{1}\bar{\bm{n}}_{2}|\hat{V}_{i,\tau\tau}^{\rm SD}+\hat{V}_{i,\tau\tau}^{\rm SE}|\bm{n}_{3}\bar{\bm{n}}_{4}\rangle+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau}^{\rm SD}|\bm{n}_{3}\bm{n}_{4}\rangle\big{)}\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma^{\prime}\sigma^{\prime}}$$

$$\displaystyle\quad\quad+\big{(}\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}+\hat{V}_{i,\tau\tau^{\prime}}^{\rm SE}|\bm{n}_{3}\bm{n}_{4}\rangle+\langle\bm{n}_{1}\bar{\bm{n}}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}|\bm{n}_{3}\bar{\bm{n}}_{4}\rangle\big{)}\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma\sigma}$$

$$\displaystyle\quad\quad+\big{(}\langle\bm{n}_{1}\bar{\bm{n}}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}+\hat{V}_{i,\tau\tau^{\prime}}^{\rm SE}|\bm{n}_{3}\bar{\bm{n}}_{4}\rangle+\langle\bm{n}_{1}\bm{n}_{2}|\hat{V}_{i,\tau\tau^{\prime}}^{\rm SD}|\bm{n}_{3}\bm{n}_{4}\rangle)\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau^{\prime}\sigma^{\prime}\sigma^{\prime}}\Big{]},$$

$${\Gamma }_{i,{\mathbf{n}}_1{\mathbf{n}}_3}^{\tau \sigma {\sigma }^{\prime }}$$

Full text

**Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis

(III) hfbtho (v3.00): a new version of the program. **

R. Navarro Perez,a N. Schunck,a111E-mail: schunck1@llnl.gov R.-D. Lasseri,b C. Zhang,c J. Sarich,d

*a**Nuclear and Chemical Science Division, Lawrence Livermore National Laboratory, Livermore, CA 94551, USA

bInstitut de Physique Nucléaire, IN2P3-CNRS, Université Paris-Sud, Université Paris-Saclay, F-91406 Orsay Cedex, France

cNSCL/FRIB Laboratory, Michigan State University, East Lansing, MI 48824, USA

dMathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA

Abstract

We describe the new version 3.00 of the code hfbtho that solves the nuclear Hartree-Fock (HF) or Hartree-Fock-Bogolyubov (HFB) problem by using the cylindrical transformed deformed harmonic oscillator basis. In the new version, we have implemented the following features: (i) the full Gogny force in both particle-hole and particle-particle channels, (ii) the calculation of the nuclear collective inertia at the perturbative cranking approximation, (iii) the calculation of fission fragment charge, mass and deformations based on the determination of the neck (iv) the regularization of zero-range pairing forces (v) the calculation of localization functions (vi) a MPI interface for large-scale mass table calculations.

PACS numbers: 07.05.T, 21.60.-n, 21.60.Jz

NEW VERSION PROGRAM SUMMARY

Title of the program: hfbtho v3.00

Catalogue number: ….

Program obtainable from: CPC Program Library,  Queen’s University of Belfast, N. Ireland (see application form in this issue)

Reference in CPC for earlier version of program: M.V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, and S. Wild, Comput. Phys. Commun. 184 (2013).

Catalogue number of previous version: ADUI_v2_0

Licensing provisions: GPL v3

Does the new version supersede the previous one: Yes

Computers on which the program has been tested: Intel Xeon, AMD-Athlon, AMD-Opteron, Cray XT5, Cray XE6

Operating systems: UNIX, LINUX

Programming language used: FORTRAN-95

Memory required to execute with typical data: 200 Mwords

No. of bits in a word: 8

Has the code been vectorised?: Yes

Has the code been parallelized?: Yes

No. of lines in distributed program: 20 458

Keywords: Nuclear many-body problem; Density functional theory; Energy density functional theory; Self-consistent mean field; Hartree-Fock-Bogolyubov theory; Finite-temperature HArtree-Fock-Bogolyubov theory; Skyrme interaction; Gogny force; Pairing correlations; Particle number projection; Pairing regularization; Collective inertia; Nuclear matter; Constrained calculations; Potential energy surface; Harmonic oscillator; Transformed harmonic oscillator; Shared memory parallelism; Distributed memory parallelism; Tensor contractions; Loop optimization.

Nature of physical problem

hfbtho is a physics computer code that is used to model the structure of the nucleus. It is an implementation of the energy density functional (EDF) approach to atomic nuclei, where the energy of the nucleus is obtained by integration over space of some phenomenological energy density, which is itself a functional of the neutron and proton intrinsic densities. In the present version of hfbtho, the energy density derives either from the zero-range Skyrme or the finite-range Gogny effective two-body interaction between nucleons. Nuclear superfluidity is treated at the Hartree-Fock-Bogoliubov (HFB) approximation. Constraints on the nuclear shape allows probing the potential energy surface of the nucleus as needed e.g., for the description of shape isomers or fission. The implementation of a local scale transformation of the single-particle basis in which the HFB solutions are expanded provide a tool to properly compute the structure of weakly-bound nuclei.

Method of solution

The program uses the axial Transformed Harmonic Oscillator (THO) single-particle basis to expand quasiparticle wave functions. It iteratively diagonalizes the Hartree-Fock-Bogolyubov Hamiltonian based on generalized Skyrme-like energy densities and zero-range pairing interactions or the finite-range Gogny force until a self-consistent solution is found. A previous version of the program was presented in M.V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, and S. Wild, Comput. Phys. Commun. 184 (2013) 1592–1604 with much of the formalism presented in the original paper M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, P. Ring, Comput. Phys. Commun. 167 (2005) 43–63.

Summary of revisions

the Gogny force in both particle-hole and particle-particle channels was implemented; 2. 2.

the nuclear collective inertia at the perturbative cranking approximation was implemented; 3. 3.

fission fragment charge, mass and deformations were implemented based on the determination of the position of the neck between nascent fragments; 4. 4.

the regularization method of zero-range pairing forces was implemented; 5. 5.

the localization functions of the HFB solution were implemented; 6. 6.

a MPI interface for large-scale mass table calculations was implemented.

Restrictions on the complexity of the problem

Axial- and time-reversal symmetries are assumed.

Typical running time

Highly variable, as it depends on the nucleus, size of the basis, requested accuracy, requested configuration, compiler and libraries, and hardware architecture. An order of magnitude would be a few seconds for ground-state configurations in small bases $N_{\text{max}}\approx 8-12$, to a few minutes in very deformed configuration of a heavy nucleus with a large basis $N_{\text{max}}>20$.

Unusual features of the program

The user must have access to (i) the LAPACK subroutines dsyeevr, dsyevd, dsytrf and dsytri, and their dependencies, which compute eigenvalues and eigenfunctions of real symmetric matrices, (ii) the LAPACK subroutines dgetri and dgetrf, which invert arbitrary real matrices, and (iii) the BLAS routines dcopy, dscal, dgemm and dgemv for double-precision linear algebra (or provide another set of subroutines that can perform such tasks). The BLAS and LAPACK subroutines can be obtained from the Netlib Repository at the University of Tennessee, Knoxville: http://netlib2.cs.utk.edu/.

LONG WRITE-UP

1 Introduction

The program hfbtho solves the Hartree-Fock-Bogolyubov (HFB) equation for atomic nuclei by expanding the solution on the harmonic oscillator basis. The original version 1.66 of the program published in [1] had the following characteristics: the HFB equation was only solved for even-even nuclei with Skyrme two-body effective potentials in the particle-hole channel and contact pairing interactions; the code implemented the transformed harmonic oscillator basis and particle number projection; the only available constraint was on the axial quadrupole deformation and parity-conservation of the HFB solution was explicitly enforced. The release 2.00d published in [2] removed some of these restrictions and added several new capabilities: parity-breaking of the HFB solution was implemented, constraints on axial multipole moments from $\lambda=2$ to $\lambda=9$ became available and the method to handle these constraints was based on the readjustment of the Lagrange parameters based on the cranking approximation of the QRPA matrix; odd nuclei could be computed based on the equal filling approximation of the blocking prescription, and the code also implemented the finite-temperature extension of the HFB equation.

This release introduces several major features in the code. Most notably, version 3.00 of hfbtho can now solve the HFB equation for the finite-range Gogny force, in both the particle-hole and particle-particle channel. The new version also allows the computation of the collective inertia tensor and zero-point energy correction within the generator coordinate method (GCM) and adiabatic time-dependent Hartree-Fock-Bogolyubov (ATDHFB) approximations. It also contains a small fission toolkit that gives the charge, mass and deformations of the fission fragments based on position of the neck between pre-fragments – defined as the expectation value of the Gaussian neck operator. Finally, the code also implements the regularization of the pairing interaction and the calculation of localization functions.

In section 2, we review the modifications introduced in this version. In section 3, we give a few numerical benchmarks for the new capabilities. Finally, in section 4, we discuss the new options available in the input file and explain how to run the code, in particular in parallel mode.

2 Modifications introduced in version 3.00

We present in this section the major new features added to the code between version 2.00d and 3.00.

2.1 Gogny force

Version 3.00 of hfbtho implements energy functionals derived from the Gogny two-body effective interaction [3],

[TABLE]

where $\hat{P}_{\sigma}$ and $\hat{P}_{\tau}$ are the spin and isospin exchange operators. The D1, D1S, D1’ and D1N parametrizations of the Gogny force are implemented in this version; see  [3, 4, 5] for information about each of these parametrizations. Calculations can be performed using any of these parametrizations by setting the functional name in the hfbtho_NAMELIST.dat file to D1, D1S, D1p or D1N.

In previous versions of hfbtho the HFB fields in coordinate space were directly integrated using Gaussian quadrature at each iteration to calculate the HFB matrix. Given the contact nature of the Skyrme-like functionals that had been implemented so far, the precise integration of the HFB fields requires only a few dozens of points. The same level of precision for the finite range part in the Gogny functional would require several hundreds of points, which would significantly increase the computation time at each iteration. To avoid this bottleneck the finite-range contributions to the HFB fields are computed in configuration space, that is, by performing the following contractions of the matrix elements of the anti-symmetrized potential with either the density matrix $\rho$ or the pairing tensor $\kappa$

[TABLE]

Note that in these expressions, summations extend over the complete basis of the single-particle Hilbert space, that is, the basis made of the states $(|\bm{n}\rangle,|\bar{\bm{n}}\rangle)$, where $|\bar{\bm{n}}\rangle=\hat{T}|\bm{n}\rangle$. In this basis, the matrix of the pairing tensor obeys $\kappa_{\bm{n}_{3}\bm{n}_{4}}=\kappa_{\bar{\bm{n}}_{3}\bar{\bm{n}}_{4}}=0$, as a consequence of the particular structure of the Bogoliubov matrices; see section 2.2. In practice, hfbtho uses what is known as the russian convention in the pairing channel, that is, it works with the pairing density $\tilde{\rho}$ instead of the pairing tensor $\kappa$ [6]. In configuration space, the two are related through

[TABLE]

As a consequence, it is more convenient to also use the Hermitian pairing field $\tilde{h}$ instead of the antisymmetric pairing field $\Delta$,

[TABLE]

An additional simplification specific to built-in time-reversal symmetry in hfbtho is that the pairing tensor, density and pairing fields are real and symmetric [7]. Therefore, we will use the following formula in the pairing channel,

[TABLE]

In these expressions, the basis states $|\bm{n}\rangle$ are defined by the quantum numbers of the HO, $\bm{n}\equiv(n_{r},\Lambda,n_{z},\sigma_{\bm{n}},\tau_{\bm{n}})$ with the spin projection $\Sigma_{\bm{n}}=\sigma_{\bm{n}}/2=\pm 1/2$ and isospin $\tau_{\bm{n}}$ quantum numbers; $\hat{\mathcal{A}}=1-\hat{P}_{r}\hat{P}_{\sigma}\hat{P}_{\tau}$ is the anti-symmetrization operator; under time-reversal operation, $\Lambda$ and $\sigma_{\bm{n}}$ change sign, hence, $\bm{n}\equiv(n_{r},\Lambda,n_{z},\sigma_{\bm{n}},\tau_{\bm{n}})\Rightarrow\bar{\bm{n}}\equiv(n_{r},-\Lambda,n_{z},-\sigma_{\bm{n}},\tau_{\bm{n}})$. Finally, note that the matrix elements of the two-body potential are independent of the HFB iterations. This allows computing them only once at the beginning of the iterative loop.

2.1.1 Matrix element of a two-body Gaussian potential

The calculation of the finite-range part of the full two-body matrix element of potential (1) requires evaluating spatial integrals of the type

[TABLE]

where $\phi_{\bm{n}}(\bm{r})$ is the spatial part of the stretched harmonic oscillator basis function in the cylindrical coordinates $\bm{r}\equiv(\rho,z,\varphi)$. A convenient property of the two-body Gaussian potential and the HO wave function is that each of them can be separated into a product of axial and polar functions, namely

[TABLE]

This allows separating the matrix element (7) into

[TABLE]

Our implementation to calculate these matrix elements closely follows the method outlined in [8]. For the axial matrix element, we use the expansion

[TABLE]

which preserves numerical precision at large number of harmonic oscillator shells. The coefficient $T_{n_{z}^{(2)}n_{z}^{(4)}}^{n}$ is defined as

[TABLE]

and the function $\bar{F}_{n_{z}^{(1)}n_{z}^{(3)}}^{n}$ as

[TABLE]

where $\xi=(n_{z}^{(1)}+n_{z}^{(3)}+n+1)/2$, ${}_{2}F_{1}$ is the hyper-geometric function and $\Gamma$ is the usual gamma function.

Since there is no analogous expansion to (10) for the polar matrix element, we expand the two-dimensional HO wave function into a sum of products of two one-dimensional HO functions (equivalent to a change from polar to Cartesian coordinates),

[TABLE]

where $n_{x}=2n_{r}+|\Lambda|-n_{y}$ and

[TABLE]

This expansion transforms the polar matrix elements into a sum of products of axial matrix elements, namely

[TABLE]

Each axial matrix elements in (15) can be calculated accurately with (10) by making the substitution $b_{z}\rightarrow b_{\bot}$. Reference [8] includes other simpler and more direct expansions to calculate $V_{n_{r}^{(1)}\Lambda^{(1)},n_{r}^{(2)}\Lambda^{(2)},n_{r}^{(3)}\Lambda^{(3)},n_{r}^{(4)}\Lambda^{(4)}}$ and $V_{n_{z}^{(1)},n_{z}^{(2)},n_{z}^{(3)},n_{z}^{(4)}}$. These expansions were implemented for testing and debugging purposes only. The corresponding subroutines and functions are included in version 3.00but are not used, since these expansions loose accuracy as the basis size increases. For a complete list of subroutines and functions see the documentation for the hfbtho_gogny module.

2.1.2 Contraction with the two-body matrix elements

Once the axial and polar matrix elements have been calculated with (10) and (15) before the first HFB iteration, a contraction with the one body density and pairing matrices has to be made at each iteration.

**Mean field - ** For the calculation of the mean field $\Gamma_{\bm{n}_{1}\bm{n}_{3}}$ it is more convenient to apply the anti-symmetrization operator to the left and embed it into the potential, rather than to the right through its action on the two-body ket. This results in

[TABLE]

where the ID and IE labels stand for isospin direct and isospin exchange respectively.

In this subsection we slightly abuse the notation so that the $\bm{n}_{i}$ label includes all the quantum numbers that are not written explicitly or have not been contracted. Inserting (16) into the expression for $\Gamma$ in (2), applying the isospin exchange operator to the right and contracting the $\tau$ indices (i.e. summing over all possible values of $\tau_{\bm{n}_{2}}$ and $\tau_{\bm{n}_{4}}$) results in

[TABLE]

where $\tau\neq\tau^{\prime}$ and the orthonormality of the isospin states along with the fact that $\rho$ is block-diagonal in isospin (i.e. $\rho_{\bm{n}_{2}\bm{n}_{4}}^{\tau\tau^{\prime}}=0$) has been used. From (17) one can quickly see that $\Gamma$ is also block-diagonal in isospin, namely $\Gamma_{i,\bm{n}_{1}\bm{n}_{3}}^{\tau_{\bm{n}_{1}}\tau_{\bm{n}_{3}}}=0$ if $\tau_{\bm{n}_{1}}\neq\tau_{\bm{n}_{3}}$, and removing the second $\tau$ for simplicity we obtain

[TABLE]

Separating $\hat{V}_{i,\tau\tau}$ and $\hat{V}_{i,\tau\tau^{\prime}}$ into spin-direct and spin-exchange parts

[TABLE]

allows applying the spin exchange operator to the right and contracting the $\sigma$ indices to obtain

[TABLE]

This expression gives different results for states with $\sigma_{\bm{n}_{1}}=\sigma_{\bm{n}_{3}}$, denoted as $\sigma\sigma$, and states with $\sigma_{\bm{n}_{1}}\neq\sigma_{\bm{n}_{3}}$, denoted as $\sigma\sigma^{\prime}$

[TABLE]

At this point it is important to remember that since axial symmetry is explicitly enforced in hfbtho, $\Omega=\Lambda+\sigma/2$ is a good quantum number and therefore the $\Gamma$, $\tilde{h}$, $\rho$ and $\tilde{\rho}$ matrices are block-diagonal. Each block is characterized by a different value of $\Omega$. Furthermore, due to time reversal symmetry, hfbtho works only with positive values of $\Omega$. However, the summations in (2) include both positive and negative values of $\Lambda$ and $\sigma$. Therefore we apply the time reversal operator $\hat{T}|\bm{n}\rangle=|\bar{\bm{n}}\rangle$ to blocks with $\Omega<0$ in order to express such states in terms of states with $\Omega>0$. This results in

[TABLE]

Inserting the definitions (20), it is possible to separate the direct and exchange contributions to $\Gamma$ due to the presence of $\hat{P}_{r}$ so that $\Gamma_{\bm{n}_{1}\bm{n}_{3}}=\Gamma_{\bm{n}_{1}\bm{n}_{3}}^{\rm D}+\Gamma_{\bm{n}_{1}\bm{n}_{3}}^{\rm E}$. Strictly speaking this separation is not necessary to calculate and diagonalize the HFB Hamiltonian, but it provides a convenient tool to verify our implementation against other codes. Using the notation from (7) we find the final result

[TABLE]

where the $\Omega>0$ label in the summations is implicit. Finally, we note that $V_{\bm{n}_{1}\bm{n}_{2}\bm{n}_{3}\bm{n}_{4}}=V_{\bm{n}_{1}\bar{\bm{n}}_{2}\bm{n}_{3}\bar{\bm{n}}_{4}}$ along with $\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma\sigma^{\prime}}=\rho_{\bm{n}_{4}\bm{n}_{2}}^{\tau\sigma^{\prime}\sigma}$ results in $\Gamma_{i,\bm{n}_{1}\bm{n}_{3}}^{{\rm D}\tau\sigma\sigma^{\prime}}=0$.

**Pairing field $\bm{\tilde{h}}$ - ** For the contraction of the two-body matrix elements with the pairing density $\tilde{\rho}$, we take advantage of the fact that for local two-body potentials, the exchange and direct contributions to the pairing field are identical. This cancels out the $1/2$ factor and only the direct part needs to be calculated, i.e.,

[TABLE]

As with $\Gamma$, after contracting the $\tau$ indices $\tilde{h}$ turns out to be block diagonal in isospin since $\tilde{\rho}$ is also block diagonal in isospin and the states $|\tau\rangle$ are orthonormal,

[TABLE]

Now, the contraction with respect to $\sigma$ can be performed explicitly. Again, this results in different expressions for spin diagonal states ($\sigma_{\bm{n}_{1}}=\sigma_{\bm{n}_{2}}$) and spin off-diagonal states ($\sigma_{\bm{n}_{1}}\neq\sigma_{\bm{n}_{2}}$),

[TABLE]

Finally, as before, the summations run over both positive and negative values of $\Omega$. Therefore, transformation properties of basis states under time reversal are again used to express the matrix elements of $\tilde{\rho}$ with negative $\Omega$ in terms of matrix elements with positive $\Omega$. This results in

[TABLE]

where the summations now run only over states with $\Omega>0$.

**Loop optimization - ** The separability of the matrix elements in axial and polar components is exploited during the calculation of the HFB fields. The contractions (2.1.2)-(2.1.2) and (2.1.2)-(35) correspond to contractions of the type

[TABLE]

where the index $\bm{n}_{k}$ contains only the quantum numbers from the axially-symmetric harmonic oscillator basis. Naively performing this summation with $N$ harmonic oscillator shells to calculate all $\Gamma_{\bm{n}_{1}\bm{n}_{3}}$ matrix elements results in an nested loop in which the number of operations is of the order of $N^{12}$. The number of operations can be significantly reduced by taking advantage of the separability of $V_{\bm{n}_{1}\bm{n}_{2}\bm{n}_{3}\bm{n}_{4}}$ into axial and polar components. For example, the mean field can be computed as

[TABLE]

The calculation of all the $\mathsf{Z}$ objects results in a nested loop where the number of operations increases like $N^{8}$. The subsequent contraction of the polar component of the two-body matrix element with $\mathsf{Z}$ requires a loop of order $N^{10}$. This simple rearrangement in the calculation of $\Gamma$ and $\tilde{h}$ with $16$ shells reduces the calculation time at each iteration by $97.7\%$.

2.2 Collective inertia

The code hfbtho in version 3.00 calculates the collective inertia tensor using either the generator coordinate method (GCM) or the adiabatic time-dependent Hartree-Fock-Bogolyubov (ATDHFB) approximation. In the case of the GCM, the calculation is performed using the local approximation of the Gaussian overlap approximation (GOA), while in the case of the ATDHFB theory, we use the perturbative cranking approximation; We give below only a short description of what is implemented in the code; see [9] for an recent overview of both methods and Chapter 10.7 in [10] for a textbook presentation.

**GCM-GOA - ** Recall that under the GOA, the many-body Schrödinger equation can be recast into a collective Schrödinger-like equation of the following form,

[TABLE]

where we note $\bm{q}=(q_{1},\dots,q_{N})$ a vector of collective variables (for instance, each of the $q_{i}$ would stand for a multipole moment operator). The collective inertia tensor is $\mathsf{B}\equiv B_{kl}(\bm{q})$, and $V_{\mathrm{coll}}(\bm{q})$ is the collective potential. It is given by

[TABLE]

with $V(\bm{q})$ the HFB energy for the values $\bm{q}$ of the collective variables and $\epsilon_{\mathrm{ZPE}}$ a zero-point energy correction accounting for the quantum fluctuations of the collective variables. Neglecting derivative terms, we have

[TABLE]

where $\mathsf{\Gamma}$ is the GCM metric. The term $\mathsf{\Gamma}\mathsf{B}$ is a contraction of two rank-2 tensors, that is, $\mathsf{\Gamma}\mathsf{B}\equiv\sum_{ab}\Gamma_{ab}B_{ba}$.

In the HFB theory with constraints on the collective variables $\bm{q}$, one can write the action of the collective momentum $\bm{P}\equiv-i/\hbar\partial/\partial\bm{q}$ in terms of quasiparticle creation and annihilation operators by using the Ring and Schuck theorem [11]. By invoking linear response theory, the matrix elements of the collective momentum can also be expressed as a function of the matrix of the collective variables. This provides a closed-form for the GCM metric, the collective inertia tensor and the zero-point energy correction, which involve only the characteristics of the HFB solution at the point $\bm{q}$ and the QRPA matrix at the same point. Using the cranking approximation for the QRPA matrix provides simple expressions for all these terms. The metric becomes

[TABLE]

with the moments $\mathsf{M}^{(K)}$ being

[TABLE]

where $E_{\mu}$ is the quasiparticle energy of quasiparticle $\mu$ and $|\mu\nu\rangle=\beta^{\dagger}_{\mu}\beta^{\dagger}_{\nu}|\Phi\rangle$ is a two-quasiparticle excitation. The term $\langle\Phi|\hat{Q}^{\dagger}_{\alpha}|\mu\nu\rangle$ corresponds to the block matrix $\tilde{Q}_{\alpha}^{12}$ of the matrix of the operator $\hat{Q}_{\alpha}$ in the quasiparticle basis,

[TABLE]

The collective mass tensor $\mathsf{M}=\mathsf{B}^{-1}$ reads

[TABLE]

and the zero-point energy can be written

[TABLE]

**ATDHFB - ** The ATDHF theory is presented in details in [12]. A quick summary of its extension to include pairing correlations is recalled in [9]. The approximation relies on writing the full time-dependent density of the TDHFB theory as the result of a unitary transformation $e^{-i\hat{\chi}(t)}$ of a time-dependent, time-even generalized density $\mathcal{R}_{0}(t)$, which plays the role of generalized coordinate for the collective motion of the nucleus. By subsequently expanding the TDHFB energy up to second order in $\hat{\chi}$, one obtains a collective, Schrödinger-like equation of motion. Further assuming that the time dependence of the generalized density $\mathcal{R}_{0}(t)$ is constrained in the collective space spanned by the variables $\bm{q}=(q_{1},\dots,q_{N})$, such that

[TABLE]

the collective equation takes nearly the same form as (39). The main difference is the absence of a zero-point energy correction from the collective potential (since the ATDHFB theory is a classical approximation). In addition, the expression for the inertia tensor $\mathsf{B}$ is slightly different from the GCM one. Using both the cranking approximation of the QRPA matrix and writing the derivatives $\partial\mathcal{R}_{0}/\partial q_{\alpha}$ in terms of the matrix of a collective momentum operator, we arrive at

[TABLE]

Although conceptually not very well-founded, it is customary to compute a zero-point energy correction for the ATDHFB energy by using (46) with the mass tensor (48).

**Implementation in hfbtho - ** The key element of the implementation is the calculation of the matrix $\tilde{Q}_{\alpha}^{12}$. Let us recall that the matrix of a one-body operator in the doubled single-particle basis reads

[TABLE]

The transformation into the quasiparticle basis to obtain $\tilde{Q}_{\alpha}$ is performed with the Bogolyubov transformation

[TABLE]

and yields

[TABLE]

At this point, we must consider specifically the structure of the single particle basis. The full basis should be complete under time-reversal symmetry, that is, for each basis state $|\bm{n}\rangle$, the state $|\bar{\bm{n}}\rangle$ is also a basis state. In hfbtho, recall that we have $\bm{n}\equiv(n_{r},\Lambda,n_{z},\Sigma)$ and

[TABLE]

with $\Sigma=\pm 1/2$, $\sigma=2\Sigma$ and $(\rho,z,\varphi)$ the cylindrical coordinates; see [1] for the full definition of each wave function. The time-reversed wave function then reads

[TABLE]

Recall that the HFB spinor for quasiparticle $\mu$ read

[TABLE]

Owing to the fact that time-reversal symmetry is conserved in hfbtho, and that the functions $U^{\pm}$ and $V^{\pm}$ are real, the spinor

[TABLE]

is also a solution of the HFB equation with the same eigenvalue $E_{\mu}$. These various relations determine the structure of the matrices of $U$ and $V$ in the full s.p. basis. Indeed, we have

[TABLE]

We can then easily show that $U_{\bm{n}\bar{\mu}}=U_{\bar{\bm{n}}\mu}=0$ and $U_{\bm{n}\mu}=U_{\bar{\bm{n}}\bar{\mu}}$, and $V_{\bm{n}\mu}=V_{\bar{\bm{n}}\bar{\mu}}=0$ and $V_{\bm{n}\bar{\mu}}=-V_{\bar{\bm{n}}\mu}$. Therefore, the matrices have the following generic structure

[TABLE]

Taking into account that the matrices $U$ and $V$ are real and that the matrix $Q_{\alpha}$ has the same generic structure as $U$, it then follows that

[TABLE]

In practice, the moments (43) are thus computed through the formula

[TABLE]

where the factor 2 comes from the particular structure of the $\tilde{Q}_{\alpha}^{12}$ matrix.

2.3 Fission toolkit

We imported in version 3.00 of the code hfbtho a few routines implemented in version 2.73 of hfodd to compute the characteristics of fission fragments. These include

• •

the capability to set the number of particles in the neck by using as a constraint the expectation value of the Gaussian neck operator

[TABLE]

where $z_{N}$ gives the position of the neck (along the $z-$axis of the intrinsic reference frame) between the two nascent fragments. As in hfodd, it is defined as the point near the origin of the intrinsic reference frame where the density is the lowest. The range $a_{N}$ gives the spatial extent of the neck and is fixed at $a_{N}=1$ fm in the code.

• •

the charge and mass of each fragment. Expectation values of the multipole moments are also computed, both with respect to the center of mass of the fissioning nucleus and that of each of the individual fragments.

The various routines are coded in the module hfbtho_fission.f90. The characteristics of the fission fragments are computed at the convergence of the HFB iterations by setting the parameter fission_fragments = T in the namelist HFBTHO_FEATURES. A new namelist for constraints on the size of the neck called HFBTHO_NECK is available; see section 4.1 for a description.

We recall that the Gaussian neck operator (60) is purely spatial and does not depend on spin or isospin degrees of freedom. Its expectation value is computed in coordinate space on the Gauss-Hermite quadrature grid used for the $z$-coordinate in hfbtho.

[TABLE]

where $\rho(\bm{r})$ is the isoscalar density. The integral is computed by Gauss quadrature. Since hfbtho works in cylindrical coordinates $(\rho,z,\varphi)$ and uses the stretched coordinates $\xi_{z}=b_{z}z$ and $\eta=b_{\perp}^{2}\rho^{2}$. The code tabulates the value of the Gaussian neck operator on the Gauss-Laguerre and Gauss-Hermite meshes, $Q_{N}(\xi_{h},\eta_{\ell})$. The expectation value is then given by

[TABLE]

with $\omega_{h}$, $\omega_{\ell}$ and $\xi_{h},\eta_{\ell}$ the weights and nodes of the Gauss-Hermite and Gauss-Laguerre quadrature, respectively.

The expectation value of the neck operator can be set by imposing a linear constraint in the HFB equation. The Lagrange parameter controlling this constraint is readjusted based on the cranking approximation of the QRPA matrix; see section 2.4 in [2] and section 2.3 in [13] for details about the method. This technique requires the matrix elements of the operator in the HO basis. In hfbtho, the matrix elements are computed directly by first defining the value of the Gaussian neck operator on the Gauss-Laguerre and Gauss-Hermite meshes and then computing the matrix elements in each $\Omega$ block as

[TABLE]

where $\phi_{\bm{n}}^{*}(\bm{r},\sigma)$ are the basis spinor (52); see [1] for additional definitions.

Given the position of the neck, we calculate the mass of the two fission fragment classically, that is, by integrating the density up to the neck position,

[TABLE]

A similar formula applies to the charge of the fragments. Since the integral over $z$ does not extend to $+\infty$, Gauss-Hermite quadrature cannot be used. In our implementation, we interpolate the value of the density along the $z$-axis using spline functions, and compute the integral over $z$ using the Simpsons’s rule on a mesh of 1000 points. Although more accurate integration schemes could be implemented, our choice is numerically accurate enough to within the 4th digit of the value of $A_{1}$.

To compute the deformations of the fission fragments in their intrinsic frame, we first compute the positions of the centers of mass of each fragment,

[TABLE]

Recall that the multipole moments operators are defined in spherical coordinates by

[TABLE]

where $Y_{\lambda\mu}(\theta,\varphi)$ are the spherical harmonics. Since hfbtho assumes axial symmetry, only $\mu=0$ multipole moments are non-zero. It is easy to show that they can be expressed in terms of $\rho$ and $z$ only. The expectation value of the multipole moment of order $\lambda$ in the intrinsic frame of the fission fragment $1$ is thus given by

[TABLE]

The expression for the second fragment is obtained by integrating over $z$ from $z_{N}$ to $+\infty$. As with the charge and mass numbers of the fission fragments, we choose to interpolate the values of the density along the $z$-axis with Spline functions and carry out all integrals over $z$ with the Simpson’s 3/8 rule.

2.4 Regularization of zero-range pairing force

The code hfbtho now implements the regularization of zero-range pairing forces according to the original idea introduced by Bulgac in [14]. Our implementation follows most closely [14, 15]. It is limited to the case of functionals of the local pairing density $\tilde{\rho}(\bm{r})$ only. Moreover, following arguments given in [15], and explicit benchmarks of [16], the Hartree-Fock field resulting from the variations of the Fermi momentum $k_{F}$ and cut-off momentum $k_{c}$ with respect to the density are not implemented in version 3.00.

Under these conditions, the pairing energy is written as an integral of a functional of the local pairing density,

[TABLE]

and the local pairing field is simply

[TABLE]

The regularization procedure consists in computing an effective pairing strength function $g_{\text{eff}}(\bm{r})$ such that these expressions simply become

[TABLE]

For each type of particles $q$ (protons or neutrons), the position-dependent effective pairing strength $g_{\text{eff}}(\bm{r})$ is given by the table below, based on the real or imaginary nature of the two cut-off momenta $k_{\text{cut}}(\bm{r})$ and $\ell_{\text{cut}}(\bm{r})$ and the Fermi momentum $k_{F}(\bm{r})$:

[TABLE]

In these expressions, the scalar effective mass for particle $q$ is defined as

[TABLE]

where $C_{t}^{\rho\tau}$, $t=0,1$ are the isoscalar and isovector coupling constants of the kinetic density part of the Skyrme functional, $\rho_{0}$ is the isoscalar density and $\rho_{q}$ the density for particle $q$. The Fermi momentum $k_{F}(\bm{r})$ is given by

[TABLE]

where $\mu$ is the Fermi energy and $U(\bm{r})$ is the self-consistent mean field potential; cf. Eq.(18) in [1]. The cutoff momenta read

[TABLE]

with $\epsilon_{\text{cut}}$ the cut-off in quasiparticle energies (relative to the Fermi energy), $\epsilon_{\text{cut}}=E_{\text{max}}-\mu$. Recall that $E_{\text{max}}$ is an input parameter of hfbtho defined in the namelist under the keyword pairing_cutoff.

2.5 Localization functions

The code hfbtho in version 3.00 calculates the nucleon spatial localization with fixed spin and isospin in subroutine hfbtho_localization.f90. The spatial localization was originally introduced in atomic and molecular physics to characterize chemical bond structures in electron systems [17, 18, 19, 20, 21, 22], and then was applied to identify cluster structures in light nuclei [23] and fragment formation in fissioning heavy nuclei [24].

The definition of spatial localization is based on the probability of finding a nucleon within a distance $\delta$ from a given nucleon at position $\bm{r}$ with the same spin $\sigma$ ($=\uparrow$ or $\downarrow$) and isospin $q$ ($=n$ or $p$). By applying a Taylor expansion with respect to $\delta$, this probability can be written as

[TABLE]

where $\rho_{q\sigma}$, $\tau_{q\sigma}$, $\bm{j}_{q\sigma}$, and $\bm{\nabla}\rho_{q\sigma}$ are the particle density, kinetic energy density, current density, and density gradient, respectively. Through the canonical HFB orbits $\psi_{\alpha}$, they can be expressed as

[TABLE]

with $v^{2}$ being the canonical occupation probability. Thus, the expression in the parentheses of (75) can serve as a localization measure. Since the spatial localization and $R_{q\sigma}$ are in a reverse relationship, we can introduce a dimensionless expression for the spatial localization by normalizing with the Thomas-Fermi kinetic energy density $\tau_{q\sigma}=\frac{3}{5}(6\pi^{2})^{2/3}\rho_{q\sigma}^{5/3}$,

[TABLE]

In hfbtho, time reversal symmetry is always conserved, therefore, $\bm{j}_{q\sigma}$ vanishes.

2.6 MPI capabilities

Two new large scale production modes have been implemented in version 3.00, DO_MASSTABLE and DRIP_LINES. To activate any of these two new modes, the corresponding pre-processor directive must be set to 1 in the makefile before compilation. To activate MPI for any of these two modes the pre-processor directive USE_MPI must be set to 1.

The DO_MASSTABLE mode allows performing, serially or in parallel, a list of hfbtho calculations with different number of protons, number of neutrons, expectation value for $\hat{Q}_{2}$ and basis deformation. The list of calculations is read from the hfbtho_MASSTABLE.dat file. The second line of the file must indicate the number of calculations to be made and the list must start at the fourth line. An example file is provided with the current version. Once all the calculations have been finished a summary of the obtained binding energies is sent to the standard output. An output file with a unique label is produced for every calculation in the mass table.

The DRIP_LINES mode performs a full scan of all even-even nuclei in the nuclear chart by looking for the two-neutron and two-proton drip-lines. This mode requires USE_MPI = 1 and the number of MPI tasks must be a multiple of 11. In future versions the second requirement could be removed although significant load unbalances might occur as a result. In this mode, MPI teams are formed to run over isotopic or isotonic chains calculating the binding energy of every nucleus with 11 different deformations. The list of isotopic and isotonic chains to be explored is read from the hfbtho_STABLELINE.dat. For every nucleus the deformation with the lowest binding energy is considered the ground-state. After the ground-state is found the MPI team will move into the next even-even nuclei on the isotopic (or isotonic) chain. The drip-line is considered to be reached when the two-neutron (-proton) separation energy becomes negative. The two-neutron separation energy is defined as

[TABLE]

and an analogous expression is used for the two-proton separation energy. This process is presented schematically in figure 1.

For bookkeeping purposes and to facilitate any subsequent review of every HFB calculation, an output file with a unique label is created for each calculation. Furthermore, each MPI team creates a log file named TeamTableXXX.dat that indicates the nuclei and deformation corresponding to each label for the output files.

2.7 Other changes

**Input/output - ** We have changed the way the code records output on disk. In the past two releases of hfbtho, the name of the file where the program was recording data, in binary format, contained the number of protons and neutrons explicitly. In addition, the program was recording the matrix elements of the HFB matrix, which is dependent on the definition of the HO basis used in the run. This format was not especially adapted to large-scale calculations of, e.g., full mass tables or potential energy surfaces (PES), where it is often advantageous to use a solution at a given point to initialize the calculation of another point – the two points having possibly different basis characteristics. Therefore, we have instead adopted the same convention as in hfodd, where the program records the values of the mean field and pairing field in the quadrature mesh. In practice, the same quadrature mesh can be used in a full PES calculation without significant loss of precision. As a result of this format change, we also modified the name of the binary file to give it the generic name hfbtho_output.hel. The new format is much more flexible and enables to initialize a new calculation even when the HO basis, the nucleus, the interaction, etc., are different.

**Charge radius - ** We have adopted the following definition of the charge radius

[TABLE]

where $\langle r_{\text{p}}\rangle^{2}$ is the expectation value on the HFB vacuum of the proton radius; $\langle R^{2}_{p}\rangle=0.769$ fm2 is the proton charge radius; $\langle R^{2}_{n}\rangle=-0.1161$ fm2 is the neutron charge radius taken from [25]; $3/4M_{p}^{2}=0.033$ fm2 is the Darwin-Foldy term from [26].

**Kick-off mode - ** We have improved the handling of what was called kick-off constraints in the previous version of hfbtho. These constraints are imposed only for the first 10 iterations and are subsequently released. This was introduced originally to push the HFB solution towards the correct minimum. In version 3.00, we allow the use of regular and kick-off constraints side-by-side. For example, the user can set a constraint on the expectation value of $\hat{Q}_{10}$ and $\hat{Q}_{30}$ while using the kick-off mode for $\hat{Q}_{20}$.

3 Benchmarks and Accuracy

We give in this section several numerical benchmarks for the new features added in the current version of hfbtho. Input data files corresponding to the benchmarks for the Gogny force, the collective inertia and the fission toolkit are included in the code release. For the regularization method, we provide the input data file corresponding to the value of the cutoff of $E_{\text{pair}}=80$ MeV.

3.1 Gogny force

The implementation of the Gogny interaction has been benchmarked against the symmetry-unrestricted DFT solver hfodd. We provide two main benchmarks:

• •

An unconstrained HFB calculation in 120Sn using a deformed basis characterized by a spherical-equivalent oscillator length of $b_{0}=2.1339661000$ fm, an axial deformation $\beta=0.4$, a maximum number of shells $N_{\rm max}=20$ and a maximum number of states of $N_{\rm states}=1771$. Such a large basis allowed probing the numerical accuracy of the two-body potential matrix elements with a large number of shells.

• •

A constrained HFB calculation in 152Dy using a deformed basis characterized by a spherical-equivalent oscillator length of $b_{0}=2.12100948454$ fm, an axial deformation $\beta=0.3$, a maximum number of shells $N_{\rm max}=14$ and a maximum number of states of $N_{\rm states}=500$. We imposed only a constraint of $\langle\hat{Q}_{20}\rangle=20$ b on the expectation value of the axial quadrupole moment.

Both benchmarks were made without including the Coulomb interaction in order to test only the implementation of the D1S Gogny functional: the treatment of the Coulomb potential is slightly different in the two codes, which can lead to small differences of the order of the keV; see [27, 2]. The benchmarks are shown in Tables 1-2.

We see that enforcing strictly identical HO basis in both codes lead to agreement between kinetic energies of less than 100 eV, and of less than 500 eV on the direct part of the Gogny interaction energy. Note that we use the recently published version 2.73y of hfodd, where the implementation of the Gogny force was also benchmarked against a spherical code.

3.2 Large-amplitude collective motion

We show in table 3 the results for the components of the collective mass tensor, both in the GCM approximation of (45) and in the ATDHFB approximation of (48), as well as the zero-point energy correction. Calculations were performed in 240Pu for the SkM* parametrization of the Skyrme potential and a surface-volume pairing force characterized by $V_{0}^{(n)}=-265.2500$ MeV, $V_{0}^{(p)}=-340.0625$ MeV and a cut-off on quasiparticle energies of $E_{\mathrm{max}}=60$ MeV; see [28] for additional details. The basis was characterized by $b_{0}=2.5$ fm, $\beta=0.9$, and $N_{\text{max}}=20$. An additional cutoff on the total number of states was imposed with $N_{\text{states}}=800$. We imposed the constraints, $\langle\hat{Q}_{20}\rangle=320$ b and $\langle\hat{Q}_{30}\rangle=25$ b3/2, which correspond to a region not far from the scission point. Calculations were verified against an unpublished branch of version 2.73 of hfodd. The collective inertia module was also benchmarked separately against the triaxial Gogny code of L. Robledo and results on each component agree up to the third digit.

Next, we illustrate the calculation of fission fragments properties. The characteristics of the run were identical as for the collective inertia. Although the discrepancy between the two codes is a little larger than for previous benchmarks, especially for multipole moments, note that this is simply a consequence of performing integration over $]-\infty,z_{N}]$ or $[z_{N},+\infty[$, for which there is no exact quadrature rule. Numerical differences between the codes are also much smaller than, e.g., basis truncation errors.

3.3 Regularization of zero-range pairing force

We have tested our implementation of the pairing regularization by focusing on the nucleus 120Sn. We use a full spherical basis of $N=20$ shells (hence $N_{\text{states}}=1771$) and oscillator length $b_{0}=2.03901407814296$ fm. We solved the HFB equation for the SLy4 parametrization of the Skyrme functional. In the pairing channel, we use a surface-volume pairing force characterized by the same pairing strength for protons and neutrons, $V_{0}=-370.2$ MeV. This should give us conditions very similar to those reported in [15]; see column 4 of Table 1 and line 3 of Table 2.

We performed the calculation for several values of the cut-off energy, here denoted $E_{\text{max}}$. As in [15], we show in figure 2 the deviation between the total energy and neutron pairing gaps from their respective mean values. The mean values were computed in the range $60\leq E_{\text{max}}\leq 80$ MeV of cut-off energy, and are equal to $\langle E\rangle=-1018.396$ MeV and $\langle\Delta\rangle=1.251$ MeV. As expected, these values are a little different from calculations performed in coordinate space. In particular, the 136 keV difference in energy is compatible with systematics shown, e.g., in [29].

Overall, the variations of the energy and of the pairing gaps are about 3-4 times larger than for the hfbrad solver. This result was expected and is discussed in more details in [16]. Note that we could not reproduce the hfbtho results of [16] for three reasons: the version of the code used back then did not properly take into account the pairing rearrangement energy, resulting in small errors of typically 10-20 keV on the total energy; the treatment of the direct Coulomb energy was numerically inaccurate, as discussed extensively in [2]; finally, the value of the pairing strength used in the studies of [16] is not listed.

3.4 Multi-threading

Finally, we illustrate the advantage of multi-threading for HFB calculations. The left panel of figure 3 shows the average time per iteration of a spherical HFB calculation as a function of the size of the basis. Here, calculations are done in a full spherical basis, hence the number of states is proportional to the number of shells $N$ according to $N_{\mathrm{states}}=(N+1)(N+2)(N+3)/6$. The figure shows results obtained in 50Cr with either the SLy6 parametrization of the Skyrme potential or the D1S parametrization of the Gogny force. In the case of the Skyrme force, a surface-volume interaction was used in the pairing channel. Threaded calculations were done with 36 threads and run on the Quartz supercomputer at LLNL. As expected, OpenMP multi-threading can drastically accelerate code execution, by about an order of magnitude for the Gogny force and a factor 5 for Skyrme forces. Interestingly, multi-threaded HFB calculations with the Gogny force are only about twice slower than multi-threaded HFB calculations with the Skyrme force.

To better understand the impact of shared memory parallelism, the right panel of figure 3 shows the ratio between serial and threaded execution time for the calculation of the finite-range contribution to the mean field and pairing field only. For instance, a ratio of 10 for 10 threads indicates that the threaded code is 10 times faster, which is perfect strong scaling. For the largest basis considered here ($N_{\mathrm{max}}=16$), the code achieves an impressive 95% of maximum strong scaling. Conversely, performance saturates for small basis sizes whenever the number of threads is (approximately) twice the number of shells in the basis.

4 Input data file

The input data file format remains similar to version 2.00d and only contains a few additional namelists.

4.1 Sample input file

&HFBTHO_GENERAL number_of_shells = 10, oscillator_length = -1.0, basis_deformation = 0.0, proton_number = 24, neutron_number = 26, type_of_calculation = 1 / &HFBTHO_ITERATIONS number_iterations = 100, accuracy = 1.E-5, restart_file = -1 / &HFBTHO_INITIAL beta2_deformation = 0.0, beta4_deformation = 0.0 / &HFBTHO_FUNCTIONAL functional = ’SLY4’, add_initial_pairing = F, type_of_coulomb = 2 / &HFBTHO_PAIRING user_pairing = F, vpair_n = -300.0, vpair_p = -300.0, pairing_cutoff = 60.0, pairing_feature = 0.5 / &HFBTHO_CONSTRAINTS lambda_values = 1, 2, 3, 4, 5, 6, 7, 8, lambda_active = 0, 0, 0, 0, 0, 0, 0, 0, expectation_values = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 / &HFBTHO_BLOCKING proton_blocking = 0, 0, 0, 0, 0, neutron_blocking = 0, 0, 0, 0, 0 / &HFBTHO_PROJECTION switch_to_THO = 0, projection_is_on = 0, gauge_points = 1, delta_Z = 0, delta_N = 0 / &HFBTHO_TEMPERATURE set_temperature = F, temperature = 0.0 / &HFBTHO_FEATURES collective_inertia = F, fission_fragments = F, pairing_regularization = F, localization_functions = F / &HFBTHO_NECK use_constrain = F, neck_value = 0.5 / &HFBTHO_DEBUG number_Gauss = 40, number_Laguerre = 40, number_Legendre = 80, compatibility_HFODD = F, number_states = 500, force_parity = T, print_time = 0 /

4.2 Description of input data

We now define the new classes of input used in version 3.00.

Keyword: HFBTHO_INITIAL

Keyword: HFBTHO_INITIAL

$\bullet$ beta2_deformation = 0.0: The axial quadrupole deformation $\beta_{2}$ of the Woods-Saxon potential used to initialize the calculation. Default: 0.0

$\bullet$ beta4_deformation = 0.0: The axial hexadecapole deformation $\beta_{4}$ of the Woods-Saxon potential used to initialize the calculation. Default: 0.0

Keyword: HFBTHO_FEATURES

Keyword: HFBTHO_FEATURES

$\bullet$ collective_inertia = F: Logical switch that activates the calculation of the collective inertia tensor and zero-point energy corrections in the ATDHFB and GCM approach (both within the cranking approximation. The calculation of these quantities is only carried out at convergence. Default: F

$\bullet$ fission_fragments = F: Logical switch that activates the calculation of the characteristics of the fission fragments. Note that in order to perform such a calculation, the code must not enforce the conservation of parity (through the input force_parity). If the code detects that force_parity = T, it will override its value to False to guarantee that parity can be broken during the self-consistent loop. As a result of this condition, iterations may be noticeably slower, even though fission fragment properties are only computed at convergence. Default: F

$\bullet$ pairing_regularization = F: Logical switch that activates the regularization of the pairing strength for zero-range, density-dependent pairing forces. The value of this keyword has no effect on calculations with the finite-range Gogny force. Default: F

$\bullet$ localization_functions = F: Logical switch that activates the calculation of the localization functions at convergence. Default: F

Keyword: HFBTHO_NECK

Keyword: HFBTHO_NECK

$\bullet$ use_constrain = F: Logical switch that activates the constraint on the size of the neck. Default: F

$\bullet$ neck_value = 0.5: Value of the constraint on the expectation value of the Gaussian neck operator. Default: 0.5

5 Program hfbtho

5.1 Structure of the code

Compared with version 2.00d, the program hfbtho has been broken into several different files implementing specific Fortran modules:

• •

hfbtho_bessel.f90: modified Bessel functions of order 0 and 1;

• •

hfbtho_collective.f90: collective inertia tensor and zero-point energy correction at the perturbative cranking approximation;

• •

hfbtho_elliptic_integrals.f90: complete elliptic integral of the second kind;

• •

hfbtho_fission.f90: charge, mass and axial multipole moments of fission fragments and Gaussian neck operator;

• •

hfbtho_gauss.f90: Gauss-Hermite, -Laguerre, and -Legendre quadrature meshes;

• •

hfbtho_gogny.f90: matrix elements of the Gogny force;

• •

hfbtho_large_scale.f90: MPI parallel interface for mass table calculations;

• •

hfbtho_linear_algebra.f90: various mathematical routines;

• •

hfbtho_localization.f90: spatial localization functions;

• •

hfbtho_main.f90: main calling program;

• •

hfbtho_multipole_moments.f90: expectation value and matrix elements of axial multipole moments;

• •

hfbtho_read_functional.f90: user-defined parameters of a generalized Skyrme-like energy functional;

• •

hfbtho_solver.f90: HFB self-consistent loop, expectation value, HF and pairing field and input/output routines;

• •

hfbtho_tho.f90: transformed harmonic oscillator basis;

• •

hfbtho_unedf.f90: parameterizations of the Skyrme and Gogny functionals and density-dependent coupling constants and fields of generalized Skyrme energy functionals;

• •

hfbtho_utilities.f90: definition of integer and real types and constants;

• •

hfbtho_variables.f90: list of global variables used throughout the code;

• •

hfbtho_version.f90: version number and old history;

The programming language of most of the code is Fortran 95, while legacy code is still written, in part or totally, in Fortran 90 and Fortran 77. The code hfbtho requires an implementation of the BLAS and LAPACK libraries to function correctly. Shared memory parallelism is available via OpenMP pragmas.

The new version also comes with a built-in Doxygen documentation. To benefit from this feature, the user should install the doxygen software available at www.doxygen.org. The documentation is built by typing

make doc

By default, Doxygen generates only an on-line HTML documentation. The main page is located at ./doc/html/index.html. A PDF documentation can also be generated by going into ./doc/latex and typing

make

The PDF file is named refman.pdf.

5.2 Running the code

The program ships with a Makefile that is preset for a number of Fortran compilers. The user should choose the compiler and set the path for the BLAS and LAPACK libraries. Assuming an executable named hfbtho_main and a Linux system, execution is started by typing

./hfbtho_main < /dev/null >& hfbtho_main.out

The program will attempt to read the file named hfbtho_NAMELIST.dat in the current directory. The user is in charge of ensuring this file is present, readable, and has the proper format. The code will automatically generate a binary file of the form named hfbtho_output.hel.

As shown in section 3, HFB calculations are greatly accelerated when OpenMP multi-threading is activated. However, the user should keep in mind that this requires setting additional environment variables. In Linux/Unix machines, the default stack size is not large enough to run the code and must be increased. This can be achieved by instructions such as

ulimit -s unlimited

export OMP_STACKSIZE=32M

The value of ulimit defines the amount of stack size for the main OpenMP thread. OpenMP supports control over the stack size limit of all additional threads via the environment variable OMP_STACKSIZE. The value given above should be sufficient for all applications. Note that this value does not affect the stack size of the main thread. For completeness, note that the GNU OpenMP run-time (libgomp) recognizes the non-standard environment variable GOMP_STACKSIZE. If set it overrides the value of OMP_STACKSIZE. Finally, the Intel OpenMP run-time also recognizes the non-standard environment variable KMP_STACKSIZE. If set it overrides the value of both OMP_STACKSIZE and GOMP_STACKSIZE.

6 Acknowledgments

We are very grateful to Luis Robledo for having helped us with the benchmark of the collective inertia tensor. We also thank Karim Bennaceur for providing us with the latest version of hfbrad to test the regularization of the pairing channel, and Tomas Rodriguez for providing access to his axial code to test the Gogny force. Support for this work was partly provided through Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Nuclear Physics. It was partly performed under the auspices of the US Department of Energy by the Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 (code release number: LLNL-CODE-573953, document release number: LLNL-JRNL-587360). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Oak Ridge Leadership Computing Facility located in the Oak Ridge National Laboratory, which is supported by the Office of Science of the Department of Energy under Contract DE-AC05-00OR22725. It also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We also acknowledge “Fusion,”a 320-node cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.