Optimizing the Computation of Many-Pair Density Matrix in VDPC

Jin Li¹, Wanli Ma^2*

¹Department of Physics, University of Shanghai for Science and Technology, Shanghai, China.
²Department of Mathematics, University of Shanghai for Science and Technology, Shanghai, China.

DOI: 10.4236/jamp.2022.107157                

Abstract

Recently an algorithm that acts the variational principle directly to a coherent-pair condensate (VDPC) has been proposed. This algorithm can avoid time-consuming projection while maintaining particle number conservation. Quickly computation of many-pair density matrix (MPDM) is one of the keys to improve the computational efficiency of VDPC algorithm. In this work, we propose a scheme that limits the energy range of block particles to the vicinity of the Fermi surface, which reduces the time complexity of computing the MPDM without losing physical details. The results show that by appropriately limiting the energy range, we can greatly reduce the number of matrix elements that need to be computed, and reducing the time required for the computation.

Keywords

Nuclear Structure, Mean-Field Theory, BCS, Density Matrix, VDPC

1. Introduction

Mean-field theory is one of the most widely used theories for the nuclear structure [1]. In order to introduce pairing correlation [2] [3] [4] between nuclei into mean-field theory, it is common to introduce quasiparticles, such as Hartree-Fock (HF) [5] [6] [7] [8] + Bardeen-Cooper-Schrieffer (BCS) [9] and Hartree-Fock-Bogoliubov (HFB) methods [10]. However, the drawback of the BCS and HFB methods is that since the BCS (or HFB) state is not an eigenstate of the particle number operator, they can only guarantee the average particle number conservation, but not the exact particle number conservation [4].

To solve this problem, one needs to project many-body states onto good particle number. [11] depending on the order of projection and variation, there are two ways, projection after variation (PAV) and the variation after projection (VAP). PAV has the advantage of restoring the correct particle number [12] and being more efficient, but it cannot give the minimum energy. VAP can give the optimal energy [13] [14] [15] [16], but the method becomes very time consuming because it requires multiple execution of numerical projection by integration [17].

Recently, [18] and [19] has proposed a method that applies the variational principle directly to coherent-pair condensate (VDPC). This method does not require performing the time-consuming numerical projection and also maintains particle number conservation. The coherent-pair condensate is the result of projecting the quasiparticle vacuum onto good particle number [20]. In order to efficiently calculate the energy or other observables thereof, an efficient calculation of the many-body density matrix (MPDM) of the coherent-pair condensate is an important part of the process.

The MPDM is actually equivalent to the reduced density matrix of coherent-pair condensate states. The coherent-pair condensate is essentially a many-particle wave function that introduces pairwise interactions between particles, and it is sometimes referred to as the number projected BCS wavefunction in nuclear and condensed matter physics, or the antisymmetrized geminal power (AGP) wave function in quantum chemistry. Thus, in addition to describing nucleonic pairings in nuclei, it can be applied to a variety of systems with significant pairing interactions, such as bond-breaking processes, transition metal complexes, and superconductors. And the efficient calculation of MPDM can help to develop the reduced density matrix functional theory (RDMFT) for this class of pairing systems.

The purpose of this work is to optimize the computation of the 3-order MPDM in the VDPC algorithm. The MPDM can be expressed in the form of normalized factor in the block subspace and calculated by the recurrence relations between the normalized factors [21] [22]. There is a one-to-one correspondence between the MPDM elements and the normalization factors, and the two differ by only a coefficient. The number of normalization factors directly depends on the number of its block indices (3 block indices are required in VDPC) and the value space of each block index. In general, the space of block index is almost unrestricted. However, the vast majority of physical phenomena actually occur only in a relatively small range around the Fermi surface. This means that restricting the block indicators to the vicinity of the Fermi surface will not have a large impact on the calculated values of the energy.

In this work, we discuss the number of normalization factors that need to be computed when the block index of the normalization factors is restricted. Since the time complexity of computing each normalization factor is O(1), the total number of normalization factors to be computed is the time complexity of the entire algorithm. In order to reduce the computational effort while capturing most of the physical details, we restrict each of the three block indices of the normalization factors to a different range around the Fermi surface. Then we compare our results with the method in [19], the results show that an appropriate selection of the range of block indices can greatly reduce the number of normalization factors to be computed, thus reducing the computation time of the MPDM.

