HFODD (v2.40h) User's Guide

We describe the input data and installation procedures of the code HFODD (v2.40h). The present write-up contains complete and comprehensive information that has originally been given in six independent publications. It is enhanced by the subject index and indexes of variables, input-data keywords, subroutines, and files that are used in this user guide.

, IV [4], V [5], and VI [6]. The present user's guide summarizes articles I-VI and provides practical details presented in the form that is suitable for a rapid reference on how to run the code. It also contains details updated to those corresponding to version (v2.40h). Information pertaining to the physics issues and numerical methods are not repeated here, and the user is referred to the original publications and to the references given therein.
Together with consulting the present user's guide, the reader is invited to visit the hfodd home page (http://www.fuw.edu.pl/∼dobaczew/hfodd/hfodd.html), where many useful files are available for downloading and which contains information on future versions and extensions being developed.
The present user's guide overrides the previous version published in Ref. [7], which is now obsolete.

INPUT DATA FILE
The code reads the input data from the standard FORTRAN input file, and a general structure of the input data file is defined by the following rules: 1. Input data file is an ASCII file composed of independent input items.
2. An item is composed of the keyword line, and of the data line which follows immediately the keyword line. Two items contain no data lines, see point 6 below.
3. The keyword line contains in columns 1 through 10 a keyword which is a specific text defining the item. If the keyword has less than 10 characters it has to be obligatorily padded with appropriate number of trailing spaces. Columns beyond 10 are ignored and can be used to place arbitrary comments or texts.
4. The data line contains a prescribed number of data values which are read in the free FORTRAN FORMAT. The type of data (REAL, INTEGER, or CHARACTER) should match the definition of the given data line. The CHARACTER data should be provided starting from the 13-th column of the data line.
5. Items can be separated by an arbitrary numbers of comment lines. A line is treated as a comment and ignored provided it does not contain in columns 1 through 10 any of the valid keywords.
6. Two items contain no data lines. The first one is defined by the keyword EXECUTE, and starts a calculation for the currently defined input parameters. The second one is defined by the keyword ALL DONE, and terminates the program. The latter item is not required if the compiler is able to properly recognize the end of input data file. When the end of input data file is encountered, the program terminates as if the keyword ALL DONE was found.
7. A given item may occur any number of times in different places of the input data file. Only the last one of the same items occurring before a given EXECUTE keyword is taken into account.
8. The items appearing between two consecutive items EXECUTE define the data set, i.e, the set of input parameters, for which the calculation will start at the moment when the latter of the two such items is found. Between two consecutive items EXECUTE, the order of other items is arbitrary.
9. All input parameters have the default values defined in the code (subroutine predef). Therefore, if the input data file contains only the item EXECUTE, and no other lines, a calculation for the default values is performed. In the present version, this results in performing 50 iterations for the superdeformed state of 152 Dy at the angular frequency of hω=0.5 MeV.
10. Only those items which define the values of input parameters which are different from the default values have to be included in the input data file. The values of input parameters defined by any item stay in effect till another occurrence of the same item, or till the code terminates. This rule is valid irrespective of how many EXECUTE items follow the given item.
In the following subsections we list all the valid keywords and their corresponding data lines. Titles of the subsections are introduced only to organize the keywords into groups pertaining to specific categories, but are otherwise irrelevant. The default values are given in data lines following the keyword lines. In the data lines we also show names of the corresponding variables, which are referred to in the text.

General data
This section lists keywords pertaining to definitions of particle numbers, iterations, and phase space. In the latter case, the odd and odd-odd nuclei are described within the false-vacuum approximation unless the quasiparticle blocking approximation is used for IPAHFB=1, see Secs. 2.7-2.7. For IPABCS=1, the blocking approximation is not yet implemented. For IN FIX=0 or IZ FIX=0, calculations are performed for one type of particles only, which allows to treat the atomic gases.

= NOITER
For NOITER>0, the specified number of iterations is performed. Unless ICONTI=1 or IREAWS=1, the code starts with iteration number 0 by constructing the initial mean fields based on the Nilsson potential. Specifying NOITER=0 requests only this initial phase. If the iteration is restarted from the previously recorded potentials (ICONTI=1), then the counting of iterations continues from the number of iterations performed in the previous run. Along with requesting the initial D T 2h transformations (INIINV>0 or INIKAR>0) or initial rotation (INIROT=1) one must set NOITER=1.

= NULAST
These two parameters govern the termination of the HF iterations according to the achieved stability of solutions. The stability of the HF energy has been defined in (I-37) as the difference between the total energies calculated from the single-particle energies and from the Skyrme functional. The HF iterations continue until the absolute value of the stability is smaller than EPSITE (in MeV) over NULAST consecutive iterations. When this condition is fulfilled, iteration procedure terminates and the final results are printed. This allows for an automated adjustment of the number of iterations that are required to achieve a given level of convergence. The number of iterations NOITER, can now be set to a large value, at which the iterations terminate if a stable solution is not found.
The default value of EPSITE=0.0 deactivates this option, and then the code continues up to NOITER iterations regardless of the achieved stability. If a non-zero value of EPSITE is used, a non-zero value of NULAST has to be used too. In practice, a value of NULAST=5 prevents an accidental termination of iterations in all cases when the stability energy gradually goes through zero, but the solution is not yet self-consistent.
The code is able to detect the "ping-pong" divergence described in Sec. III-2.6, i.e., the situation when the HF iteration procedure gives oscillating results in every second iteration. Upon continuing the iteration, both sequences of results, i.e., those which correspond to the iteration numbers being even and odd, stay different but perfectly stable, and hence the correct selfconsistent solution is never attained.
The code recognizes such a situation by calculating the averages and variances of the stability energy (I-37), separately in the even and in the odd sequences of results, over the last NUPING pairs of iterations. The "ping-pong" divergence condition occurs when both variances become a factor EPSPNG smaller than the absolute value of the difference of the corresponding averages, i.e., when where and Here, n denotes the number of the last accomplished HF iteration, and δE n ′ denotes the stability energy (I-37) obtained in the n ′ -th iteration. The default value of EPSPNG=0.0 deactivates this option, and then the code continues up to NOITER iterations regardless of the "ping-pong" divergence. If a non-zero value of EPSPNG is used, a value of NUPING>1 has to be used too. In practice, values of EPSPNG=0.01 and NUPING=3 allow for an efficient detection of the "ping-pong" divergence condition.
Upon discovering the "ping-pong" divergence, the HF iterations are terminated and a table of absolute values of maximum differences of single-particle observables between the two sequences of iterations is printed, see Sec. III-2.6. These maximum differences are determined for states in each of the charge-parity-signature, charge-simplex, parity, or charge blocks, depending on the conserved symmetries, and separately for particle and hole states. Whenever such a maximum difference is found for a particle state and for a hole state with adjacent indices, such a pair is proposed as a candidate for the diabatic blocking calculation.

= NUCHAO
The code is able to detect the "chaotic" divergence which occurs when the HF iterations give results which chaotically vary from one iteration to another one. The code recognizes such a divergence by finding the local maxima M k , k=1,2,. . . , in the sequence of absolute values of the stability energies (I-37), obtained in the entire series of the HF iterations performed. The "chaotic" divergence condition occurs when the code finds NUCHAO positive differences M k −M k−1 . When this condition occurs, iteration procedure terminates and the final results are printed.
For NUCHAO=0 (the default value) the code does not check whether the "chaotic" divergence occurs or not. In practice, a value of NUCHAO=5 allows for an efficient detection of the "chaotic" divergence condition. However, for a small value of NUCHAO and a small value of EPSPNG, the "ping-pong" divergence can sometimes be mistaken for the "chaotic" divergence. If one is interested in the diabatic-blocking data, printed after the "ping-pong" divergence, the recommended value of NUCHAO=5 should be increased to 10 or more. Numbers of the lowest mean-field eigenstates which are kept after the diagonalization of the mean-field Hamiltonians in the four charge-simplex blocks: (s, q)=(+i, n), (−i, n), (+i, p), (−i, p). All other eigenstates are discarded. If any of these numbers is equal to zero (the default value), the code sets it equal to the number of neutrons IN FIX (for q=n) or protons IZ FIX (for q=p). When no symmetry is conserved, the numbers of states kept in charge blocks are determined by sums NUMBSP(0,0) + NUMBSP(1,0) and NUMBSP(0,1) + NUMBSP(1,1) For calculations without pairing, i.e. for IPAIRI=0, the user is responsible for using values of NUMBSP large enough to accommodate all wave functions which might be useful for the required vacuum and particle-hole configurations, see Sec, II-3.4. In practice, the use of the default values described above is recommended as a safe option. The size of the matrices defined by the NDSTAT parameter can be reduced by the user when small values of NUMBSP are used.
For calculations with the BCS pairing, i.e. for IPAIRI=1 and IPAHFB=0, values of NUMBSP define the pairing window in which the BCS equations are solved.
For calculations with the HFB pairing, i.e. for IPAIRI=1 and IPAHFB=1, values of NUMBSP define the numbers of canonical states kept after diagonalizing density matrices. They are irrelevant for the pairing window, which is then defined by the cutoff energy ECUTOF.

Interaction
This section lists keywords pertaining to definitions of coupling constants and interaction parameters.  Parameters of several standard Skyrme forces are encoded within the program. Calculation for a given standard force can be requested by specifying in the input data file its acronym SKYRME. Along with the force parameters, for each force there is encoded information on how the given force should be used, i.e., with which value of the parameterh 2 /2m, and with which treatment of tensor, spin-orbit, and center-of-mass terms. For ISTAND=1, calculations are performed with these features set in the way specific for the given force, and the rest of the switches read on the same line is ignored. For ISTAND=0, the switches define nonstandard features in the following way: • For KETA J=1 the code takes into account in the functional, and for KETA J=0 neglects, the so-called tensor J 2 terms. The second option is equivalent to setting in Eq. (I-12) coupling constants C J t =0 and C T t =0.
• For KETA W=0 or 1 the code uses the traditional (W ′ 0 ≡ W 0 ) or generalized (W ′ 0 = W 0 ) spin-orbit term, respectively, where For forces that use generalized spin-orbit term, option KETA W=0 sets W ′ 0 equal to W 0 , while for forces that use traditional spin-orbit, option KETA W=1 has no effect. For KETA W=2, strengths W 0 and W ′ 0 are set equal to the values of W0 INP and W0PINP, respectively, which are read under the keyword of SPIN ORBIT. Note that if KETA W=2 is used, and keyword SPIN ORBIT is not specified in the input data file, the default values of W0 INP and W0PINP supersede those that are encoded within the program for the given Skyrme force.
• For KETACM=0 the code uses the traditional one-body center-of-mass correction before variation, i is the one-body kinetic-energy operator. For KETACM=1 the code uses the two-body center-of-mass correction after variation, whereP =−ih A i=1 ∇ i is the total linear momentum operator. Value of KETACM=2 is reserved for a future implementation of the two-body center-of-mass correction before variation. For KETACM=3 the center-of-mass correction is neglected.
• For KETA M=0 the code uses the value ofh 2 /2m=20.73620941 MeV fm 2 , which was encoded in version (v1.75r). For KETA M=1 the code uses the value specific for the given Skyrme force. For KETA M=2 the code uses the value specified in the input data file under keyword HBAR2OVR2M.

