Source-linked AI summary
Axially deformed solution of the Skyrme-Hartree-Fock-Bogolyubov equations using the transformed harmonic oscillator basis (II) HFBTHO v2.00d: a new version of the program
M. V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, S. Wild
TL;DR
Nuclear Skyrme HF/HFB calculations require a basis that handles deformation, pairing, and weak binding while supporting broader numerical treatments. HFBTHO v2.00d extends the axial transformed-harmonic-oscillator solver with new physical capabilities, generalized functionals, constraints, blocking, and parallelism, and provides benchmark comparisons with other solvers. Its numerical accuracy remains bounded by quadrature aliasing, which can make Coulomb energies depend on grid and quadrature choices.
Problem
Weakly bound paired nuclei require accurate asymptotic quasiparticle wave functions while calculations also need treatments for deformation, pairing, and broader nuclear-DFT applications.
Method
HFBTHO v2.00d expands an axial transformed-harmonic-oscillator HFB solver with generalized Skyrme-like energy densities, new HFB capabilities, validation benchmarks, and OpenMP parallelism.
Results
The release provides the listed finite-temperature, constraint, blocking, generalized-functional, multipole, symmetry-breaking, Broyden, and shared-memory capabilities, with benchmark agreement reaching the 1 eV level in specified spherical comparisons.
Takeaways & Limitations
The extended solver supports broader nuclear-structure calculations and supplies validation and performance evaluations for subsequent HFBTHO and related DFT-code development.
Takeaways & Limitations
Quadrature aliasing leaves Coulomb energies dependent, potentially significantly, on length scales and the numbers of Gauss-Legendre, Gauss-Hermite, and Gauss-Laguerre points.
Abstract
from arXiv · showhide
We describe the new version 2.00d of the code HFBTHO that solves the nuclear Skyrme Hartree-Fock (HF) or Skyrme 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 modified Broyden method for non-linear problems, (ii) optional breaking of reflection symmetry, (iii) calculation of axial multipole moments, (iv) finite temperature formalism for the HFB method, (v) linear constraint method based on the approximation of the Random Phase Approximation (RPA) matrix for multi-constraint calculations, (vi) blocking of quasi-particles in the Equal Filling Approximation (EFA), (vii) framework for generalized energy density with arbitrary density-dependences, and (viii) shared memory parallelism via OpenMP pragmas.
NEW VERSION PROGRAM SUMMARY
HFBTHO v2.00d is a FORTRAN-95, parallelized solver for self-consistent nuclear HF/HFB calculations using an axial transformed harmonic-oscillator basis. It adds finite-temperature, generalized-functional, constraint, blocking, multipole, symmetry-breaking, and OpenMP capabilities.
- Physical problem: The transformed harmonic-oscillator basis supports deformation effects and pairing correlations for weakly bound nuclei near particle drip lines.The basis is intended to reproduce appropriate asymptotic behavior of quasiparticle wave functions.
- Program scope: HFBTHO v2.00d solves self-consistent HF/HFB equations using an axial transformed harmonic-oscillator basis and generalized Skyrme-like energy densities.The solver iteratively diagonalizes the HFB Hamiltonian until self-consistency, including zero-range pairing interactions.
- New capabilities: The release implements modified Broyden iterations, reflection-symmetry breaking, axial multipoles up to λ = 8, finite temperature, multi-constraint calculations, quasi-particle blocking, generalized functionals, and OpenMP parallelism.These are the major listed revisions distinguishing version 2.00d.
- Execution: The code assumes axial and time-reversal symmetries and has been parallelized using shared-memory OpenMP pragmas.It is available for UNIX, LINUX, and Windowsxp systems and requires standard BLAS/LAPACK routines.
- Execution: Typical runtime ranges from a few seconds for small-basis ground states to a few minutes for highly deformed heavy nuclei with Nmax > 20.Runtime depends on the nucleus, basis size, requested accuracy and configuration, compiler, libraries, and hardware.
1 Introduction
The paper presents HFBTHO v2.00d as an extended solver release and establishes benchmark comparisons with other nuclear DFT codes. It also notes compatibility changes affecting users of version 1.66.
- Purpose: Version 2.00d extends the earlier transformed-harmonic-oscillator HFB method with reflection breaking, multipoles, constraints, blocking, finite temperature, and generalized functionals.The paper describes these as capabilities added to the solver.
- Purpose: The paper compares HOSPHE, HFBTHO v2.00d, and HFODD to provide validation and performance benchmarks for future code development.The comparisons cover spherical, axially deformed, and symmetry-unrestricted DFT solvers.
- Compatibility: Backward compatibility of input and output files is broken between version 1.66 and version 2.00d.The change supports future HFBTHO development and deeper integration with later HFODD releases.
- Organization: The paper reviews the new code capabilities before presenting numerical benchmarks between HFBTHO and the other DFT solvers.The benchmarks are presented as important for future program development.
2 Modifications introduced in version 2.00d
Version 2.00d expands HFBTHO’s numerical methods, physical treatments, functional flexibility, and computational implementation. The revisions include finite-temperature and constrained HFB, blocking, generalized functionals, multipole calculations, and shared-memory parallelism.
- Self-consistent iterations: The modified Broyden method updates HFB matrix elements at each iteration instead of the linear mixing used in version 1.66.The method is introduced for nonlinear self-consistent iterations.
- Axial multipole moments: Axial multipole moments are computed through lmax = 8 and can also be used as constraints in the harmonic-oscillator basis.Constraint matrix elements are evaluated numerically on Gauss-Laguerre and Gauss-Hermite integration nodes.
- Finite temperature: At finite temperature, the HFB density matrices depend on Fermi-Dirac quasiparticle occupations, and particle-number iterations include ∂fµ/∂λ.The finite-temperature equations retain the HFB form while modifying occupations and the particle-number derivative.
- Linear constraints: Multi-constraint calculations readjust Lagrange parameters using a first-order perturbation of the generalized density matrix in the cranking approximation.The procedure infers parameter changes from deviations between calculated and requested multipole moments and operates in the quasiparticle basis.
- Quasi-particle blocking: Equal-filling quasi-particle blocking enables calculations for odd-even and odd-odd nuclei while retaining time-reversal symmetry.The blocked quasiparticle is selected from overlaps between quasiparticle and single-particle wave functions.
- Generalized energy densities: The rewritten kernel supports generalized second-order Skyrme functionals with arbitrary density dependence, although more advanced functionals must be coded by the user.The implementation requires the functional and its first density derivatives, with second derivatives needed for nuclear-matter properties.
- Shared-memory parallelism: OpenMP parallelism distributes selected time-intensive operations, including HFB-matrix diagonalization across Ω-blocks and configuration-space reconstruction.The code separates the DFT solver kernel from its calling program to facilitate large-scale applications.
- Performance: With one thread, version 2.00d is slightly slower per iteration than version 1.66 because generalized-functional densities and Broyden updates add work.The paper cautions that overall performance comparisons are difficult because the versions use different calculation procedures.
3 Benchmarks and Accuracy
The paper motivates a comprehensive validation and performance comparison of HFBTHO with other density-functional solvers. The supplied passage frames these benchmarks as references for subsequent code development.
- Benchmark purpose: The benchmark suite compares HFBTHO with other DFT solvers for even-even and odd nuclei, covering both solution approaches and solver validation.The authors aim to collect these comparisons in one reference set for later developments.
3.1 Benchmarks in spherical nuclei: 208Pb and 120Sn
The spherical benchmarks compare HFBTHO with HOSPHE and HFODD for Hartree-Fock in 208Pb and HFB in 120Sn under matched basis and interaction settings. Agreement is at the eV level without Coulomb terms, while Coulomb inclusion produces a discrepancy of about 100 eV in the 208Pb test.
- Benchmark setup: HFBTHO was benchmarked against HOSPHE and HFODD for spherical Hartree-Fock calculations in 208Pb and spherical HFB calculations in 120Sn.The 208Pb calculation used SLy5, while the 120Sn calculation used UNEDF0 with density-dependent delta pairing.
- 208Pb Hartree-Fock: The 208Pb Hartree-Fock benchmark used Nmax = 16 spherical oscillator shells with oscillator length b = 2.0 fm and matched quadrature settings.The corresponding oscillator-frequency factor was f = 1.49831558 in 208Pb, with NLeg = 80 and NGH = NGL = 40.
- 208Pb Hartree-Fock: Without Coulomb potentials, HFBTHO and HFODD differ by no more than 2 eV in ESkyrme, while radii agree through at least the sixth digit.HOSPHE differs by less than 1 eV in energy, with a proton-radius difference of 0.0002 fm.
- 208Pb Hartree-Fock: Including direct and exchange Coulomb potentials increases the 208Pb total-energy discrepancy to approximately 100 eV.The kinetic-energy contribution is identified as particularly sensitive to numerical implementation details.
- 120Sn HFB: The 120Sn HFB benchmark used UNEDF0 with mixed volume-surface density-dependent delta pairing and a sharp quasiparticle cutoff of Ecut = 60 MeV.The pairing strengths were fitted with the UNEDF parameterization and therefore were not adjusted independently.
3.2 Benchmarks in even-even deformed nuclei: 240Pu
The 240Pu benchmark tests HFBTHO and HFODD with a strongly deformed, relatively small basis where basis-selection differences are amplified. Although the prescriptions are analytically equivalent for a prolate ellipsoid, the resulting energy difference can reach 300 keV without Coulomb terms and returns to the 1 eV level when oscillator lengths match.
- Benchmark purpose: The deformed benchmark compares HFBTHO and HFODD basis-truncation schemes for a strongly deformed even-even 240Pu configuration.Accurate calculations require a suitably deformed or stretched harmonic-oscillator basis characterized by directional oscillator frequencies and a finite number of retained states.
- HFBTHO stretched basis: The HFBTHO stretched basis is defined from an axially deformed prolate ellipsoid, with its harmonic-oscillator lengths uniquely fixed by spherical length b0 and basis deformation β.The construction uses deformation parameters and volume-conservation relations to determine the basis scales.
- Basis comparison: For a prolate ellipsoid with β > 0 and γ = 0°, the HFODD and HFBTHO prescriptions are analytically identical, but numerical radius determination produces slightly different oscillator frequencies.Those frequency differences alter which basis states are selected.
- Numerical benchmark: 300 keV is the maximum energy difference between the codes without Coulomb terms in the extreme 240Pu test.The discrepancy is attributed to slightly different oscillator lengths, magnified by large deformation and the relatively small Nmax = 16, Nstates = 500 basis; Coulomb terms do not qualitatively change the picture.
- Numerical benchmark: 1 eV-level agreement is recovered without Coulomb when oscillator lengths are manually made identical, or when a spherical basis is used.At β = 0.7, the reported transverse and axial oscillator lengths differ between HFODD and HFBTHO.
3.3 Benchmark in deformed odd nuclei: 159Ba
The deformed odd-nucleus benchmark evaluates HFBTHO equal-filling blocking against HFODD for 159Ba using a stretched harmonic-oscillator basis. Agreement is about 1 keV, with maximum deviations of 3.6 keV, because the less heavy and less deformed case reduces basis-truncation effects.
- Blocking framework: HFBTHO supports blocking calculations in odd-even and odd-odd nuclei through the equal filling approximation because time-reversal symmetry is conserved.Under this prescription, the time-odd Skyrme fields are identically zero.
- Benchmark setup: The 159Ba benchmark used 158Ba as the even-even core, SLy4, mixed surface-volume pairing with V0 = −300 MeV, Ecut = 60 MeV, β = 0.2, and NStates = 500.The calculation used Nmax = 16 shells, oscillator length b = 2.2 fm, and NGL = NGH = 40.
- Numerical agreement: Approximately 1 keV agreement was obtained between HFBTHO and HFODD, with maximum deviations of 3.6 keV.Both codes used the same basis-selection rules in their respective implementations, while Coulomb terms were neglected in the related numerical tests.
- Interpretation: The smaller mass and weaker deformation of 159Ba mitigate basis-truncation effects, so oscillator-scale differences have less impact than in the extreme 240Pu benchmark.The reported agreement is therefore consistent with the preceding benchmark results.
3.4 Transformed harmonic oscillator basis: 90Ni
HFBTHO generates a transformed harmonic oscillator basis through a local scale transformation designed to enforce proper density asymptotics. In 90Ni, its THO results agree very well with an independently developed spherical implementation, with only tiny deviations beyond r = 8 fm.
- Basis construction: The THO basis is generated by applying a local scale transformation to harmonic-oscillator single-particle functions.The transformation acts on each point (ρ, z).
- Basis construction: The local scale transformation is chosen to enforce exponential decay of the density at large distances.
- Benchmark setup: The benchmark compares HFBTHO with an independently developed spherical implementation using the same transformed-basis construction.HFBTHO generated the scale-transformation function and its partial derivatives for the comparison.
- Benchmark setup: The 90Ni benchmark uses the SLy4 functional, surface pairing, 20 spherical shells, b = 2.0 fm, and 40 Gauss-Laguerre and Gauss-Hermite points.The pairing cutoff is Ecut = 60 MeV, and Coulomb terms are neglected.
- Benchmark result: Agreement between the two implementations is very good, with tiny neutron-density deviations beyond r = 8 fm.The deviations arise from quantizing the transformation in HFBTHO and using that numerical function in the spherical code.
3.5 Benchmark at finite temperature: 50Cr
HFBTHO v2.00d implements finite-temperature HFB and is benchmarked against HFODD for 50Cr. The comparison uses a spherical basis and a calculation at T = 1.5 MeV, with an entropy correction noted for HFODD v2.49t.
- Finite-temperature implementation: HFBTHO v2.00d implements the finite-temperature HFB equations and compares them with HFODD in 50Cr.
- Benchmark setup: The benchmark uses a full spherical basis of Nmax = 12 shells, the SLy4 interaction, surface-volume pairing, and Ecut = 60.0 MeV.The oscillator length is b0 = 1.7622146 fm, equivalent to f = 1.2.
- Benchmark setup: The temperature was set at T = 1.5 MeV for the finite-temperature comparison.
- Benchmark caveat: HFODD version 2.49t contains an entropy bug: its reported entropy should be multiplied by a factor 2.
3.6 Precision of the Coulomb term
HFBTHO’s direct Coulomb calculation is affected by aliasing from finite quadrature, making the energy depend on quadrature settings rather than converging straightforwardly with Legendre points alone. Increasing spatial quadrature precision mitigates this dependence, while recommended defaults balance accuracy and speed.
- Origin of the numerical error: Finite quadrature of the Coulomb kernel and spatial coordinates introduces aliasing, so the calculated Coulomb energy depends on L and the quadrature-point counts.The quantized 1/|r − r′| representation has a singularity-related systematic error that remains nonzero for practical grids.
- Legendre integration: The Coulomb energy shows no asymptotic convergence with increasing NLeg alone; instead, it develops an L-dependent plateau.For 208Pb with Nmax = 16 and SLy5, the dependence on L is clearly visible when varying Gauss-Legendre points.
- Spatial quadrature: Increasing NGH and NGL at fixed L = 5 fm clearly reduces the Coulomb energy’s dependence on Legendre quadrature.This indicates that the dominant practical error is associated with finite spatial integration rather than the Legendre integration itself.
- Integrand diagnostic: Coarse spatial meshes make the integrand bend upward as ξ → 1, although it should tend to 0; finer meshes mitigate this nonphysical behavior.The upward bend is identified as a manifestation of aliasing in the close-up of the integrand.
- Recommended settings: HFBTHO defaults to L = 50 fm and NLeg = 80; NGH = NGL = 40 is sufficient for ground-state properties, while very deformed configurations require finer Gauss integration.The authors also identify a future Fock-space Coulomb treatment using Moshinsky transformations as a way to eliminate aliasing errors.
4 Input data file
Version 2.00d replaces the old input format with a flexible Fortran namelist and exposes settings for basis construction, functionals, pairing, constraints, blocking, and convergence. The documented keywords specify both calculation choices and their defaults or allowed behaviors.
- File format: The input format changed entirely from version 1.66 to a more flexible Fortran namelist structure.The change accommodates the larger number of features introduced in version 2.00d.
- Core calculation settings: Core namelist entries select shell count, oscillator length, basis deformation, proton and neutron numbers, calculation type, iterations, accuracy, restart mode, and functional.The calculation type distinguishes standard HFB from HFB+LN, while convergence uses the norm difference between successive HFB matrices.
- Constraints: Multipole constraints activate selected moments Q_l for l = 1,…,8 and assign their target values in internally scaled units.A negative activation value triggers a kickoff mode that releases the constraints after up to 10 iterations.
- Blocking and deformation: Blocking inputs encode a Nilsson state and determine whether the resulting proton calculation corresponds to the (Z + 1, N) or (Z − 1, N) nucleus.The sign of 2Ω selects the neighboring proton number, while the basis supports axial prolate or oblate deformation through β2.
5 Program HFBTHO v2.00d
The HFBTHO v2.00d program is organized as a self-contained solver plus a calling program, with Fortran modules covering transformations, energy-density calculations, integration, and shared numerical utilities. It is compiled with BLAS/LAPACK support and run from a namelist-driven executable workflow.
- Program structure: HFBTHO v2.00d is distributed as a main solver file and a separate calling-program file.The solver file contains the principal Fortran modules and the calling program launches calculations.
- Build requirements: The code is primarily Fortran 95, retains some Fortran 90 and Fortran 77 code, requires BLAS and LAPACK, and provides shared-memory parallelism.Compilation uses a supplied Makefile after selecting the compiler and library paths.
- Modules: Modules implement the UNEDF energy density and fields, THO transformation, integration meshes, interpolation, numerical types, Bessel functions, and elliptic integrals.The HFBTHO module stores public variables used throughout the code, while some routines remain outside modules.
- Execution: The executable is started on Linux with ./main < /dev/null >& main.out and reads NAMELIST.dat from the current directory.The user must ensure that the namelist file exists and is readable.
- Restart workflow: Restart files encode spherical, prolate, or oblate shapes and distinguish regular HO runs from THO runs through filename extensions.The shape markers are s, p, and o, while the extensions are hel and tel.