The manuscript is organized as follows. In Section 2, we briefly review the analytic form of the MDPM and the normalized factors in the block subspace. In Section 3 we discuss how to restrict the block metric of the normalized factors and the time complexity of computing the MPDM under this restriction. In Section 4, we give a brief summary and outlook of the work.

2. Analytical Expressions

For an N-pair system containing 2N particles, the ground state wave function can be expressed as

$|{\varphi }_N\rangle =\frac{1}{\sqrt{{\chi }_N}}{\left(P^{†}\right)}^N|0\rangle \mathrm{,}$ (1)

where

${\chi }_N=\langle 0|P^N{\left(P^{†}\right)}^N|0\rangle$ (2)

is the normalization factor, and $|0\rangle$ is vacuum state. The coherent pair-creation operator is

$P^{†}={\displaystyle \underset{p=1}{\overset{M}{\sum }}}{\eta }_p{P}_p^{†},$ (3)

where

$P_p^{†}=c_p^{†}{c}_{\bar{p}}^{†}$ (4)

is the pair creation operator that generates a pair particles on p and $\bar{p}$ orbitals. And ${\eta }_p$ is the expansion coefficient.

[21] introduced the many-pair density matrix (MPDM) $t_{{\alpha }_1{\alpha }_2\cdots {\alpha }_p\mathrm{<mo>\; ;\; </mo>}{\beta }_1{\beta }_2\cdots {\beta }_q}^{\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]\mathrm{,}M}$, and [22] later expressed it in the form of normalization factor

$\begin{array}{c}{t}_{{\alpha }_1{\alpha }_2\cdots {\alpha }_p\mathrm{<mo>\; ;\; </mo>}{\beta }_1{\beta }_2\cdots {\beta }_q}^{\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]\mathrm{,}M}\equiv \langle 0|{\left(P\right)}^{M-p}{P}_{{\gamma }_1}{P}_{{\gamma }_2}\cdots P_{{\gamma }_r}{P}_{{\alpha }_1}{P}_{{\alpha }_2}\cdots P_{{\alpha }_p}\\ \times P_{{\beta }_1}^{†}{P}_{{\beta }_2}^{†}\cdots P_{{\beta }_q}^{†}{P}_{{\gamma }_1}^{†}{P}_{{\gamma }_2}^{†}\cdots P_{{\gamma }_r}^{†}{\left(P^{†}\right)}^{M-q}|0\rangle \\ =\frac{\left(M-p\right)!\left(M-q\right)!}{{\left[\left(M-p-q\right)!\right]}^2}{\eta }_{{\alpha }_1}{\eta }_{{\alpha }_2}\cdots {\eta }_{{\alpha }_p}{\eta }_{{\beta }_1}{\eta }_{{\beta }_2}\cdots {\eta }_{{\beta }_q}\\ \times {\chi }_{M-p-q}^{\left[{\alpha }_1{\alpha }_2\cdots {\alpha }_p{\beta }_1{\beta }_2\cdots {\beta }_q{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]}.\end{array}$ (5)

The total number of pair creation operators in the formula is $M+r$. We use ${\gamma }_1\mathrm{,}{\gamma }_2\mathrm{,}\cdots \mathrm{,}{\gamma }_r$ to mark those creation and annihilation operators with the same subscript. And there will be no identical subscript between ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\cdots \mathrm{,}{\alpha }_p$ and ${\beta }_1\mathrm{,}{\beta }_2\mathrm{,}\cdots \mathrm{,}{\beta }_q$. Due to the restriction of the Pauli principle, there will be no identical indicators in ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\cdots \mathrm{,}{\alpha }_p$, otherwise the matrix element will be zero. The same is true for ${\beta }_1\mathrm{,}{\beta }_2\mathrm{,}\cdots \mathrm{,}{\beta }_q$ and ${\gamma }_1\mathrm{,}{\gamma }_2\mathrm{,}\cdots \mathrm{,}{\gamma }_r$. The result is that there will be no identical indicator in ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\cdots \mathrm{,}{\alpha }_p\mathrm{,}{\beta }_1\mathrm{,}{\beta }_2\mathrm{,}\cdots \mathrm{,}{\beta }_q\mathrm{,}{\gamma }_1\mathrm{,}{\gamma }_2\mathrm{,}\cdots \mathrm{,}{\gamma }_r$.

The definition of ${\chi }_N^{\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]}$ is as follows,