= HBMINP
Value of theh 2 /2m parameter, which is used if KETA M=2 is set under keyword SKYRME-STD, and ignored otherwise. Strengths W 0 and W ′ 0 of the generalized spin-orbit interaction (4), which are used if KETA W=2 is set under keyword SKYRME-STD, and ignored otherwise. Three-digit steering switch LANODD is used to specify the type of operations performed in order to define selected time-odd coupling constants based on values of the Landau parameters. The switch is followed by values of x 0 , x 1 , g 0 , g ′ 0 , g 1 , and g ′ 1 , as given by equations: where β = (3π 2 /2) 2/3 and For LANODD=000, coupling constants are not defined by the Landau parameters, and the rest of parameters read on the same data line is ignored. For LANODD=111, Eqs. (7) are solved, and coupling constants C s For other values of LANODD, the following three steps are performed in sequence:

1.
Step one: If the rightmost digit of LANODD is equal to 1 then the coupling constants C T 0 and C T 1 are determined form g 1 and g ′ 1 , Eqs. (7e) and (7f); otherwise these two coupling constants are determined from the Skyrme force parameters.

2.
Step two: If the middle digit of LANODD is equal to 1 then the coupling constants C s 0 [ρ sat ] and C s 1 [ρ sat ] are determined form g 0 and g ′ 0 , Eqs. (7c) and (7d); otherwise these two coupling constants are determined from the Skyrme force parameters.

3.
Step three: If the leftmost digit of LANODD is equal to 1 then the coupling constants C s 0 [0] and C s 1 [0] are determined form x 0 and x 1 , Eqs. (7a) and (7b); otherwise these two coupling constants are determined from the Skyrme force parameters.
All combinations of zeros and ones are allowed in LANODD. Values of parametersh 2 /2m, ρ sat , and m * /m, respectively, which are used when solving Eq. (7). If a negative number is read for any of these parameters, then the program uses the corresponding value calculated from parameters of the given Skyrme force.

Keyword: EVE SCA TS
1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. SRHO T, SRHO S,SRHODT, SRHODS,SLPR T, SLPR S, STAU T, STAU S, SSCU T, SSCU S, SDIV T, SDIV S By using this item, the coupling constants corresponding to a given Skyrme parameter set can be arbitrarily scaled. This allows calculations with modified Skyrme functionals. The time-even coupling constants in the total-sum representation, are multiplied by the 12 numbers from SRHO T to SDIV S. Variables with names ending with T and S multiply the "total" and "sum" coupling constants, respectively. Variables with names containing the acronyms RHO, LPR, TAU, SCU, and DIV, multiply the coupling constants with superscripts ρ, ∆ρ, τ , J, and ∇J, respectively, and those with RHOD multiply the densitydependent part of C ρ . Similar name convention is used for many other variables in the code hfodd.
which is used in the HFB pairing or BCS pairing calculations for neutrons. In case the values of V 1 and α allow it, the code calculates the value of the saturation density ρ 0 that gives the equivalent form of the pairing form-factor (11). In case ρ 0 cannot be calculated, the codes set its value to 1; ρ 0 is calculated only for the purpose of information, while internally the code uses only the value of V 1 . For IDEFIN=1, or for IPAHFB =1 and IPABCS =3, parameters PRHO N, PRHODN, and POWERN are ignored.  Variables YUKAGT, YUKAG0, YUKAG1, and YUKAG2 correspond, respectively, to the coupling constants g, g 0 ,ḡ 1 , andḡ 2 , see Eq. (VI-59). For I YUKA=2 or 3, the direct matrix elements of the Yukawa interaction (VI-59) are, in addition, added to the self-consistent mean field. For I YUKA=2 or 4, the exchange matrix elements are added. For IYUTYP=1, expression (VI-61) is used, while for IYUTYP=2, an analogous six-Gaussian expression is used without the short-range correction (VI-58), that is, for f (r) = 1. For I YUKA=0, all these input data are ignored and the Yukawa interaction is not taken into account. Figures 1 and 2 show the actual CPU times required in calculations using the Skyrme and Yukawa interactions, respectively. At N 0 = 10, the latter take more than two orders of magnitude longer and, moreover, their CPU times scale like N 7 0 rather than N 4 0 . A local EDF clearly makes for easier computing. shows that these times scale as N 4 0 .

Keyword: COULOMBPAR
7, 1, 1 = ICOTYP, ICOUDI, ICOUEX For ICOUDI=0, 1, or 2, the Coulomb direct energy and Coulomb direct mean field are neglected, calculated by using the Green-function method, see Section I-5, or calculated by using the Gaussian-expansion method, see Section VI-2.10, respectively. Similarly, for ICOUEX=0, 1, or 2, the Coulomb exchange energy and Coulomb exchange mean field are neglected, calculated by using the Slater approximation (I-19) or calculated exactly by using the Gaussian-expansion method, see Section VI-2.10, respectively. For the Gaussian-expansion method, that is, for ICOUDI=2 or ICOUEX=2, positive values of ICOTYP denote the numbers of Gauss-Legendre nodes used in the integral of Eq. (VI-103). For ICOUDI=2 or ICOUEX=2, the iteration can later be smoothly continued (IFCONT=1, see the keyword CONTFIELDS) only by saving the matrix elements of the mean field, that is, by requesting IWRIFI=1, see the keyword FIELD SAVE. Therefore, ICOUDI=2 or ICOUEX=2 and ICONTI=1 requires IFCONT=1. In Fig. 3,the error in the exact Coulomb exchange energy is plotted as function of the number of Gauss-Legendre nodes N C . A quite precise estimate is obtained for N C = 7 Gauss-Legendre nodes and the machine precision is obtained by doubling this number (along with doubling the CPU time). For N C = 7, these CPU times are shown in Fig. 4, which illustrates the fact that these times scale with the number of HO shells as N 7 0 .

Keyword: INSERT HO
, an external HO potential is added to the self-consistent mean field. Parameters of the potential are identical to those defining the HO basis.

Symmetries
This section lists keywords pertaining to selection of conserved, broken, and restored symmetries.

= ISIGNY
For ISIGNY=1, calculation with signature conserved are performed, while for ISIGNY=0 the signature is broken. Value of ISIGNY must be consistent with switches ISIMPY, IPARTY, ISIMTX, and ISIMTZ, see Tables 1 and 3.

= IROTAT
For IROTAT=1, calculation with time-reversal breaking are performed, while for IROTAT=0 the time-reversal symmetry is conserved. In the latter case the calculations are performed only for one value of the simplex, s=+i, which gives almost twice shorter CPU times. IROTAT=0 is incompatible with providing a non-zero value of the angular frequency or with attempting calculation for odd or odd-odd nuclei without pairing, IPAIRI=0. IROTAT=1 is incompatible with IPABCS>0, i.e., the rotating solutions can be obtained for the HFB pairing but not for the BCS pairing.   For ISIMTY=1, calculation with y-simplex T conserved are performed, while for ISIMTY=0 this symmetry is broken. Value of ISIMTY must be consistent with switches ISIMPY and IROTAT, see Table 2. For ISIMTY=−1 (the compatibility mode), the code sets ISIMTY to 1-IROTAT and requires that ISIMPY=1.

1, 1 = ISIMTX, ISIMTZ
For ISIMTX=1 and/or ISIMTZ=1, calculation with conserved symmetries given by x-simplex T and/or z-simplex T are performed, respectively. These symmetries are broken for ISIMTX=0 and ISIMTZ=0. Values of ISIMTX and ISIMTZ must be consistent with that of ISIGNY, see Table 3.

Keyword: TSIMPLEX3D
1, −1, 1 = ISIMTX, ISIMTY, ISIMTZ This keyword allows to simultaneously input all the three switches that define the three Tsimplexes. It is equivalent to using keywords TSIMPLEXES and TSIMPLEX Y together.

