DIRQFAM v2.0.0

The corresponding DIRQFAM v1.0.0 paper [1] is located in doc directory.

The corresponding DIRQFAM v2.0.0 paper [2] is located in doc directory.

The corresponding DIRHB paper [3] is located in doc directory.

Tested versions of the DIRQFAM code are located at releases page.

How to use

The DIRQFAM code is built upon the DIRHB program package [3] for the solution of the stationary relativistic Hartree-Bogoliubov equations for even-even open-shell nuclei with axially symmetric quadrupole deformation. DIRQFAM complements DIRHBZ code with the QFAM solver to calculate the multipole response for systems with axially symmetric quadrupole deformation.

Input

The input data are provided via the dirqfam.dat file, and are separated into two parts: first part is the same as in Ref. [3] and determines the input parameters for the ground state calculation, while the second part serves as an interface for QFAM parameters.

Ground state parameters (same as in Ref. [3])

• n0f, n0b:
  Number of oscillator shells used in expanding the large component of Dirac spinor (small component is expanded in n0f+1 shells) and wave functions of meson fields respectively. Both n0f and n0b must be even numbers. Recommended value of n0f depends on the nucleus and one should in principle rerun the calculation with larger n0f and compare the difference in output to establish whether the convergence is satisfying. Recommended value of n0b is at least 2(n0f+1). Since one very rarely uses more than n0f=24 shells, we recommend fixing the value of n0b=50, in which case one doesn't have to worry about the n0b parameter.

• beta0, betai:
  Deformation parameter of the oscillator basis beta0 and of the initial Woods-Saxon potentials betai respectively. In order to improve accuracy, these parameters should be close to the self-consistent ground state quadrupole deformation. For example, when dealing with a nucleus having deformation parameter β=+0.550, one should use beta0=+0.550 and betai=+0.550.

• inin:
  The starting parameters for the initial potentials and initial pairing field respectively. If set to 1, the code starts with Woods-Saxon model as initial guess for the self-consistent potentials and with diagonal pairing field respectively. Otherwise, if set to 0, the code uses the data from previous run stored in dirhb.wel and/or dirhb.del as starting potentials and pairing field respectively. For casual users, we recommend using the value of 1.

• Even-even nuclide to be computed. Element name, followed by the mass number. If the element name has only one character, it should begin with an underscore, eg. _C 12, _O 16 and _U 238. Otherwise, for two character elements simply type e.g. Zr 100.

• Init.Gap:
  Initial pairing gap (in MeV) of diagonal pairing field for protons and neutrons, relevant if inin for pairing field is set to 1. We recommend the generic value of 1 MeV.

• Force:
  Acronym of the parameter set of the selected energy density functional. In current version of the code, DD-PC1 and DD-ME2 are available.

• icstr, betac, cquad:
  The quadrupole deformation constraint control parameters. If icstr is set to 0, the quadrupole constraint is not included. If icstr is set to 1, then the constrained value of expected deformation betac is imposed with the stiffness constant cquad. We recommend the value of cquad = 0.010, but if the code fails to constrain the quadrupole moment, one should increase it keeping in mind that too large stiffness constant may disrupt the convergence of iterations. The contrained value of deformation betac is defined as: $$\beta = \sqrt{\frac{5\pi}{9}} \frac{1}{A R_0^2} \int \rho_v(\boldsymbol{r}) (2z^2-r_\perp^2) d\boldsymbol{r} = \frac{4\pi}{3} \frac{1}{A R_0^2} \int \rho_v(\boldsymbol{r}) |\boldsymbol{r}|^2 Y_{20}(\theta,\varphi) d\boldsymbol{r},$$ where $\rho_v(\boldsymbol{r})$ is axially symmetric self-consistent ground state isoscalar-vector density, and $R_0=1.2 A^{1/3}$ fm.

Mind the alignment of input parameters, for example, parameter n0f should be written in 5 character width after the equality sign. An example of dirqfam.dat file is provided and the user must follow the same alignment pattern.

QFAM parameters interface

• Calculation type:
  Value 0: Free response is calculated for a given range of energies. Value 1: Self-consistent response is calculated for a given range of energies. Value 2: Self-consistent response is calculated for a given energy and various data are printed. Value 3: Self-consistent solution is calculated along a circular contour and the contour integral is calculated.

• Include Coulomb, Include pairing:
  If set to 0/1, the Coulomb interaction or pairing is omitted/included both in ground state and in the QFAM calculation respectively.

• NGH, NGL:
  Parameters defining the size of the Gaussian quadrature grid. NGH is the number of Gauss-Hermite nodes in z>0 direction and NGL is the number of Gauss-Laguerre nodes in r direction. One should use at least NGH=max(n0f+1,n0b) and NGL=max(2(n0f+1),n0b). We recommend fixing these values to NGH=25 and NGL=50, since one rarely uses more than n0f=24 and n0b=50 shells.