${\chi }_N^{\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]}=\langle 0|{\left(P\right)}^N{P}_{{\gamma }_1}{P}_{{\gamma }_2}\cdots P_{{\gamma }_r}{P}_{{\gamma }_1}^{†}{P}_{{\gamma }_2}^{†}\cdots P_{{\gamma }_r}^{†}{\left(P^{†}\right)}^N|0\rangle \mathrm{.}$ (6)

Because of the operation of $P_{\gamma }$ and $P_{\gamma }^{†}$, $\gamma$ orbits are inactive in the whole system, and we call these orbits block orbits. Thus we can interpret Equation (5) as the amplitude in the block subspace $\left\{\left[{\gamma }_1{\gamma }_2\cdots {\gamma }_r\right]\right\}$, and Equation (6) is called the normalization factor in this subspace.

[18] states that the normalization factor can be calculated quickly by the following recurrence relation,

${\chi }_N=N{\displaystyle \underset{\alpha }{\sum }}{\left({\eta }_{\alpha }\right)}^2{\chi }_{N-1}^{\left[\alpha \right]},$ (7)

${\chi }_N-{\chi }_N^{\left[\alpha \right]}={\left(N{\eta }_{\alpha }\right)}^2{\chi }_{N-1}^{\left[\alpha \right]}.$ (8)

By using the initial value ${\chi }_0^{\left[\alpha \right]=1}$, we can derive ${\chi }_1$ with Equation (7), and then we can derive ${\chi }_1^{\left[\alpha \right]}$ with Equation (8). Repeating the above iterations we can obtain ${\chi }_1\cdots {\chi }_N$ and ${\chi }_1^{\left[\alpha \right]}\cdots {\chi }_N^{\left[\alpha \right]}$.

The above two equations can be written in the $\left[\beta \right]$ subspace as

${\chi }_N^{\left[\beta \right]}=N{\displaystyle \underset{\alpha \ne \beta }{\sum }}{\left({\eta }_{\alpha }\right)}^2{\chi }_{N-1}^{\left[\alpha \beta \right]},$ (9)

${\chi }_N^{\beta }-{\chi }_N^{\left[\alpha \beta \right]}={\left(N{\eta }_{\alpha }\right)}^2{\chi }_{N-1}^{\left[\alpha \beta \right]}.$ (10)

And using Equation (10), it is easy to derive

${\chi }_N^{\left[\alpha \beta \right]}=\frac{{\left({\eta }_{\alpha }\right)}^2{\chi }_N^{\left[\alpha \right]}-{\left({\eta }_{\beta }\right)}^2{\chi }_N^{\left[\beta \right]}}{{\left({\eta }_{\alpha }\right)}^2-{\left({\eta }_{\beta }\right)}^2}.$ (11)

Similarly, in the $\left[\beta \gamma \right]$ subspace, we can also write a decomposition formula like the above equation. The difference is that in the $\left[\beta \gamma \right]$ subspace, we can have various ways to construct ${\chi }_N^{\left[\alpha \beta \gamma \right]}$, as follows

$\begin{array}{c}{\chi }_N^{\left[\alpha \beta \gamma \right]}=\frac{{\left({\eta }_{\alpha }\right)}^2{\chi }_N^{\left[\alpha \gamma \right]}-{\left({\eta }_{\beta }\right)}^2{\chi }_N^{\left[\beta \gamma \right]}}{{\left({\eta }_{\alpha }\right)}^2-{\left({\eta }_{\beta }\right)}^2}\\ =\frac{{\left({\eta }_{\alpha }\right)}^2{\chi }_N^{\left[\alpha \beta \right]}-{\left({\eta }_{\gamma }\right)}^2{\chi }_N^{\left[\beta \gamma \right]}}{{\left({\eta }_{\alpha }\right)}^2-{\left({\eta }_{\gamma }\right)}^2}\\ =\frac{{\left({\eta }_{\beta }\right)}^2{\chi }_N^{\left[\alpha \beta \right]}-{\left({\eta }_{\gamma }\right)}^2{\chi }_N^{\left[\alpha \gamma \right]}}{{\left({\eta }_{\beta }\right)}^2-{\left({\eta }_{\gamma }\right)}^2}.\end{array}$ (12)