= IPAIRI
For IPAIRI=0, calculation without pairing correlations (pure HF) are performed, while for IPAIRI=1 the pairing correlations are included. In the latter case, for IPAHFB=0 or 1 pairing is included within the BCS or HFB method, respectively. IPAIRI=1 requires either IPAHFB=1 or IPABCS>0 and is incompatible with simultaneously setting IPAHFB=1 and IPABCS>0, i.e., the HFB pairing and BCS pairing cannot be used simultaneously. IPAIRI=1 is incompatible with IROTAT=1, unless IPAHFB=1, i.e., the rotating solutions can be obtained for the HFB pairing but not for the BCS pairing. Value of IPABCS=−1 is allowed for the sake of compatibility with earlier versions of the code, before (v2.13f). Then, IPABCS=0 is set for IPAHFB=1 and IPABCS=IPAIRI is set for IPAHFB=0.
whereP is the space inversion,T is the time-reversal symmetry,R k is the signature (rotation by angle π about the k = x, y, or z axis), andP T =PT , . For INIKAR=INIINV=0, no transformation is performed and this option is inactive.
Transformations are performed at the level of the densities, after the first iteration. As a security measure, nonzero values of INIKAR and INIINV require NOITER=1. Such values also require SLOWEV=SLOWOD=SLOWPA=0.0, so as not to mix the old and new potentials corresponding to the original and transformed densities, respectively. They are also incompatible with IPRGCM =0.
A given initial D T 2h transformation must be accompanied by the correspondingly broken symmetries, that is, INIINV=1 or 3 requires IROTAT=1 and INIINV=2 or 3 requires IPARTY=0. A rotation about the z axis must be accompanied by the broken signature, that is, ALPINI =0 or GAMINI =0 requires ISIGNY=0. For IPRGCM=1 and 2, the code calculates the diagonal and non-diagonal GCM kernels, respectively. In addition, for NUAKNO =1 or NUBKNO =1, the AMP kernels are calculated, as described in Section VI-2.1, although in the present version (v2.40h) this option is not yet available for non-diagonal kernels of IPRGCM=2. The angular momentum projection (AMP) is performed for doubled angular momenta from IPROMI to IPROMA and requires a spherical HO basis of hω x =hω y =hω z . For even (odd) particle number IN FIX+IZ FIX, IPROMI, IPROMA, and KPROJE must be even (odd). For IPRGCM=1 and 2, the present version (v2.40h) also requires IROTAT=1 and IPAIRI=0.
The AMP has been tested up to the values of angular momenta of 70h, and therefore, IPROMA must not be larger than 2*70=140. After further tests, higher values could be used by increasing the parameter JMAX=70 in function dsmalg, which calculates the Wigner functions, see Eqs. (VI-1), (VI-4), and (VI-5).
NUAKNO is the number of Gauss-Tchebyschev nodes, which are used to perform integrations over the α and γ Euler angles. For NUAKNO=1, these integrations are not performed (1D AMP) and the states are assumed to be axial with the doubled projection of the angular momentum on the z axis equal to KPROJE. NUBKNO is the number of Gauss-Legendre nodes, which are used to perform integrations over the β Euler angle. For NUAKNO>1 and NUBKNO>1, a full 3D AMP is performed and the value of KPROJE is ignored. NUAKNO>1 requires ISIMPY=0 and ISIGNY=0. IPROMA must be larger than the absolute value of KPROJE.
For IPRGCM=2, the code calculates the GCM kernels between the states labeled by three-digit indices from "000" to "999". Indices of the "left" states vary between IFRWAV and ITOWAV, and these states are read from the disc. The index of the "right" state equals to ITOWAV, and this state is equal to the current state. In addition, for IWRWAV=1, the current state is saved on the disc with the index of ITOWAV. This allows for a simultaneous calculation of the "right" state along with calculating its kernels with all previously calculated and stored "left" states. The states can also be stored on disc without calculating kernels in the given run, that is, by setting the value of ITOWAV along with IWRWAV=1 and IPRGCM=0. Names of files on the disc are composed by concatenating the three-digit index, "-", and FILWAV.

Keyword: SAVEKERNEL
For ISAKER=1, the code attempts reading the kernel files Nxxx-Lyyy-Rzzz-//FILKER, where // denotes concatenated strings. The three-digit indices are: • xxx is the consecutive index of the kernel file, • yyy is the index of the left wave function, • zzz is the index of the right wave function.
In the work directory, the file names for all indices xxx are scanned, starting from 001. The kernels stored in these files are read into memory and are not recalculated. The kernels that have not been found in the kernel files are calculated and stored in the kernel file with the lowest available index xxx. In this way, one can submit many parallel jobs, see the keyword PARAKERNEL, that calculate kernels for different values of the Euler angles α, β, and γ. The results are then collected in different kernel files with indices xxx attributed automatically. If any of the jobs is terminated before completing its task, the same input data can be resubmitted and the calculation automatically continues from the point where it has been interrupted. Once all the kernels will have been calculated, which requires a large CPU time, the AMP can be performed within a very small CPU time by reading, again automatically, all the created kernel files. ISAKER=1 requires IPRGCM>0 2.3.18. Keyword: CHECKERNEL 1 = ICHKER The names of kernel files are saved within these files. As a security measure, when reading the kernel files, their names are cross-checked against the saved information. This cross-checking can be deactivated by using ICHKER=0. This option is useful whenever the kernel files have been renamed for any reason.
2.3.20. Keyword: TRANSITION 2, 1, 0 = NMURED, NMARED, NSIRED Maximum numbers of transition electric, transition magnetic, and transition surface or Schiff moments, respectively, for which proton kernels and reduced matrix elements are calculated, printed, and stored in the kernel files. NMARED and NSIRED must not exceed NMURED. For NMURED=0, NMARED=0, or NSIRED, the corresponding proton kernels and reduced matrix elements are not calculated.

Configurations with no conserved symmetries
This section lists keywords pertaining to definitions of configurations in the case when no symmetries are conserved, and the single-particle or single-quasiparticle space is divided into two charge blocks (neutrons and protons).

Keyword: PHNONE NEU
1, 00, 00 = NUPAHO, KPNONE, KHNONE Neutron particle-hole excitations pertaining to the situation when all neutrons are in one common block (no simplex, signature, or parity is conserved). NUPAHO is the consecutive number from 1 to 5 (up to five sets of excitations can be specified in separate items). Particles are removed from the KHNONE-th state and put in the KPNONE-th state. At every stage of constructing excitations the Pauli exclusion principle has to be respected (particle removed from an occupied state and put in an empty state). Values equal zero have no effect. Note that for all neutrons sitting in one common block the reference configuration from which the particle-hole excitations are counted is defined by the total number of neutrons. These particle-hole excitations are ignored unless ISIMPY=0, IPARTY=0, and IPAIRI=0.

Keyword: PHNONE PRO
1, 00, 00 = NUPAHO, KPNONE, KHNONE Same as above but for the proton particle-hole excitations. The diabatic blocking of neutron single-particle configurations in the situation when all neutrons are in one common block. The rules to define the blocking are analogous to those described for the conserved parity, keyword DIAPAR NEU. For |IDSIZN|=1, the code performs the neutron quasiparticle blocking calculations in the case when no symmetries are conserved, see Section VI-2.7. For IDSIZN=+1 or −1, the blocked quasiparticle state is selected by having the largest overlap with the INSIZN-th neutron singleparticle eigenstate of the HFB mean-field Routhian or with its time-reversed partner, respectively. Note that for rotating states, the time-reversed eigenstate is not necessarily an eigenstate of the Routhian. |IDSIZN|=1 requires ISIMPY=0, IPARTY=0, IPAHFB=1, and IROTAT=1.
2.4.6. Keyword: BLOCKSIZ P 1, 0 = INSIZP, IDSIZP Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and proton quasiparticles can be simultaneously blocked.

Configurations with conserved simplex
This section lists keywords pertaining to definitions of configurations in the case when simplex is conserved, and the single-particle or single-quasiparticle space is divided into four chargesimplex blocks (neutrons and protons with simplexes s=±i). Numbers of lowest neutron states occupied in the two simplex blocks, denoted by (+) and (−), of given simplexes, s=+i and s=−i, respectively. These numbers define the simplex reference configuration from which the particle-hole excitations are counted. The definitions of simplex reference configuration and excitations are ignored unless ISIMPY=1, ISIGNY=0, and IPAIRI=0. Same as above but for the numbers of proton states.

Keyword: PHSIMP NEU
1, 00, 00, 00, 00 = NUPAHO, KPSIMP, KPSIMM, KHSIMP, KHSIMM Neutron particle-hole excitations in the simplex blocks. NUPAHO is the consecutive number from 1 to 5 (up to five sets of excitations can be specified in separate items). Particles are removed from the KHSIMP-th state in the (+) block and from the KHSIMM-th state in the (−) block, and put in the KPSIMP-th state in the (+) block and in the KPSIMM-th state in the (−) block. At every stage of constructing excitations the Pauli exclusion principle has to be respected (particle removed from an occupied state and put in an empty state). Values equal zero have no effect. In practice, reasonable excitations can only be constructed by consulting the printed lists of single-particle states with their consecutive numbers in given blocks. These particle-hole excitations are ignored unless ISIMPY=1, ISIGNY=0, and IPAIRI=0.  The diabatic blocking of neutron single-particle simplex configurations. Matrices KPFLIM contain the indices of the particle states in the two simplex blocks denoted by (+) and (−), of given simplex values, i.e., s=+i and −i, respectively. Matrices KHFLIM contain analogous indices of the hole states, and matrices KOFLIM define types of blocking in analogy to those shown in Table 4: Same as above but for the diabatic blocking of proton single-particle simplex configurations.

Configurations with conserved parity
This section lists keywords pertaining to definitions of configurations in the case when parity is conserved, and the single-particle or single-quasiparticle space is divided into four charge-parity blocks (neutrons and protons with parities π=±1). Numbers of lowest neutron states occupied in the two parity blocks, denoted by (+1) and (−1), of given parities, π=+1 and π=−1, respectively. These numbers define the parity reference configuration from which the particle-hole excitations are counted. The definitions of parity reference configuration and excitations are ignored unless IPARTY=1, or IPARTY=−1 and ISIMPY=ISIGNY=1. They are also ignored for ISIMPY=1 and IPAIRI=1. Same as above but for the numbers of proton states.

Keyword: PHPARI NEU
1, 00, 00, 00, 00 = NUPAHO, KPSIQP, KPSIQM, KHSIQP, KHSIQM Neutron particle-hole excitations in the parity blocks. NUPAHO is the consecutive number from 1 to 5 (up to five sets of excitations can be specified in separate items). Particles are removed from the KHSIQP-th state in the (+1) block and from the KHSIQM-th state in the (−1) block, and put in the KPSIQP-th state in the (+1) block and in the KPSIQM-th state in the (−1) block. At every stage of constructing excitations the Pauli exclusion principle has to be respected (particle removed from an occupied state and put in an empty state). Values equal zero have no effect. These particle-hole excitations are ignored unless ISIMPY=0, IPARTY=1, and IPAIRI=0.  The diabatic blocking of neutron single-particle parity configurations. Matrices KPFLIQ contain the indices of the particle states in the two parity blocks denoted by (+1) and (−1), of given parities, i.e., π=+1 and −1, respectively. Matrices KHFLIQ contain analogous indices of the hole states, and matrices KOFLIQ define the type of blocking according to definitions shown in Table 4:

