Project 3: Developing a Hartree-Fock program

Here we will demonstrate how do develop a Hartree-Fock program that will allow us to compute the binding energy of nuclei like $^{16}$O or $^{40}$Ca. The single-particle energies obtained by solving the Hartree-Fock equations can then be directly related to experimental separation energies. Since Hartree-Fock theory is the starting point for several many-body techniques (density functional theory, random-phase approximation, shell-model etc), the aim here is to develop a computer program to solve the Hartree-Fock equations in a given single-particle basis, here the harmonic oscillator.

The Hamiltonian for a system of \( A \) nucleons confined in a harmonic potential reads $$ \hat{H} = \sum_{i=1}^{A} \frac{\hat{p}_{i}^{2}}{2m}+\sum_{i=1}^{A} \frac{1}{2} m\omega {r}_{i}^{2}+\sum_{i < j} \hat{V}_{ij}, $$ with \( \hbar^{2}/2m = 20.73 \) fm$^{2}$, \( mc^{2} = 938.926 \) MeV, and \( \hat{V}_{ij} \) is the two-body interaction potential whose matrix elements are precalculated and to be read in by you, see weblink below.

The Hartree-Fock algorithm can be broken down as follows. We recall that our Hartree-Fock matrix is $$ \hat{h}_{\alpha\beta}^{HF}=\langle \alpha \vert\hat{h}_0 \vert \beta \rangle+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$ Normally we assume that the single-particle basis \( \vert\beta\rangle \) forms an eigenbasis for the operator \( \hat{h}_0 \) (this is our case), meaning that the Hartree-Fock matrix becomes $$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$ The Hartree-Fock eigenvalue problem $$ \sum_{\beta}\hat{h}_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{\mathrm{HF}}C_{i\alpha}, $$ can be written out in a more compact form as $$ \hat{h}^{HF}\hat{C}=\epsilon^{\mathrm{HF}}\hat{C}. $$