In the general case, these three decompositions are equivalent.

With all ${\chi }_N^{\left[\alpha \right]}$, we can construct each ${\chi }_N^{\left[\alpha \beta \right]}$ in $O\left(1\right)$ time complexity using Equation (11), and then construct each ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ in $O\left(1\right)$ time complexity using Equation (12). In order to make the computation efficient even in a fairly large orbit space, we would like to perform further optimization of the above method.

3. Restricted Orbit

In order to calculate energy and other observables more quickly, efficient calculation of the MPDM or normalization factor is necessary. To achieve this goal, we can restrict the block index of the normalization factor to a small region near the Fermi surface. This is because most of the physically meaningful behavior occurs near the Fermi surface, while the region far from the Fermi surface has little particle activity. But in order to preserve as much physical detail and possibility as possible (even if only qualitatively), we do not want to restrict all the orbits to a very small region.

Therefore we can choose to restrict each of the three block indices to different ranges around the Fermi surface. Specifically, we can consider restricting the values of $\alpha \mathrm{,}\beta$, and $\gamma$ to the sets ${\Omega }_1$, ${\Omega }_2$, and ${\Omega }_3$, respectively. The number of orbitals contained in ${\Omega }_1$, ${\Omega }_2$, ${\Omega }_3$ are $M_1$, $M_2$, $M_3$, respectively, and these orbitals are numbered $1~M_1$, $1~M_2$ and $1~M_3$. Where ${\Omega }_1$, and ${\Omega }_2$ have orbitals whose energies are restricted to some energy range centered on the Fermi surface, while ${\Omega }_3=\Omega$ denotes the set of orbitals contained in the full single-particle orbit space. And ${\Omega }_1\mathrm{,}{\Omega }_2\mathrm{,}{\Omega }_3$ satisfy condition

${\Omega }_1\subset {\Omega }_2\subset {\Omega }_3=\Omega .$ (13)

Therefore the number of orbits satisfies

$M_1<M_2<M_3.$ (14)

Such a restriction method can greatly reduce the number of normalization factor that need to be computed, while the computation using these normalization factors can minimize the loss of physical details. And in the case of that restriction method, we would like to discuss exactly how much the number of computations is reduced.

Note here that since the normalized factor $\chi$ are symmetric with respect to the arrangement of the superscripts, for convenience we assume in the subsequent discussion that $\alpha \mathrm{,}\beta \mathrm{,}\gamma$ satisfy $\alpha <\beta <\gamma$.

We first discuss the number of second-order normalization factor ${\chi }_N^{\left[\alpha \beta \right]}$ that need to be computed to constitute all third-order normalization factor ${\chi }_N^{\left[\alpha \beta \gamma \right]}$. In the case where the range of values of the orbit is restricted, the number of ${\chi }_N^{\left[\alpha \beta \right]}$ that need to be computed also depend on the way ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ are decomposed.

Specifically, since superscript of ${\chi }_N^{\left[\alpha \beta \right]}$, ${\chi }_N^{\left[\beta \gamma \right]}$ and ${\chi }_N^{\left[\alpha \gamma \right]}$ in Equation (12) are restricted to the direct product ${\Omega }_1{\Omega }_2$, ${\Omega }_2{\Omega }_3$ and ${\Omega }_1{\Omega }_3$, and ${\Omega }_1{\Omega }_2\subset {\Omega }_1{\Omega }_3\subset {\Omega }_2{\Omega }_3$, so that the second-order normalization factors required for the three decomposition methods in Equation (12) are restricted to ${\Omega }_2{\Omega }_3$, ${\Omega }_2{\Omega }_3$ and ${\Omega }_1{\Omega }_3$, respectively. Obviously, the third decomposition method only requires the least number of second-order $\chi$ to be computed. Assuming that there are $M_1$ orbits in ${\Omega }_1$, $M_2$ orbits in ${\Omega }_2$, and the number of orbits $M_3$ in ${\Omega }_3$ is equal to the number of orbits M in the full space, the number of second-order $\chi$ to be computed is