KOFLIQ=0
⇐⇒ No diabatic blocking in the given parity block. KOFLIQ=+1 ⇐⇒ The state which has the larger y-alignment is occupied. KOFLIQ=−1 ⇐⇒ The state which has the smaller y-alignment is occupied. KOFLIQ=+2 ⇐⇒ The state which has the larger y-intrinsic spin is occupied. KOFLIQ=−2 ⇐⇒ The state which has the smaller y-intrinsic spin is occupied. KOFLIQ=+3 ⇐⇒ The state which has the larger x-alignment is occupied. KOFLIQ=−3 ⇐⇒ The state which has the smaller x-alignment is occupied. KOFLIQ=+4 ⇐⇒ The state which has the larger x-intrinsic spin is occupied. KOFLIQ=−4 ⇐⇒ The state which has the smaller x-intrinsic spin is occupied. KOFLIQ=+5 ⇐⇒ The state which has the larger z-alignment is occupied. KOFLIQ=−5 ⇐⇒ The state which has the smaller z-alignment is occupied. KOFLIQ=+6 ⇐⇒ The state which has the larger z-intrinsic spin is occupied. KOFLIQ=−6 ⇐⇒ The state which has the smaller z-intrinsic spin is occupied. KOFLIQ=+7 ⇐⇒ The state which has the larger o-alignment is occupied. KOFLIQ=−7 ⇐⇒ The state which has the smaller o-alignment is occupied. KOFLIQ=+8 ⇐⇒ The state which has the larger o-intrinsic spin is occupied. KOFLIQ=−8 ⇐⇒ The state which has the smaller o-intrinsic spin is occupied. KOFLIQ=+9 ⇐⇒ The state which has the larger multipole moment Q 20 is occupied. KOFLIQ=−9 ⇐⇒ The state which has the smaller multipole moment Q 20 is occupied. KOFLIQ=+10 ⇐⇒ The state which has the larger multipole moment Q 22 is occupied. KOFLIQ=−10 ⇐⇒ The state which has the smaller multipole moment Q 22 is occupied.
Here, the x-, y-, and z-alignments or intrinsic spins denote projections of the total angular momentum or spin, respectively, on the x, y, and z axes. Similarly, the o-alignment or o-intrinsic spin denotes analogous projections on the direction of the angular frequency ω J .
Within the diabatic blocking procedure one does not predefine whether the particle or the hole state is occupied (as is the case when the particle-hole excitations are defined, see Section 3.4 of II). In each iteration the code calculates the average alignments (or average intrinsic spins, or average quadrupole moments) of both states (those defined by KPFLIM and KHFLIM), and occupies that state for which a larger, or a smaller value is obtained. Therefore, the order of both states in the Routhian spectrum is irrelevant.
The user is responsible for choosing the particle-state indices (in KPFLIM) only among those corresponding to empty single-particle states, and the hole-state indices (in KHFLIM) only among those corresponding to occupied single-particle states. Same as above but for the diabatic blocking of proton single-particle parity configurations.

Configurations with conserved parity and signature
This section lists keywords pertaining to definitions of configurations in the case when parity and signature are conserved, and the single-particle or single-quasiparticle space is divided into eight charge-parity-signature blocks (neutrons and protons with parities π=±1 and signatures r=±i).  Table 4. Other rules described for the parity case apply here analogously. Same as above but for the diabatic blocking of proton single-particle parity-signature configurations.

Configurations
This section lists keywords pertaining to definitions of configurations in the case when any combination of parity and signature symmetry is conserved.
For IFIBLN=1, the neutron quasiparticle blocking is based on calculating overlaps with a fixed single-particle wave function. The method is based on beginning the iteration by defining the number of the single-particle state, that is , INSIGN, INSIMN, INSIQN, or INSIZN, depending on the selected symmetry conditions. Then, for INIBLN=1 this one single-particle wave function is stored in memory and on the record file, and in consecutive iterations the overlaps are calculated with respect to this fixed single-particle wave function. Therefore, subsequent changes in the ordering and structure of single-particle states do not affect the blocking mechanism. For INIBLN=0, in the first iteration this fixed single-particle wave function is not initialized, but it is read from the record file. IFIBLN=1 requires |IDSIGN|=1, |IDSIMN|=1, |IDSIQN|=1, or |IDSIZN|=1, depending on the selected symmetry conditions. IFIBLN=1 and INIBLN=0 requires ICONTI=1.

Miscellaneous parameters
This section lists keywords pertaining to definitions of parameters specifying various numerical conditions of the code execution. The standard prescription to calculate the HF potential in the next iteration is to mix a given fraction ǫ of the HF potentials from the previous iteration, with the fraction 1−ǫ of potentials calculated in the current iteration. SLOWEV and SLOWOD give the values of the mixing fractions ǫ, separately for the time-even and time-odd potentials. This mixing slows down the convergence, and is necessary for the convergence of the iterative procedure. For IBROYD=1, the Broyden method is used to accelerate the convergence, see Section VI-2.8. For IBROYD=0, or IBROYD=1 and N ITER=0, variables N ITER, ALPHAM, and BROTRI, which are read on this line, are ignored, that is, the previous values are kept. N ITER is the number of iterations used to approximate the inverse Jacobian, ALPHAM is the value of the parameter α of the linear mixing to which is added the Broyden correction, and BROTRI triggers an automatic switch to the Broyden method, that is, when the absolute value of the stability energy becomes lower than BROTRI, iterations are changed to the Broyden scheme. A large value of BROTRI ensures that the Broyden method is used from the very first iteration. IBROYD=1 is incompatible with I YUKA>0, ICOUDI=2, or ICOUEX=2. The Broyden method is implemented only in the FORTRAN-90 version of hfodd.

Parameters of the HO basis
This section lists keywords pertaining to definitions of parameters of the HO basis. The HO basis is composed of states having not more than N 0 =NOSCIL quanta in either of the Cartesian directions, and not more than M=NLIMIT states in total. The states are added to the basis according to the increasing energy of the deformed harmonic oscillator, ǫ nxnynz =hω x (n x + 1 2 ) +hω y (n y + 1 2 ) +hω z (n z + 1 2 ).
In case of degenerate HO states (e.g., for an axially deformed HO) the complete multiplets are included, so the actual number of states LDBASE can be slightly larger than NLIMIT. If NLIMIT=0, it is set to the number of states up to NOSCIL spherical HO shells, that is, NLIMIT = ((NOSCIL+1)*(NOSCIL+2)*(NOSCIL+3))/6. If NLIMIT≥0, the cutoff energy ENECUT is ignored; otherwise all states having HO energy smaller then ENECUT are included in the basis, and NLIMIT is ignored.

= FCHOM0
The code uses the standard value of the spherical HO frequencyhω 0 =41 MeV/A 1/3 multiplied by the scaling factor f =FCHOM0. The code hfodd calculates parameters of the HO basis, and the zero-iteration Nilsson potential, by defining the standard nuclear surface Σ, corresponding to the volume 4 3 π*R0PARM 3 *(INNUMB+IZNUMB), where c(α) is a function of α λµ such that the volume enclosed by the surface Σ does not depend on α, i.e., in particular c(0)=R 0 =R0PARM*(INNUMB+IZNUMB) 1/3 . The code defines frequencies of the deformed HO in three directions by using relation ω x R x = ω y R y = ω z R z , where R µ are the lengths of principal axes of the nuclear surface Σ (13) defined by real deformation parameters α λµ =ALPHAR(LAMBDA,MIU). The overall factor is defined by (ω x ω y ω z ) 1/3 =ω 0 .

Constraints on angular momentum
This section lists keywords pertaining to definitions of constraints on angular momentum. The complete cranking term in the energy density functional is given by:

= IMOVAX
For IMOVAX=1 or −1, the isoscalar angular frequency vector ω J0 is in each iteration set in the direction of the total angular momentum J 0 , or opposite to this direction, respectively, while its length is kept equal to that fixed by the sum of squares of OMEGAX, OMEGAY, and OMEGAZ. For IMOVAX=0, vector ω J0 is not changed during the iteration. IMOVAX=1 or −1 requires ITILAX=1, ITISAX=0, and IFLAGA=0. For IFLAGI(2)=1, the quadratic constraint on spin is used together with the linear constraint. Values of STIFFI(2) and ASKEDI(2) correspond respectively to C Jy andJ 0y in Eq. (14). For IFLAGI(2)=0, there is no quadratic constraint on spin. This keyword is maintained only for compatibility reasons; in version (v2.40h) of the code hfodd it has been replaced by SPICON XYZ.

Constraints on multipole moments
This section lists keywords pertaining to definitions of constraints on multipole moments and surface multipole moments. The complete multipole constraint term in the energy density functional is given by: The multipole and surface multipole moments are defined as and respectively, where Y λµ are the standard spherical harmonics in the convention of Ref. [10]. For λ≤2, factors a λµ are listed in Table 5, while for λ>2 a λµ =1.  For IFLAGQ=1, the total multipole moment of the given multipolarity λ and µ is constrained. Values of LAMBDA, MIU, STIFFQ, and QASKED correspond respectively to λ, µ, C λµ , andQ λµ in Eq. (15). For IFLAGQ=0, there is no constraint in the given multipolarity. For IFLAGS=1, the total surface multipole moment of the given multipolarity λ and µ is constrained. Values of LAMBDA, MIU, STIFFS, and SASKED correspond respectively to λ, µ, C S λµ , and Q S λµ in Eq. (15). For IFLAGS=0, there is no constraint in the given multipolarity. For IFLAGS=2 or 3, only the neutron or proton, respectively, surface multipole moments are constrained. For all constrained multipolarities, values of IFLAGS must be the same. For IFLAGG=1, the total scalar multipole moment of the given multipolarity λ is constrained. Values of LAMBDA, STIFFG, and GASKED correspond, respectively, to λ, C λ , andQ λ in Eq. (VI-80). For IFLAGG=0, there is no scalar constraint in the given multipolarity. The constrained multipolarity LAMBDA cannot be equal to 0 or exceed NMUCON. For conserved parity IPARTY=1, only even moments can be constrained.

Output-file parameters
This section lists keywords pertaining to definitions of parameters specifying character of information printed on the output file. The code prints results for the first, middle, and/or last iteration if the corresponding parameters IPRSTA, IPRMID, and IPRSTO equal 1.

= INTRIP
The code prints the values of moments (multipole moments, surface multipole moments, or magnetic moments) and angular momenta in the intrinsic frame only if INTRIP=1. The code prints tables of single-particle properties for states with values of the single-particle Routhians between EMINAL and EMAXAL. No table is printed unless EMINAL ≤ EMAXAL.

= EMAXQU
The code prints tables of quasiparticle properties for the HFB states with values of the quasiparticle energies smaller than EMAXQU. No table is printed unless EMAXQU>0.
2.13.7. Keyword: PRINT AMP 0, 999, 1, 1, 1, 0, 0 ISLPRI, ISUPRI, IENPRI, ISRPRI, IMIPRI, IKEPRI, IRMPRI These parameters govern the printing of the AMP results. If not restricted by the values of the doubled angular momenta for which the calculations are performed, IPROMI and IPROMA, the results are printed only for the doubled angular momenta between ISLPRI and ISUPRI. For IENPRI=1, the AMP energies are printed. In addition, for IENPRI=2, the AMP kernels are also printed. For ISRPRI=1, the sum rules are printed and compared with the HF average values. For IMIPRI=1, the energies of the K-mixed states are printed. For IKEPRI=1 and/or IRMPRI=1, the proton reduced kernels and/or reduced matrix elements, respectively, are printed.