The equations are often rewritten in terms of a so-called density matrix, which is defined as $$ \begin{equation} \rho_{\gamma\delta}=\sum_{i=1}^{A}\langle\gamma|i\rangle\langle i|\delta\rangle = \sum_{i=1}^{A}C_{i\gamma}C^*_{i\delta}. \tag{18} \end{equation} $$ It means that we can rewrite the Hartree-Fock Hamiltonian as $$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{\gamma\delta} \rho_{\gamma\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$ It is convenient to use the density matrix since we can precalculate in every iteration the product of two eigenvector components \( C \).

Note that \( \langle \alpha\vert\hat{h}_0\vert\beta \rangle \) denotes the matrix elements of the one-body part of the starting hamiltonian. For self-bound nuclei \( \langle \alpha\vert\hat{h}_0\vert\beta \rangle \) is the kinetic energy while \( \langle \alpha \vert \hat{h}_0 \vert \beta \rangle \) represents the harmonic oscillator hamiltonian since the system is confined in a harmonic trap. If we are working in a harmonic oscillator basis with the same \( \omega \) as the trapping potential, then \( \langle \alpha\vert\hat{h}_0 \vert \beta \rangle \) is diagonal.

An example of a function written in python which performs the Hartree-Fock calculation is shown here. In setting up your code you will need to write a function which sets up the single-particle basis, the matrix elements \( t_{\alpha\gamma} \) of the one-body operator (called \( h0 \) in the function below) and the antisymmetrized TBMEs (called nninteraction in the code link below) and the density matrix elements \( \rho_{\beta\delta} \) (called densityMatrix below). The python program shows how one can, in a brute force way read in matrix elements in \( m \)-scheme and compute the Hartree-Fock single-particle energies for four major shells. The interaction which has been used is the so-called N3LO interaction of Machleidt and Entem using the Similarity Renormalization Group approach method to renormalize the interaction, using an oscillator energy \( \hbar\omega=10 \) MeV.

The codes used to obtain the matrix elements are available at the CENS site.

The python program shows how one can, in a brute force way read in matrix elements in \( m \)-scheme and compute the Hartree-Fock single-particle energies for four major shells. The interaction which has been used is the so-called N3LO interaction of Machleidt and Entem using the Similarity Renormalization Group approach method to renormalize the interaction, using an oscillator energy \( \hbar\omega=10 \) MeV.

The nucleon-nucleon two-body matrix elements are in \( m \)-scheme and are fully anti-symmetrized. The Hartree-Fock programs uses the density matrix discussed above in order to compute the Hartree-Fock matrix. Here we display the Hartree-Fock part only, assuming that single-particle data and two-body matrix elements have already been read in.

	""" Star HF-iterations, preparing variables and density matrix """

        """ Coefficients for setting up density matrix, assuming only one along the diagonals """
	C = np.eye(spOrbitals) # HF coefficients
        DensityMatrix = np.zeros([spOrbitals,spOrbitals])
        for gamma in range(spOrbitals):
            for delta in range(spOrbitals):
                sum = 0.0
                for i in range(Nparticles):
                    sum += C[gamma][i]*C[delta][i]
                DensityMatrix[gamma][delta] = Decimal(sum)
        maxHFiter = 100
        epsilon =  1.0e-10 
        difference = 1.0
	hf_count = 0
	oldenergies = np.zeros(spOrbitals)
	newenergies = np.zeros(spOrbitals)
	while hf_count < maxHFiter and difference > epsilon:
		print "############### Iteration %i ###############" % hf_count
   	        HFmatrix = np.zeros([spOrbitals,spOrbitals])		
		for alpha in range(spOrbitals):
			for beta in range(spOrbitals):
                            """  If tests for three-dimensional systems, including isospin conservation """
                            if l[alpha] != l[beta] and j[alpha] != j[beta] and mj[alpha] != mj[beta] and tz[alpha] != tz[beta]: continue
                            """  Setting up the Fock matrix using the density matrix and antisymmetrized NN interaction in m-scheme """
     		            sumFockTerm = 0.0
                            for gamma in range(spOrbitals):
                                for delta in range(spOrbitals):
                                    if (mj[alpha]+mj[gamma]) != (mj[beta]+mj[delta]) and (tz[alpha]+tz[gamma]) != (tz[beta]+tz[delta]): continue
                                    sumFockTerm += DensityMatrix[gamma][delta]*nninteraction[alpha][gamma][beta][delta]
                            HFmatrix[alpha][beta] = Decimal(sumFockTerm)
                            """  Adding the one-body term, here plain harmonic oscillator """
                            if beta == alpha:   HFmatrix[alpha][alpha] += singleparticleH[alpha]
		spenergies, C = np.linalg.eigh(HFmatrix)
                """ Setting up new density matrix in m-scheme """
                DensityMatrix = np.zeros([spOrbitals,spOrbitals])
                for gamma in range(spOrbitals):
                    for delta in range(spOrbitals):
                        sum = 0.0
                        for i in range(Nparticles):
                            sum += C[gamma][i]*C[delta][i]
                        DensityMatrix[gamma][delta] = Decimal(sum)
		newenergies = spenergies
                """ Brute force computation of difference between previous and new sp HF energies """
                sum =0.0
                for i in range(spOrbitals):
                    sum += (abs(newenergies[i]-oldenergies[i]))/spOrbitals
                difference = sum
                oldenergies = newenergies
                print "Single-particle energies, ordering may have changed "
                for i in range(spOrbitals):
                    print('{0:4d}  {1:.4f}'.format(i, Decimal(oldenergies[i])))
		hf_count += 1

Running the program, one finds that the lowest-lying states for a nucleus like \( {}^{16}\mbox{O} \), we see that the nucleon-nucleon force brings a natural spin-orbit splitting for the \( 0p \) states (or other states except the \( s \)-states). Since we are using the \( m \)-scheme for our calculations, we observe that there are several states with the same eigenvalues. The number of eigenvalues corresponds to the degeneracy \( 2j+1 \) and is well respected in our calculations, as see from the table here.

The values of the lowest-lying states are (\( \pi \) for protons and \( \nu \) for neutrons)

We can use these results to attempt our first link with experimental data, namely to compute the shell gap or the separation energies. The shell gap for neutrons is given by $$ \Delta S_n= 2BE(N,Z)-BE(N-1,Z)-BE(N+1,Z). $$ For \( {}^{16}\mbox{O} \) we have an experimental value for the shell gap of \( 11.51 \) MeV for neutrons, while our Hartree-Fock calculations result in \( 25.65 \) MeV. This means that correlations beyond a simple Hartree-Fock calculation with a two-body force play an important role in nuclear physics. The splitting between the \( 0p_{3/2}^{\nu} \) and the \( 0p_{1/2}^{\nu} \) state is 4.88 MeV, while the experimental value for the gap between the ground state \( 1/2^{-} \) and the first excited \( 3/2^{-} \) states is 6.08 MeV. The two-nucleon spin-orbit force plays a central role here. In our discussion of nuclear forces we will see how the spin-orbit force comes into play here.

• «
• 1
• ...
• 59
• 60
• 61
• 62
• 63
• 64
• 65
• 66
• 67