$\left(\begin{array}{c}{M}_3\\ 2\end{array}\right)-\left(\begin{array}{c}{M}_3-M_1\\ 2\end{array}\right)\mathrm{.}$ (15)

The first term indicates the number of second-order normalization factors that do not impose restrictions on the range of values of the orbits. The second term denotes the elimination of all terms with $\alpha \notin {\Omega }_1$. Thus all the remaining terms satisfy $\alpha \in {\Omega }_1$ and $\gamma \in {\Omega }_3$.

Next we discuss the number of third-order normalized factors ${\chi }_N^{\left[\alpha \beta \gamma \right]}$.

Before giving a specific form for the number of ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ to be computed, we first need to decide how we are going to pick the normalization factors with restricted orbits from those with unrestricted orbits. In principle, the results are equivalent using any of the picking methods. But the appropriate method can give some formal convenience to the results. Our picking method is clearly shown in the Figure 1. The reason for choosing this method is that it allows the form of the results to merge better with Equation (15) in the summation.

Following the classification shown in Figure 1, we can directly write that the number of ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ to be calculated under the restriction is

Figure 1. The method of selecting the normalization factors ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ that satisfy $\alpha \in {\Omega }_1$, $\beta \in {\Omega }_2$, and $\gamma \in {\Omega }_3$ from among all possible normalization factors. Note that the preset $\alpha <\beta <\gamma$. We divide all possible terms into three regions: the terms in region 1 satisfy $\alpha \in {\Omega }_1$, $\beta \in {\Omega }_2$, and $\gamma \in {\Omega }_3$ ; the terms in region 2 satisfy $\alpha \in {\Omega }_1$, $\beta \notin {\Omega }_2$ ; and the terms in region 3 satisfy $\alpha \notin {\Omega }_1$.

$\left(\begin{array}{c}{M}_3\\ 3\end{array}\right)-\left(\begin{array}{c}{M}_3-M_1\\ 3\end{array}\right)-\left(\begin{array}{c}{M}_1\\ 1\end{array}\right)\left(\begin{array}{c}{M}_3-M_2\\ 2\end{array}\right)\mathrm{.}$ (16)

The first term in the formula denotes the number of ${\chi }_N^{\left[\alpha \beta \gamma \right]}$ without restrictions. The second term indicates the exclusion of all terms where $\alpha \notin {\Omega }_1$ (region 3 in Figure 1. $\alpha ,\beta ,\gamma$ are all in ${\Omega }_3-{\Omega }_1$ ). The third term indicates the exclusion of all terms where $\alpha \in {\Omega }_1$ and $\beta \notin {\Omega }_2$ (region 2 in Figure 1. $\beta$, $\gamma$ are in ${\Omega }_3-{\Omega }_2$ ). Thus all remaining terms satisfy the condition $\alpha \in {\Omega }_1$, $\beta \in {\Omega }_2$, and $\gamma \in {\Omega }_3$.

Recall that when calculating the third-order normalized coefficients using Equations (11) and (12), the time complexity of each term is $O\left(1\right)$. This means that the computer’s calculation time is proportional to the number of computations. Therefore we can use the ratio of the number of computations required before and after imposing the restriction as a reference for the total computation time.

Combining Equations (15) and (16) easily knows that when we restrict the values of the block orbits, the number of computations is

$\begin{array}{l}\left(\begin{array}{c}{M}_3+1\\ 3\end{array}\right)-\left(\begin{array}{c}{M}_3-M_1+1\\ 3\end{array}\right)-\left(\begin{array}{c}{M}_1\\ 1\end{array}\right)\left(\begin{array}{c}{M}_3-M_2\\ 2\end{array}\right)\mathrm{,}\hfill \end{array}$ (17)

and when the value of block orbits are not restricted, the number of computations is

$\begin{array}{l}\left(\begin{array}{c}{M}_3+1\\ 3\end{array}\right)\mathrm{.}\hfill \end{array}$ (18)

Note that this number is also the number of normalization factors that need to be calculated in [19].

Figure 2 shows the ratio of Equation (17) to Equation (18) for different values of $M_1$ and $M_2$. The results show that an appropriate range of $M_1$ and $M_2$

Figure 2. The ratio of the number of normalization coefficients before and after restricting the range of orbit taking. $M_3=500$. The values of $M_2$ are taken as 200, 300, 400. The values of $M_1$ is taken at intervals of 20 in the range 1 to $M_2-19$.