= NILXYZ
This option used for NILXYZ=1, 2, or 3, allows for printing the Nilsson labels defined with respect to the x, y, or z axis, respectively. This feature is useful when the axial symmetry axis of a nucleus and its angular momentum are aligned along the x or y axis and not along the z axis, for which the standard Nilsson labels are defined. This is particularly important for analyzing configurations of band heads of odd nuclei within the conserved y-signature-symmetry limit.

Files
This section lists keywords that define files used by the code.
2.14.1. Keyword: REVIEWFILE HFODD.REV = FILREV CHARACTER*68 file name of the review file. Must start at the 13-th column of the data line. The ASCII review file is written after calculating every data set (i.e., once per every EXECUTE item), provided IREVIE>0. The file is not rewound, so the results for several data sets can be accumulated in a single file. This is so provided the filename FILREV is not changed between the EXECUTE items. The file contains all relevant parameters and results of calculation in a form suitable for reading by another program. It is meant as an interface to programs which analyze and/or plot the results. The file contains sections defined by keywords (different than keywords used in the input data file described here). The detailed structure of the review file is not documented in the present write up, and can be inferred from inspecting the specimen produced by the sample run, and from the subroutine review.

= IREVIE
The review file is not written if IREVIE=0. For IREVIE=2, the table of single-particle properties is included in the review file in addition to other results, which are written for IREVIE=1.

Keyword: REPLAYFILE
HFODD.REP = FILREP CHARACTER*68 file name of the replay file. Must start at the 13-th column of the data line. The binary replay file with the name defined in FILREP must exist if ICONTI=1, and will be read. If the filenames FILREP and FILREC are identical, the replay file will be subsequently overwritten as a new record file. These feature is implemented to facilitate chaining of jobs which follow one another.

Keyword: RECORDFILE
HFODD.REC = FILREC CHARACTER*68 file name of the record file. Must start at the 13-th column of the data line. If IWRIRE=1, binary record file is written after each HF iteration. It contains complete information which allows restarting the iteration in another run of the code. To restart, one has to specify ICONTI=1 and provide the name of the file by defining FILREP. In case of the computer crash, or upon a successful completion of the given input data set, the record file contains the results of the last performed iteration. The file is always rewound before it is written, so the results for consecutive iterations do not pile up.

= IWRIRE
For IWRIRE=1, the record file is saved on disc after each iteration is completed. For IWRIRE=0, it is saved only once, after all iterations are completed. This option is useful on systems where the speed of writing large amounts of data to disc could hamper the system performance. For IWRIRE=−1, the record file is never saved, and the run cannot be later continued.  :  ICOULI ICOULO  Action  0  0  neither read nor write the Coulomb file  1  0  read, but do not write the Coulomb file  0  1  do not read, but write the Coulomb file  1 1 use automated handling of the Coulomb file The default values of ICOULI=ICOULO=0 ensure that the code hfodd (v2.40h) behaves as that in version (v1.60r). However, unless it is required by special circumstances, a use of the automated handling of the Coulomb file (ICOULI=ICOULO=1) is recommended. Within the automated mode, the code checks whether the Coulomb file exists, and whether it contains data which match the current values of the input parameters NUMCOU, NUMETA, and FURMAX. If this is the case, the code reads the data from the Coulomb file. If this is not the case, the code calculates the Coulomb auxiliary results and stores them in the Coulomb file FILCOU. In the automated mode, the user is informed by appropriate messages printed on the output file about what type of the action has been taken in the given run of the code. 2.14.11. Keyword: YUKAWASAVE 1 = IWRIYU For IWRIYU=1 and I YUKA≥2, the Yukawa file is saved on disc after each iteration is completed. The file contains the matrix elements of the Yukawa mean field. For IWRIYU=0 and I YUKA≥2, the file is saved only once, after all iterations are completed. For IWRIYU=−1, the file is never saved. The Yukawa calculations can be restarted by using the Yukawa file. Therefore, IYCONT=1 requires setting IWRIYU=0 or 1 in the run that is to be continued. Equivalently, the Yukawa calculations can be restarted by using the field file, which requires setting IWRIFI=0 or 1 in the run that is to be continued.

Keyword: REPYUKFILE
HFODD.YUP = FILYUP CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. The binary Yukawa file with the name defined in FILYUP must exist if IYCONT=1, and will be read. If the filenames FILYUP and FILYUC are identical, the Yukawa file will be subsequently overwritten as a new Yukawa file.
2.14.13. Keyword: RECYUKFILE HFODD.YUC = FILYUC CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. If IWRIYU=1, binary Yukawa file is written after each HF iteration. It contains complete information that allows restarting the Yukawa calculations in another run of the code. To restart, one has to specify IYCONT=1 and provide the name of the file by defining FILYUP.
2.14.14. Keyword: LIPKINSAVE 1 = IWRILI For IWRILI=1 and LIPKIN=1 or LIPKIP=1, the Lipkin-Nogami file is saved on disc after each iteration is completed. The file contains the matrix elements of the particle density matrices. For IWRILI=0 and LIPKIN=1 or LIPKIP=1, the file is saved only once, after all iterations are completed. For IWRILI=−1, the file is never saved. The Lipkin-Nogami calculations can be restarted by using the Lipkin-Nogami file. Therefore, ILCONT=1 requires setting IWRILI=0 or 1 in the run that is to be continued. Equivalently, the Lipkin-Nogami calculations can be restarted by using the field file, which requires setting IWRIFI=0 or 1 in the run that is to be continued.  2.14.17. Keyword: FIELD SAVE −1 = IWRIFI For IWRIFI=1, the field file, is saved on disc after each iteration is completed. The file contains the matrix elements of the mean field. For IWRIFI=0, the file is saved only once, after all iterations are completed. For IWRIFI=−1, the file is never saved. To restart calculations of the exact Coulomb exchange energy, the field file must exist. IFCONT=1 requires setting IWRIFI=0 or 1 in the run that is to be continued.
2.14.18. Keyword: FIELD OLD 0 = IWRIOL For IWRIOL=1, in the last iteration the mean fields are not updated, that is, the mixing fractions ǫ are set equal to 1, irrespective of values of the SLOWEV, SLOWOD and SLOWPA parameters. In this way, information stored on the record file corresponds to the last but one iteration. Only then, restarting of the calculation from the field file leads to a smooth continuation of iterations.

Keyword: REP FIELDS
HFODD.FIP = FILFIP CHARACTER*68 file name of the field file. Must start at the 13-th column of the data line. The binary field file with the name defined in FILFIP must exist if IFCONT=1, and will be read. If the filenames FILFIP and FILFIC are identical, the field file will be subsequently overwritten as a new field file.
2.14.20. Keyword: REC FIELDS HFODD.FIC = FILFIC CHARACTER*68 file name of the field file. Must start at the 13-th column of the data line. If IWRIFI=1, binary field file is written after each HF iteration. It contains complete information that allows restarting the calculations from the matrix elements of fields. To restart, one has to specify IFCONT=1 and provide the name of the file by defining FILFIP.

Starting the iteration
This section lists keywords that pertain to starting and restarting the code.

= ICONTI
For ICONTI=1, results stored in the replay file (written as a record file in a previous run) are used to start the iteration. The code assumes (without checking) that the parameters of the HO basis are the same as those used in the previous run. The replay file name should be provided by defining FILREP. If the previous run was done with IREAWS=1, the current run must also use IREAWS=1, and the same Woods-Saxon file must be provided. This so because the Woods-Saxon file contains not only the information about the starting potential (which is ignored for ICONTI=1) but also defines the HO basis.

= IPCONT
For IPCONT=1, results stored in the replay file are used to define the HFB pairing properties in the first iteration; otherwise values read in matrices FERINI and DELINI are used. When the HFB pairing correlations are taken into account (i.e., for IPAHFB=1), and if a smooth restart and continuation of iterations from previously stored results is required, value of IPCONT=1 must be used. IPCONT=1 is incompatible with either of IPAHFB=0 or ICONTI=0.

= IOCONT
For IOCONT=1, the angular momentum vector J 0 stored in the replay file is used to define in the first iteration the direction of the angular frequency vector ω J0 , see Eq. (14); otherwise values read in OMEGAX, OMEGAY, and OMEGAZ are used (see keyword OMEGA XYZ). IOCONT=1 requires IMOVAX=1 or −1.

= IREAWS
For IREAWS=1, the results stored by the previously performed Woods-Saxon calculation are used to start the iteration and to define the HO basis. For ICONTI=1 and IREAWS=1 the Woods-Saxon file must also be provided, and is used only to define the HO basis.  Table 2.3 of Ref. [11], while frequencies of the deformed HO are the same as those of the HO basis defined under keyword SURFAC DEF. For NILDAT=0, the Nilsson parameters read from the input file are ignored.

= IYCONT
For IYCONT=1, results stored in the Yukawa file are used to define the Yukawa fields in the first iteration; otherwise the Yukawa fields are set equal to zero. When the Yukawa fields are taken into account (i.e., for I YUKA=2), and if a smooth restart and continuation of iterations from previously stored results is required, value of IYCONT=1 must be used. IYCONT=1 is incompatible with either of I YUKA<2 or ICONTI=0.

= ILCONT
For ILCONT=1, results stored in the Lipkin-Nogami file are used to define the density matrix required for the Lipkin-Nogami calculations in the first iteration; otherwise the required density matrix is set equal to zero. For the Lipkin-Nogami calculations (i.e., for LIPKIN=1 or LIPKIP=1), and if a smooth restart and continuation of iterations from previously stored results is required, value of ILCONT=1 must be used. ILCONT=1 is incompatible with either of LIPKIN=LIPKIP=0 or IPCONT=0.

= IFCONT
For IFCONT=1, results stored in the field file are used to define the matrix elements of fields in the first iteration; otherwise the matrix elements are recalculated. For the Gaussian-expansion method used to calculate the Coulomb energy and Coulomb mean field, that is for ICOUDI=2 or ICOUEX=2, and if a smooth restart and continuation of iterations from previously stored results is required, value of IFCONT=1 must be used. IFCONT=1 is incompatible with ICONTI=0.