• Smearing gamma:
  The imaginary part of the complex frequency (smearing width) given in MeV. Reasonable value is around 0.05-0.50 MeV. This parameter is relevant only if calculation type is se to 0, 1 or 2.

• Solver tolerance:
  Relative residual error tolerance for GMRES solver. We recommend using the value of 1.e-5, which has shown to give the strength function accurate up to 4 most significant digits.

• Arnoldi vectors:
  Maximum number of Arnoldi vectors used in GMRES solver. This is the limit on the number of QFAM solver steps. We recommend using the value of 70. If the GMRES solver fails to satisfy the relative residual error tolerance, we recommned increasing this value, however keep in mind that this means larger memory consumption of the program.

• J multipolarity, K multipolarity, Isospin:
  Values of J, K, and T, that define the multipole excitation operator. In the current version, J value is restricted to 0 <= J <= 5. Multipolarity K should be 0 <= K <= J. Isospin selects isoscalar/isovector excitation if set to 0/1 respectively.

• Omega start, Omega end, Delta omega:
  Parameters (MeV) that control the starting point, the ending point and the increment of the energy range over which the response is calculated. Relevant only if the calculation type is set to 0 or 1.

• Omega print:
  The energy for which the self-consistent solution is calculated if calculation type is set to 2.

• Omega center, Omega radius, No. points:
  Circular contour parameters used if calculation type is set to 3. The contour is a circle centered at Omega center MeV with radius Omega radius MeV. Number of integration points used for contour integration is selected via No. points parameter.

Mind the input format, all selected values should be aligned between the equality sign = and the sentinel |.

Output

The output of the calculation is divided into two parts. The first output file dirhb.out, located in the GS_output directory, contains the information on the ground state calculation. Detailed description of this file can be found in [3].

The second part of the output relevant for the QFAM calculation is located in the QFAM_output directory. The calculated strength function is written in the strength.out file. If calculation type is set to 3, the value of the contour integral is also printed in strength.out. If calculation type is set to 2, QFAM_output contains additional details and information about the self-consistent solution obtained for a given energy Omega print.

How to run

The programming language of the DIRQFAM code is Fortran and the user should provide an implementation of the BLAS and LAPACK (version 3.6.0. or higher) linear algebra libraries. Since the code depends heavily on zgemm, dgemm and dgemv subroutines, the user should provide an efficient implementation of the BLAS library. We recommend an open source implementation OpenBLAS, or freely available Intel^® oneAPI Math Kernel Library as a part of the Intel^® oneAPI Base Toolkit.

The code is compiled by standard Makefile build automation which is set to work with the GFortran compiler.

If the user invokes make command, the compilation of the code will produce the executable file run. The code is then executed by invoking the ./run command. If the user invokes make dbg, the code is compiled in debug mode in which various additional checks are performed, and the executable file dbg is generated. Since the executable produced in debug mode is considerably slower, this mode should be used only for testing and developing. If OpenBLAS is employed, the command export OPENBLAS_NUM_THREADS=4 can be invoked to select the number of threads used by OpenBLAS. If Intel^® oneAPI Math Kernel Library is employed, the command export MKL_NUM_THREADS=4 can be invoked to select the number of threads used by the Intel^® oneAPI Math Kernel Library.

Test examples

In test directory we provide a set of input files and expected output of the program. Each test contains an input and output subdirectory which contain the input and expected output data. It is advisable that the user first tries to reproduce these values.

Performance benchmark

Benchmark was done using Intel^® Core^® i7-9750H @ 2.60GHz machine (laptop). BLAS and LAPACK are provided via Intel^® oneAPI Math Kernel Library, which are forced to run using a single thread, i.e. the entire benchmark is performed using a single thread.

We select the isoscalar J=5, K=3 excitation with Gaussian quadrature grid: NGH=25, NGL=50. The DD-ME2 parametrization is used with n0b=50 shells, and 70 Arnoldi vectors stored in the memory are used by the GMRES solver.

The following table shows running time per QFAM iteration together with total memory usage. It takes roughly 30-60 iterations to reach self-consistency for a given excitation energy, depending on the self-consistency tolerance.

n0f	Memory[GB]	Time[s]
10	0.66	0.25
12	0.90	0.52
14	1.31	1.04
16	1.94	1.98
18	2.99	3.65
20	4.54	6.48

References

[1] A. Bjelčić, T.Nikšić, Comp. Phys. Comm. 253, 107184 (2020).

[2] A. Bjelčić, T. Nikšić, Comp. Phys. Comm. 287, 108689 (2023).

[3] T. Nikšić, N. Paar, D. Vretenar, P. Ring, Comp. Phys. Comm. 185, 1808 (2014).