Figure 3. Mean of calculation times of normalization factors when the values of the block indices are restricted. The range of values of $\alpha$, $\beta$ and $\gamma$ is $1~\left(M_3-2i\right)$, $1~\left(M_3-i\right)$ and $1~M_3$, respectively. The black line is constant c multiplied by the Equation (17), where $c=1.86\times {10}^{-6}$ is obtained by least square fitting.

values can greatly reduce the computation cost. Taking $M_2=300$ and $M_1=141$ as an example, the number of normalization factors is reduced by about 50%. Correspondingly, the total computation time should also be reduced by half due to the positive relationship between the total computation time and the number of normalization coefficients.

In Figure 3 we show the actual calculation time of the code when the values of the block indices are restricted, and then we fit the results using Equation (17). Every point in the plot is the mean of 10 calculations, and the black line denotes Equation (17) multiplied by a constant c, with c derived from by least squares fitting. The three block indices take values in the ranges $1~\left(M_3-2i\right)$, $1~\left(M_3-i\right)$, and $1~M_3$, respectively. i is taken from 10 to 50 with an interval of 5. The result shows that Equation (17) can well describe the variation of the computation time with the value range of the block indices.

4. Conclusions

In most independent particle models, the calculation of the reduced density matrix is a crucial step. However, in some models that include pairwise interactions, such as the number projected BCS method, instead of calculating the reduced density matrix, the calculation of the MPDM is performed. The calculation of MPDM is almost the most time-consuming part in large-scale computations, so efficient calculation of MPDM is very important for the development of this class of pairwise models.

Since the vast majority of physical processes occur near the Fermi surface. The many-body states far from the Fermi surface, due to their small probability amplitude, have little impact on the energy. Therefore, in many cases, the matrix elements in MPDM of states far away from the Fermi surface are not necessary for the calculation of the energy. Moreover, in the calculation of some physical quantities, only a few matrix elements are needed. In this work, we want to reduce the computational effort by excluding those states far from the Fermi surface during the calculation. This goal can be accomplished by limiting the index of the normalization factor in the block subspace corresponding to the MPDM. In addition, if only a few matrix elements in MPDM are needed in some calculations, one can quickly derive these matrix elements by selecting a suitable restricted range.

We discussed which second-order and third-order normalization factors are needed to be computed in the case where the value range of the block indices is restricted, and we discuss the number of these normalization factors. This number directly reflects the time complexity of the computation. The results show that limiting the value range of the block indices can greatly reduce the number of normalization factors to be computed. In the next work, we would like to discuss the effect of different restrictions on the energy in order to find a balance between energy accuracy and time complexity.

References

1. 1. Bender, M., Heenen, P.-H., and Reinhard, P.-G. (2003) Self-Consistent Mean-Field Models for Nuclear Structure. Reviews of Modern Physics, 75, 121-180. https://doi.org/10.1103/RevModPhys.75.121


2. 2. Bohr, A., Mottelson, B.R. and Pines, D. (1958) Possible Analogy between the Excitation Spectra of Nuclei and Those of the Superconducting Metallic State. Physical Review, 110, 936-938. https://doi.org/10.1103/PhysRev.110.936


3. 3. Bohr, A. and Mottelson, B. (1975) Nuclear Structure. Benjamin, New York.


4. 4. Ring, P. and Schuck, P. (1980) The Nuclear Many-Body Problem. Springer-Verlag, Berlin.


5. 5. Muthukrishnan, R. and Baranger, M. (1965) Nuclear Hartree-Fock Calculations. Physics Letters, 18, 160-163. https://doi.org/10.1016/0031-9163(65)90698-0


6. 6. Davies, K., Krieger, S. and Baranger, M. (1966) A Study of the Hartree-Fock Approximation as Applied to Finite Nuclei. Nuclear Physics, 84, 545-568. https://doi.org/10.1016/0029-5582(66)91013-3


7. 7. Bassichis, W., Kerman, A. and Svenne, J. (1967) Unrestricted Hartree-Fock Treatment of Finite Nuclei. Physical Review, 160, 746-752. https://doi.org/10.1103/PhysRev.160.746