OUTPUT FILE
The code writes results on standard FORTRAN output file. The output file begins with the information pertaining to the general parameters of the calculation, then gives information about the starting point of the iteration, provides the convergence report, and finally contains the results calculated at the last iteration. Structure of the output file strongly depends on the type of calculation that are requested in the given run; especially on the assumed conserved and broken symmetries. As a rule, only the information relevant for the specific run is printed. Descriptions and comments printed on the output file are in most cases self-explanatory. In this section we list only those explanations that are not explicitly included in the output file.
• Section CLASSICAL NUCLEAR SURFACE lists deformation parameters used to define the nuclear surface from which the basis parameters are derived, see keyword SURFAC DEF. It also gives the basis parameters such as the OSCILLATOR FREQUENCIES: HBAROX, HBAROY, and HBAROZ corresponding tohω x ,hω y , andhω z .
• Section PHYSICAL CONSTANTS gives values ofhc in MeV fm (H BARC), ofhc/e 2 (HBCOE2), of the neutron and proton masses in MeV/c 2 (XMASSN and XMASSP), of the kinetic-energy coefficienth 2 /2m before (HBMASS) and after the center-of-mass correction (HBMRPA), of the elementary charge squared e 2 in MeV fm (ECHAR2), and of the coefficient preceding the integral in the Coulomb exchange energy (COULEX). For details see the comments in the setbas subroutine.
• Section OSCILLATOR LENGTHS, CONSTANTS, and FREQUENCIES gives the values of 1/b µ , b µ , andhω µ , respectively, which characterize the HO basis in three Cartesian directions, µ=x, y, z.
• Section BASIS CUT-OFF CONTROL PARAMETERS gives maximum numbers of the HO quanta in three directions NXMAXX, NYMAXX, and NZMAXX corresponding to N x , N y , and N z , as well as the orders of the Gauss-Hermite quadratures NXHERM, NYHERM, and NZHERM, corresponding to L x , L y , and L z , see keyword OPTI GAUSS. It also gives the number M of the HO states included in the basis as requested in the input file (NLIMIT) and as used in the calculation (LDBASE).
• Section SHAPE OF THE OSCILLATOR-BASIS DIAMOND gives the numbers of the HO quanta in a given direction for fixed numbers of the HO quanta in both remaining directions. The output is arranged in such a way that the shape of the grid of points n x n y n z is clearly visualized by projections in every of the three directions.
• Section SKYRME FORCE DEFINITION lists the name and parameters of the Skyrme force, together with parameters KETA J, KETA W, KETACM, and KETA M that define the way the given force should be used (see keyword SKYRME-STD).
• Section COEFFICIENTS DEFINING THE SKYRME FUNCTIONAL gives values of the coupling constants in the Skyrme functional. These values take into account scaling factors, which are printed in the section SCALING FACTORS, unless they are all equal to 1.
• Section CALCULATIONS WITH THE TILTED-AXIS CRANKING gives values of components of the angular frequency vector and its length. For switches IMOVAX=1 or −1 and IOCONT=0, these values are used only in the first iteration, and later ignored, because the angular frequency vector is in each iteration readjusted to be aligned or anti-aligned with the angular momentum vector.
• Sections PARITY/SIGNATURE, SIMPLEX, or PARITY CONFIGURATIONS give the vacuum and particle-hole configurations requested in the input file for the corresponding symmetry.
• Section CONVERGENCE REPORT gives the list of performed iterations. For each iteration, one line is printed showing the energy, stability energy, average values of the intrinsic quadrupole moment Q and γ deformation, total angular momentum, angular frequency, angle between the angular frequency and angular momentum vectors, and total pairing energy.
• Section SINGLE-PARTICLE PROPERTIES lists the single-particle states calculated for the Nilsson, Woods-Saxon, or Hartree-Fock Routhian operators. For every state one line is printed which gives the value of the single-particle Routhian, the consecutive numbers in the parity-signature block or simplex block, the quantum numbers [N, n z , Λ]Ω of the Nilsson asymptotic state which has the largest component in the given state, the average value of the parity operator (in %), the average values of projections of the intrinsic and total angular momenta (inh), and their ratio called the g-factor. For broken simplex symmetry, three projections of the total angular momentum and intrinsic spin are printed in the first and second line, respectively, for each single-particle state.
• Section DENSITY INTEGRALS IN THE SKYRME FUNCTIONAL gives integrals of products of densities, which appear in the Skyrme functional. Terms in the functional are identified by the acronyms described under keywords EVE SCA TS and ODD SCA TS.
• Section CONTRIBUTIONS TO ENERGY IN THE SKYRME FUNCTIONAL gives the values of various terms which appear in the Skyrme functional. These contributions are the products of the coupling constants and of the density integrals described above. The sums of time-even and time-odd contributions are also printed.
• Section EULER ANGLES OF THE PRINCIPAL-AXES FRAME gives the standard Euler angles α, β, and γ [10] in degrees, which define the orientation of the intrinsic frame (principalaxes frame) of reference with respect to the original frame.
• Section MULTIPOLE MOMENTS gives values of neutron, proton, or mass multipole moments with respect to the original frame of reference. In particular, the λ=0 moment corresponds to the number of particles. Similarly, section MULTIPOLE MOMENTS [UNITS: (10 FERMI)^LAMBDA] [INTRINSIC FRAME] gives analogous values with respect to the intrinsic frame of reference. Whenever the parity symmetry is broken, yet another analogous section gives information on values of multipole moments with respect to the center-of-mass reference frame.
• Sections ROOT-MEAN-SQUARE AND GEOMETRIC SIZES gives the rms average values of the radius and of the x, y, and z coordinates. In order to better visualize the size of the nucleus, the geometric sizes are also calculated by multiplying the rms radius by 5/3 and the rms coordinates by √ 5.
• • Section ANGULAR MOMENTA gives for ISIMPY=1 the average values of the total and intrinsic neutron, proton, and total angular momentum (inh). It also gives the corresponding values and contributions to the first moment of inertia J (1) =I/ω (inh 2 /MeV). For ISIMPY=0, this section gives values of neutron, proton, and total projections of the total angular momentum and intrinsic spin on the three Cartesian axes in the original and intrinsic frame, respectively.
• Sections NEUTRON CONFIGURATIONS and PROTON CONFIGURATIONS give a visual representation of states occupied in the parity-signature blocks, simplex blocks, or parity blocks. The lines denoted by CONF: give the configurations requested in the input data set, while those denoted by VACC: give the configurations characterizing the given HF state. By comparing the two sets one can verify whether the requested configuration has been obtained, and eventually devise a new configuration to be calculated. The consecutive numbers printed in the horizontal direction correspond to the consecutive numbers in blocks printed in the section SINGLE-PARTICLE PROPERTIES. By comparing the two sections one can effectively associate the Nilsson labels with the calculated configurations and also prepare the configuration input data.
• Section ENERGIES gives a summary of the energies calculated for the HF state. The kinetic energy, single-particle energy, and pairing energy are printed for neutrons (NEU), protons (PRO), and all particles (TOT). Then, pairing rearrangement energy P-REARR, pairing gap PAIRGAP, and Fermi energy E-FERMI are printed. The Coulomb direct energy (DIR) and Coulomb exchange energy (EXC) are printed together with their sum (TOT). The multipole (MULT), surface multipole (SURF), and cranking (SPIN) constraint energies, are printed together with the corresponding corrections energies (CORR.). Then the rearrangement energy is printed followed by the values of the Routhian and the spin-orbit and Skyrme energies, the latter two split in the time-even (EVE) and time-odd (ODD) contributions. Finally, the total energies are printed as (SP) and (FUN), respectively, while their difference is printed as the stability (STAB).
• Section CODE COMPILED WITH THE... lists the values of the most important array dimensions and switches declared in the PARAMETER statements, fixed at the compilation stage.
• Sections BOHR DEFORMATIONS (EXACT MULTIPOLE MOMENTS) give values of the Bohr deformation parameters, determined as described in Section VI-2.5.
• Section KERNELS AND AVERAGE VALUES... gives values of the diagonal norm kernels N KK (VI-5) and average AMP energies H KK /N KK for different values of I and K. Real and imaginary parts of N KK are both printed, although all imaginary parts should be equal to zero when the numerical precision is sufficiently good.
• Section SUM RULES... compares real and imaginary parts of the sum rules (VI-49) and (VI-50) with the corresponding HF average values. Results are printed for the norm kernels N I KK , which must sum up to 1. This condition is a primary test of whether a sufficient number of the angular momenta are included in the sum of Eq. (VI-48). In the example of the file ge064-a.out, the sum rule for the norm kernel NORM= 0.3248 indicates a much too small value of I max and too small numbers of the Gauss-Tchebyschev and Gauss-Legendre integration nodes.

FORTRAN SOURCE FILE
FORTRAN source of code hfodd (v2.40h) is provided in the file hfodd.f and can be modified in several places, as described in this section.

Dimensions of arrays
The code hfodd uses the arrays' dimensions declared through the PARAMETER statements. This allows changing the dimensions and adapting the size of the reserved memory to the problem being solved. Whenever too small a dimension is defined, the code aborts with a message indicating the dimension that should be increased. Substantial amount of memory is required only for arrays depending on the following PARAMETER values:

PARAMETER (NDMAIN=21)
Corresponds to the maximum number of the HO quanta. Should be larger or equal to the input parameter NOSCIL defined under keyword BASIS SIZE.

PARAMETER (NDBASE=507)
Corresponds to the maximum number of the HO basis states. Should be larger or equal to the input parameter NLIMIT defined under keyword BASIS SIZE. It should also be larger or equal to the actual size of the HO basis LDBASE, which can be larger than NLIMIT in case of degenerate HO states.

PARAMETER (NDXHRM=33,NDYHRM=33,NDZHRM=45)
Correspond to the maximum numbers of the HO Gauss-Hermite integration nodes in three Cartesian directions. Should be larger or equal to the input parameters NXHERM, NYHERM, and NZHERM, respectively, that are defined under keyword GAUSHERMIT or determined in an optimal way, see keyword OPTI GAUSS.

PARAMETER (NDSTAT=181)
Corresponds to the maximum numbers of the single-particle states or quasiparticle states, which are kept after diagonalization. Should be larger or equal to the input parameters NUMBSP(0,0) + NUMBSP(1,0) and NUMBSP(0,1) + NUMBSP(1,1), see keyword PHASESPACE. On vector machines, parameters NDBASE, NDXHRM, NDYHRM, NDZHRM, and NDSTAT should be odd integers in order to minimize the risk of bank memory conflicts. By the same token, parameter NDMAIN should be even, because it defines the array dimensions beginning with 0.

Vectorization properties
As discussed in Sec. I-4.3, the code hfodd has to operate by using seven-fold nested short loops, and this part does not perform well in a vector processor. However, it turns out that the loops can be artificially made longer in such a way that the final CPU time in a vector processor actually becomes much shorter. All the places where this trick has been applied can be identified in the source file by finding the lines beginning with CVECTOR and CSCALAR, for example: The line beginning with CSCALAR should be made active on a scalar or superscalar machine. The line immediately below the line which begins with CVECTOR should be active on a vector machine. The results of calculation do not depend on which version of loops is activated.