8. 8. Saunier, G. and Pearson, J. (1967) Effective Central Interaction for Nuclear Hartree-Fock Calculations. Physical Review, 160, 740-745. https://doi.org/10.1103/PhysRev.160.740


9. 9. Bardeen, J., Cooper, L.N. and Schrieffer, J.R. (1957) Theory of Superconductivity. Physical Review, 108, 1175-1204. https://doi.org/10.1103/PhysRev.108.1175


10. 10. Mang, H. (1975) The Self-Consistent Single-Particle Model in Nuclear Physics. Physics Reports, 18, 325-368. https://doi.org/10.1016/0370-1573(75)90012-5


11. 11. Dietrich, K., Mang, H.J. and Pradal, J.H. (1964) Conservation of Particle Number in the Nuclear Pairing Model. Physical Review, 135, B22-B34. https://doi.org/10.1103/PhysRev.135.B22


12. 12. Anguiano, M., Egido, J. and Robledo, L. (2002) Mean-Field Based Approaches to Pairing Correlations in Atomic Nuclei. Physics Letters B, 545, 62-72. https://doi.org/10.1016/S0370-2693(02)02557-1


13. 13. Dukelsky, J. and Sierra, G. (2000) Crossover from Bulk to Few-Electron Limit in Ultrasmall Metallic Grains. Physical Review B, 61, 12302-12314. https://doi.org/10.1103/PhysRevB.61.12302


14. 14. Sandulescu, N. and Bertsch, G.F. (2008) Accuracy of BCS-Based Approximations for Pairing in Small Fermi Systems. Physical Review C, 78, Article ID: 064318. https://doi.org/10.1103/PhysRevC.78.064318


15. 15. Anguiano, M., Egido, J. and Robledo, L. (2001) Particle Number Projection with Effective Forces. Nuclear Physics A, 696, 467-493. https://doi.org/10.1016/S0375-9474(01)01219-2


16. 16. Hupin, G. and Lacroix, D. (2012) Number-Conserving Approach to the Pairing Problem: Application to Kr and Sn Isotopic Chains. Physical Review C, 86, Article ID: 24309. https://doi.org/10.1103/PhysRevC.86.024309


17. 17. Wang, X.B., Dobaczewski, J., Kortelainen, M., Yu, L.F. and Stoitsov, M.V. (2014) Lipkin Method of Particle-Number Restoration to Higher Orders. Physical Review C, 90, Article ID: 14312. https://doi.org/10.1103/PhysRevC.90.014312


18. 18. Jia, L.Y. (2019) Application of the Variational Principle to a Coherent-Pair Condensate: The BCS Case. Physical Review C, 99, Article ID: 014302. https://doi.org/10.1103/PhysRevC.99.014302


19. 19. Jia, L.Y. (2020) Computer Code for Applying the Variational Principle to a Coherent-Pair Condensate: The BCS Case. Computer Physics Communications, 248, Article ID: 106967. https://doi.org/10.1016/j.cpc.2019.106967


20. 20. Talmi, I. (1993) Simple Models of Complex Nuclei: The Shell Model and Interacting Boson Model. Harwood Academic, Chur.


21. 21. Jia, L.Y. (2015) Practical Calculation Scheme for Generalized Seniority. Journal of Physics G: Nuclear and Particle Physics, 42, Article ID: 115105. https://doi.org/10.1088/0954-3899/42/11/115105