FORTRAN-90 version
The code hfodd version (v2.40h) is written in FORTRAN-77. However, on several platforms, several tested FORTRAN-77 compilers did not perform well with respect to the memory management. In particular, they often distributed large arrays on computer's stack and heap in a wasteful way, and then the available memory could be insufficient for the code to execute, especially for a large HO basis. Therefore, crucial parts of the code have been promoted to FORTRAN-90, with a use of an explicit memory allocation and deallocation. In the FORTRAN source code provided in the file hfodd.f, all these FORTRAN-90 features are commented out and inactive. However, very simple modifications of the source code can easily be performed to transform code hfodd to FORTRAN-90. To this end, the user should run an automated editor script which performs the following two operations: 1. Replace all strings "C F90" by spaces, and 2. Remove two lines of code (or three lines of the file) after lines marked by the string: "C IN F90 VERSION PLEASE REMOVE THE TWO FOLLOWING LINES".
The FORTRAN-90 version should be compiled with the "fixed form" non-standard option of the FORTRAN-90 compiler. A set of c-shell and ex-editor scripts is provided within the hfodd distribution file, and allows performing the above modifications of the FORTRAN source code.
The present version of the code hfodd version (v2.40h) is the last one that works under FORTRAN-77; future releases will only use FORTRAN-90 programming.

Library subroutines 4.4.1. BLAS
The code hfodd requires an implementation of the BLAS (Basic Linear Algebra Subprogams) interface for common dense vector and matrix computations. A reference implementation of these subroutines can be downloaded from http://www.netlib.org/blas/blas.tgz, but for peak performance it is recommended that an optimized version of these routines should be installed. Optimized machine-specific BLAS (such as Sun Performance Library and Intel Math Kernel Library) are available from most hardware vendors, and there are also third-party optimized implementations such as GOTO and ATLAS available for various architectures.
The BLAS subroutines are in the REAL*8/COMPLEX*16 version, and should be compiled without promoting real numbers to the double precision. On the other hand, the code hfodd itself does require compilation with an option promoting to double precision. Therefore, the code and the BLAS package should be compiled separately, and then should be linked together.

Diagonalization subroutines
The code hfodd requires an external subroutine that diagonalizes complex hermitian matrices. Version (v1.60r) (see II) has been prepared with an interface to the NAGLIB subroutine f02axe, version (v1.75r) (see III) with an interface to the LAPACK subroutine zhpev, and version (v2.08i) (see IV) with an interface to the LAPACK subroutine zhpevx. In the present version (v2.40h), all these interfaces remain supported and can be activated as described in II-IV. However, the recommended interface is now to the LAPACK subroutine zheevr, as described in this section below.
In version (v2.08i) we have implemented interface to the LAPACK subroutine zhpevx, which can be downloaded (with dependencies) from http://netlib2.cs.utk.edu/cgi-bin/netlibfiles.pl?filename=/lapack/complex16/zhpevx.f This subroutine finds not all, but only the lowest eigenvectors, and hence performs calculations in less CPU time; the gain is particularly significant for a large HO basis. Numbers of eigenvectors to be found are defined by the size of the phasespace, see keyword PHASESPACE.
Subroutine zhpevx and its dependencies are in the REAL*8/COMPLEX*16 version, and should be compiled without promoting real numbers to the double precision. Therefore, the code and the zhpevx package should be compiled separately, and then should be linked together.
In order to activate the interface to the LAPACK zhpevx subroutine, the following modifications of the code hfodd (v2.40h) have to be made: 1. Change everywhere the value of parameter I CRAY=1 into I CRAY=0.
3. If your compiler does not support undefined externals, or subroutines called with different parameters, remove calls to subroutines cgemm, f02axe, zhpev, and zheevr.
In version (v2.40h) we have implemented interface to the LAPACK subroutine zheevr, which can be downloaded (with dependencies) from http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/lapack/complex16/zheevr.f This subroutine works with unpacked matrices and hence performs calculations in less CPU time at the expense of using a larger memory. Alternatively, the entire LAPACK library can be downloaded from http://www.netlib.org/lapack/lapack.tgz Subroutine zheevr and its dependencies are in the REAL*8/COMPLEX*16 version, and should be compiled without promoting real numbers to the double precision. Therefore, the code and the zheevr package should be compiled separately, and then should be linked together.
In order to activate the interface to the LAPACK zheevr subroutine, the following modifications of the code hfodd (v2.40h) have to be made: 1. Change everywhere the value of parameter I CRAY=1 into I CRAY=0.
3. If your compiler and linker do not support undefined externals, or subroutines called with different parameters, remove calls to subroutines cgemm, f02axe, zhpev, and zhpevx.

Matrix inversion subroutines
The code hfodd requires external subroutines that invert complex matrices and calculate their determinants. These subroutines are in the REAL*8/COMPLEX*16 version, and should be compiled without promoting real numbers to the double precision. Therefore, the code and these subroutines should be compiled separately, and then should be linked together.

Printing the execution times
After successful completion of execution, the code hfodd prints the table of CPU times spent in principal subroutines, as well as the total execution time. Since implementation of the time-checking functions widely varies between compilers and platforms, in the file hfodd.f all references to these functions are commented out and inactive. However, very simple modifications of the source code can easily be performed to activate one of the following options: 1. Replace all strings "CCPUT" by spaces in order to activate a call to the FORTRAN-90 subroutine cpu time, which returns CPU execution times.
2. Replace all strings "CETIM" by spaces in order to activate a call to function etime, which returns CPU execution times.
3. Replace all strings "CUNIX" by spaces in order to activate a call to the UNIX function time, which returns wall clock execution times.
4. Replace all strings "CSAL" by spaces in order to activate a call to the SALFORD c FORTRAN subroutine dclock@, which returns wall clock execution times.
5. Replace all strings "CF32" by spaces in order to activate a call to the MICROSOFT c F32 subroutine gettim, which returns wall clock execution times.
Other time-checking functions, available to the user, can be analogously called in the subroutine cputim.

TEST RUNS
Several examples of input and output files are provided with the released versions (v1.60r), (v1.75r), (v2.08i), and (v2.40h) of the program hfodd, or added in this user's guide, as listed in Table 6. The input data files dy152-a.dat through dy152-f.dat present the valid input items of each version and the default or recommended values of the input parameters. The input data files tb151-a.dat through tb151-d.dat are distributed with version (v1.75r), and provide examples of calculations with diabatic blocking. The input data files la132-a.dat and la132-b.dat are specific for version (v2.08i) and later, and illustrate various symmetry options of the code. The files pb128-a.dat and sn120-a.dat are added in this user's guide and give examples of the spherical HF and HFB calculations, respectively. The input data files ge064-a.dat through ge064-c.dat and sn120-b.dat are specific for version (v2.40h) and later, and illustrate the AMP calculations and Lipkin-Nogami method, respectively. In particular, running each version with its default data set corresponds to identical conditions and gives the same solution. All the input items valid for a given version are also recognized and read by the subsequent ones. In general, running a later version with an input data file prepared for an earlier version leads to the same solution. However, the structure of the output file may be modified from version to version; in particular because later versions may calculate more observables.
Running the code with the default values results in performing 50 iterations for the superdeformed state in 152 Dy at the angular frequency ofhω = 0.5 MeV. The pairing correlations are not taken into account. The HO basis is defined by an axially symmetric (with respect to the z axis) nuclear surface with physical deformation typical for superdeformed bands in the A ∼ 150 region. The parity, y-signature, x-simplex T (and their products) are imposed as conserved symmetries, while the time-reversal symmetry is broken in order to describe rotation (about the y axis). The run starts from the Nilsson potential, and an appropriate weak constraint on the quadrupole moment Q 20 is used to ensure that the iterations fall into the superdeformed minimum.
The input data files dy152-a.dat through dy152-f.dat define either the default or recommended values of the input parameters for the released versions of the code. They contain either all or some of the valid input items, as described in Table 6. File dy152-a.dat contains only keyword EXECUTE, and when used as input for any version of the code it results in performing calculations with the default data set. Files dy152-b.dat, dy152-c.dat, dy152e-dat, and dy152f-dat contain all the valid input items of versions (v1.60r), (v1.75r), (v2.08i), and (v2.40h), respectively, and the input data are identical to the default values. File dy152-d.dat contains only the input items introduced in version (v1.75r), and the corresponding new parameters are set to recommended values, which differ from the default ones in the following way.
The default values are such that the code works identically to version (v1.60r). The parameters EPSITE, NULAST, EPSPNG, and NUCHAO, control the termination of the iterations. By default, they are set so that always the number of iterations specified by NOITER is performed. The recommended values define an optimal level, at which convergence, "chaotic" divergence, or "ping-pong" divergence are detected to stop the iterations earlier. The differences between the default and recommended values of the variables governing the diabatic blocking are irrelevant because the diabatic blocking is disabled by default. Also by default, the Coulomb file is neither read from nor written to disc. With the recommended values of the parameters ICOULI and ICOULO, its reading and writing is handled automatically, which saves the CPU time.
Files dy152-a.dat through dy152-f.dat serve only as suitable patterns to modify the input parameters. However, it should be a good practice to include in the input data file only those input items which modify the input parameters with respect to the default values. In this way, the danger of accidentally modifying some intricate numerical parameters of the code is minimized.