22. 22. Jia, L.Y. (2017) Generalized Seniority on a Deformed Single-Particle Basis. Physical Review C, 96, Article ID: 034313. https://doi.org/10.1103/PhysRevC.96.034313

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Bender, M., Heenen, P.-H., and Reinhard, P.-G. (2003) Self-Consistent Mean-Field Models for Nuclear Structure. Reviews of Modern Physics, 75, 121-180. https://doi.org/10.1103/RevModPhys.75.121
[2]	Bohr, A., Mottelson, B.R. and Pines, D. (1958) Possible Analogy between the Excitation Spectra of Nuclei and Those of the Superconducting Metallic State. Physical Review, 110, 936-938. https://doi.org/10.1103/PhysRev.110.936
[3]	Bohr, A. and Mottelson, B. (1975) Nuclear Structure. Benjamin, New York.
[4]	Ring, P. and Schuck, P. (1980) The Nuclear Many-Body Problem. Springer-Verlag, Berlin.
[5]	Muthukrishnan, R. and Baranger, M. (1965) Nuclear Hartree-Fock Calculations. Physics Letters, 18, 160-163. https://doi.org/10.1016/0031-9163(65)90698-0
[6]	Davies, K., Krieger, S. and Baranger, M. (1966) A Study of the Hartree-Fock Approximation as Applied to Finite Nuclei. Nuclear Physics, 84, 545-568. https://doi.org/10.1016/0029-5582(66)91013-3
[7]	Bassichis, W., Kerman, A. and Svenne, J. (1967) Unrestricted Hartree-Fock Treatment of Finite Nuclei. Physical Review, 160, 746-752. https://doi.org/10.1103/PhysRev.160.746
[8]	Saunier, G. and Pearson, J. (1967) Effective Central Interaction for Nuclear Hartree-Fock Calculations. Physical Review, 160, 740-745. https://doi.org/10.1103/PhysRev.160.740
[9]	Bardeen, J., Cooper, L.N. and Schrieffer, J.R. (1957) Theory of Superconductivity. Physical Review, 108, 1175-1204. https://doi.org/10.1103/PhysRev.108.1175
[10]	Mang, H. (1975) The Self-Consistent Single-Particle Model in Nuclear Physics. Physics Reports, 18, 325-368. https://doi.org/10.1016/0370-1573(75)90012-5
[11]	Dietrich, K., Mang, H.J. and Pradal, J.H. (1964) Conservation of Particle Number in the Nuclear Pairing Model. Physical Review, 135, B22-B34. https://doi.org/10.1103/PhysRev.135.B22
[12]	Anguiano, M., Egido, J. and Robledo, L. (2002) Mean-Field Based Approaches to Pairing Correlations in Atomic Nuclei. Physics Letters B, 545, 62-72. https://doi.org/10.1016/S0370-2693(02)02557-1
[13]	Dukelsky, J. and Sierra, G. (2000) Crossover from Bulk to Few-Electron Limit in Ultrasmall Metallic Grains. Physical Review B, 61, 12302-12314. https://doi.org/10.1103/PhysRevB.61.12302
[14]	Sandulescu, N. and Bertsch, G.F. (2008) Accuracy of BCS-Based Approximations for Pairing in Small Fermi Systems. Physical Review C, 78, Article ID: 064318. https://doi.org/10.1103/PhysRevC.78.064318
[15]	Anguiano, M., Egido, J. and Robledo, L. (2001) Particle Number Projection with Effective Forces. Nuclear Physics A, 696, 467-493. https://doi.org/10.1016/S0375-9474(01)01219-2
[16]	Hupin, G. and Lacroix, D. (2012) Number-Conserving Approach to the Pairing Problem: Application to Kr and Sn Isotopic Chains. Physical Review C, 86, Article ID: 24309. https://doi.org/10.1103/PhysRevC.86.024309
[17]	Wang, X.B., Dobaczewski, J., Kortelainen, M., Yu, L.F. and Stoitsov, M.V. (2014) Lipkin Method of Particle-Number Restoration to Higher Orders. Physical Review C, 90, Article ID: 14312. https://doi.org/10.1103/PhysRevC.90.014312
[18]	Jia, L.Y. (2019) Application of the Variational Principle to a Coherent-Pair Condensate: The BCS Case. Physical Review C, 99, Article ID: 014302. https://doi.org/10.1103/PhysRevC.99.014302
[19]	Jia, L.Y. (2020) Computer Code for Applying the Variational Principle to a Coherent-Pair Condensate: The BCS Case. Computer Physics Communications, 248, Article ID: 106967. https://doi.org/10.1016/j.cpc.2019.106967
[20]	Talmi, I. (1993) Simple Models of Complex Nuclei: The Shell Model and Interacting Boson Model. Harwood Academic, Chur.
[21]	Jia, L.Y. (2015) Practical Calculation Scheme for Generalized Seniority. Journal of Physics G: Nuclear and Particle Physics, 42, Article ID: 115105. https://doi.org/10.1088/0954-3899/42/11/115105
[22]	Jia, L.Y. (2017) Generalized Seniority on a Deformed Single-Particle Basis. Physical Review C, 96, Article ID: 034313. https://doi.org/10.1103/PhysRevC.96.034313