Files tb151-a.dat through tb151-d.dat
Files tb151-a.dat through tb151-d.dat provide data for calculations in which single-particle configurations corresponding to two crossing bands in 151 Tb [12]  is lower in energy, and has number 23 -+ (23rd state of π=-1 and r=+i), while the level [761]3/2(r=+i) has number 24 -+. These two levels cross in the concerned frequency interval so that athω = 0.80 MeV the situation is the opposite. This is illustrated in Figure 5, which is explained in detail below.
In the files tb151-a.dat and tb151-b.dat it is required that always the level 23 -+ or 24 -+ be occupied, respectively, which amounts to following the configuration adiabatically. If, say, the level 23 -+ is occupied, the solutions at low and high frequencies correspond to physical configurations [514]9/2(r=+i) and [761]3/2(r=+i), respectively, while in the vicinity of the crossing no solution can be found, because the so-called "ping-pong" divergence occurs. This means that the physical states [514]9/2(r=+i) and [761]3/2(r=+i) change their energy order every second iteration, so that each of them once becomes the level 23 -+. Thus, putting the neutron always on the level 23 -+ results in occupying alternately the very different physical states [514]9/2(r=+i) and [761]3/2(r=+i), which cannot lead to convergence.
The code detects the "ping-pong" divergence, stops the iterations, and prints a table that helps identify the crossing orbitals. Such a table is reproduced in Table 7 for the 23 -+ adiabatic configuration athω = 0.65 MeV. For both protons and neutrons, in each parity-signature block, and separately for the particle and hole states, the code finds the level for which the difference in the angular momentum or intrinsic spin alignment (on the y axis) between the last and last but one iteration is the largest. The indices of those levels in the parity-signature blocks, and the absolute values of the differences are printed in the upper part of the table, independently for the angular momentum and intrinsic spin. If in some symmetry block, the indices of the particle and hole states found in this way differ by one, then it is probable that these two levels cross and cause the "ping-pong" divergence. When such pairs exist, they are listed in the lower part of the table. In the given example, the pair of states 23 -+ and 24 -+ is properly recognized on the basis of both the angular momentum and intrinsic spin. The crossing orbitals, [514]9/2(r=+i) and [761]3/2(r=+i), have markedly different alignments of about −0.46h and +2.23, respectively, which corresponds to the values 2.820 and 2.807h given in the table.
If one wishes to follow, say, the configuration [761]3/2(r=+i) diabatically, one should request that, of the levels 23 -+ and 24 -+, the one with physical properties of the state [761]3/2(r=+i) be occupied. One can chose, for instance, the state with higher angular momentum alignment, independent of whether it is lower or higher in energy. In order to diabatically follow the configuration [514]9/2(r=+i), one can chose the state with lower alignment. The corresponding runs are defined in the files tb151-c.dat and tb-151-d.dat. Indeed, in each case solutions corresponding to the same physical configuration is found at each angular frequency. Moreover, the "ping-pong" divergence does not occur, even very close to the crossing frequency. Figure 5 shows the negative-parity neutron single-particle Routhians obtained from the calculations of files tb151-c.dat and tb-151-d.dat.

Files la132-a.dat and la132-b.dat
File la132-a.dat defines two consecutive runs, which lead to a planar solution [13] in the triaxially deformed nucleus 132 La forhω = 0.25 MeV. The pairing correlations are not taken into account. Since the calculations involve the tilted-axis cranking, a spherically invariant HO basis (with the three HO frequencies equal and comprising only entire HO shells) is used in order not to favor any orientation of the solution in space. The first run is performed at hω = 0 with the parity, y-signature, x-simplex T (and their products) conserved, and with the time-reversal broken in order to describe odd numbers of neutrons and protons. Due to the presence of the odd nucleons, even forhω=0 one obtains a non-zero angular momentum, which is oriented along the y axis. The run starts from the Nilsson potential, and in 94 iterations finds Table 7: Example of the output printed when the "ping-pong" divergence is found, see text.
The second run starts from the solution obtained in the first run, and performs calculations with all symmetries broken except parity, with multipole constraints released, and for the angular frequency ofhω = 0.25 MeV. The single-particle configuration specified in terms of parity blocks is obtained from the configuration for the first run simply by adding the numbers of particles in the two signature sub-blocks of each parity block. All the three components of the initial angular frequency vector are non-zero,hω x =hω y =hω z = (0.25/ √ 3) MeV, which breaks all symmetries except parity in the first iteration, thus making the calculations otherwise symmetry-unrestricted. If a symmetry was not broken in the first iteration, it would artificially remain conserved as governed by the theorem about self-consistent symmetries [11].
The Kerman-Onishi condition [14] states that the angular momentum and angular frequency vectors must be parallel in self-consistent solutions. Standard iteration, with IMOVAX=0, requires several thousand iterations for the nucleus to turn appropriately in space to achieve the alignment of these vectors below the angle of 0.01 • . On the other hand, iteration with IMOVAX=1, as shown in the example of file la132-a.dat, achieves the alignment of 0.000001 • within 92 iterations. In the solution resulting from the second run, the angular momentum vector has non-zero components only on the y and z axes of the intrinsic frame, which are the short and long axes of the triaxial mass distribution. The calculations of the file la132-a.dat are performed in two steps to get a better control over what is going on and to avoid divergence or falling into undesired minima.
The input data provided in the file la132-b.dat are similar to those of the file la132-a.dat, except that the second run is performed with the x-simplex T conserved. In this way, the angular momentum and angular frequency vectors are confined within the plane y-z, which is also a principal plane of the mass distribution. After the first run, this plane contains the short and long axes. The components of the initial angular frequency vector are set tō hω y =hω z = (0.25/ √ 2) MeV andhω x =0. Since the planar solution in question corresponds to non-equal projections of the angular momentum on the short and long axes, in order to satisfy the Kerman-Onishi condition, the nucleus still has to turn in space about the x axis. For IMOVAX=0, this process is also very slow, while for IMOVAX=1, as in file la132-b.dat, it takes again only 92 iterations. The x-simplex T -imposing solution obtained by using the file la132-b.dat differs from the symmetry-unrestricted result of the file la132-a.dat only by orientation in space; apart from that, the values of all observables are identical. This shows that in the considered self-consistent solution the x-simplex T symmetry is spontaneously conserved.

File pb208-a.dat
The input data file pb208-a.dat defines two consecutive runs to perform precise calculations for the spherical ground state of the doubly-magic nucleus 208 Pb in a large basis of 26 HO shells. The pairing correlations are not taken into account. Conservation of all the symmetries considered in the code (the group D T 2h ) is imposed. For the spherical shape, the specified parity-signature configuration corresponds to occupation of the entire magic shells, up to the proton and neutron magic numbers 82 and 126, respectively. A spherical basis is used, i.e., with all the three characteristic HO frequencies equal (hω x =hω y =hω z =hω 0 ) and comprising only the entire HO shells. The so-called physical value ofhω 0 = f × 41 MeV/A 1/3 with f = 1.2 is taken, following the default settings of the code. To speed up the calculations, the program hfodd uses a diagonalization algorithm that can reduce the computational effort by finding less eigenstates than the matrix size (the lowest ones). In the present example, in each chargesimplex block only one extra state is calculated as compared to the number of occupied states in that block, as governed by the input parameter NUMBSP. This is far less than the code's default, and such a way of proceeding is possible only because no particle-hole excitations into higher states are required. Also for the sake of saving time, a solution in a small basis of N 0 =12 HO shells is obtained first, starting from the Nilsson potential, and then serves as a good starting point for the time-consuming calculations with N 0 =26. The number of states in the spherical basis comprising the HO shells from 0 to N 0 equals M=(N 0 + 1)(N 0 + 2)(N 0 + 3)/6, and hence the numbers of the HO basis states (for each projection of the intrinsic spin) equal M=455 and M=3654, respectively for N 0 =12 and N 0 =26. Since the code hfodd can restart iterations only if the number of Gauss-Hermite integration nodes is preserved, both runs are performed with 54 points in each of the Cartesian directions, which is excessive for N 0 =12, but optimal for N 0 =26. For N 0 =26, the 54 points ensure exact integration of all matrix elements apart from those of the Coulomb interaction and density-dependent term.
The energies of 208 Pb obtained from calculations with N 0 =12 and N 0 =26 equal −1630.638 and −1635.531 MeV, respectively. The latter value is very close to the result of −1635.666 MeV obtained by using a very precise one-dimensional spherical code working in spatial coordinates [15]. Although the N 0 =12 HO basis is evidently too small for a precise determination of the energy, other observables may not be so sensitive to the HO basis size. For instance, the mass root-mean-square radii obtained for N 0 =12 and N 0 =26 equal 5.5540 and 5.5554 fm, respectively, the result of coordinate-space calculations being 5.5550 fm.

File sn120-a.dat
The input data file sn120-a.dat provides an example of the HFB calculations for the spherical ground state in 120 Sn. All symmetries considered in the code (the group D T 2h ) are imposed, except the particle-number symmetry, which is broken in order to describe the pairing correlations. A spherical basis is used. The pairing correlations are included within the full HFB method, and the zero-range density-dependent interaction given by the pairing form-factor (11) is used in the particle-particle channel. The particle and pairing densities are obtained by summing up contributions from the quasiparticle states with the equivalent-spectrum energies not exceedingǭ max =60 MeV, which is the code's default pairing window. In the pairing formfactor (11), the power α=1 and the saturation density ρ 0 =0.32 fm −3 are taken. The interaction strength, V 0 =−285.88 MeV fm 3 , was adjusted for the given cutoff energy,ǭ max , to reproduce the value of∆ N =1.245 MeV for the average neutron pairing gap in 120 Sn. The same values of α, ρ 0 and V 0 are taken for both neutrons and protons. The Fermi energy of −8 MeV and the pairing gap of 1 MeV are used in the first iteration to construct the HFB Hamiltonian, while the initial mean field is given by the Nilsson potential.
Since 120 Sn is a proton-magic nucleus, the proton pairing decreases during the iteration and vanishes in the HFB solution. As expected, one obtains a significant thickness of the neutron skin (difference between the neutron and proton root-mean-square radii), which in the present case equals 0.14 fm.
5.1.6. Files ge064-a.dat and ge064-b.dat File ge064-a.dat contains the input data that allow for determining the rotating triaxial state in the nucleus 64 Ge and then performing a 3D AMP of this state. This input data file comprises four runs of the code: (i) starting from the Nilsson initial potential, the code performs 20 iterations with constrained values of quadrupole moments Q 20 and Q 22 , such that the initial deformation and orientation of the nucleus is determined, (ii) after releasing the constraints, the triaxial minimum is converged, (iii) the rotating solution for I = 6 is converged by setting the angular frequency ofhω = 0.575 MeV, and (iv) the AMP of the rotating solution is performed. Small numbers of the Gauss-Tchebyschev (10) and Gauss-Legendre nodes (10), see Section VI-2.1, which are used in this run, do not allow for a precise AMP resolution of high angular momenta. This input data file is only meant to provide an example of a rapid calculation.
File ge064-b.dat contains the input data that allow for a correct AMP, with higher numbers of the Gauss-Tchebyschev (40) and Gauss-Legendre nodes (40), but its execution requires a CPU time which is 4 3 =64 times longer (about 200h). Alternatively, one can run in parallel 40 jobs by executing the input data files ge064-c.dat with characters NUASTA replaced by integers from 1 to 40.

File sn120-b.dat
File sn120-b.dat contains the input data that allow for determining the ground state in the nucleus 120 Sn with the Lipkin-Nogami corrections taken into account.