Nuclear Structure in the Hartree Approximation¶
See class nucleus_rmf.
Nucleus in the Hartree approximation example¶
This example uses the NL3 EOS as implemented in eos_had_rmf and computes the structure of \(^{208}\mathrm{Pb}\) using nucleus_rmf. The results are stored in table_units objects and output to HDF files. It also computes the energy of three unoccupied neutron and proton levels.
/* Example: ex_nucleus_rmf.cpp
-------------------------------------------------------------------
This example uses the NL3 interaction to compute the structure
of Lead-208
*/
#include <o2scl/nucleus_rmf.h>
#include <o2scl/test_mgr.h>
#include <o2scl/hdf_io.h>
#include <o2scl/hdf_eos_io.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_const;
using namespace o2scl_hdf;
int lead_chden_exp(std::shared_ptr<table_units<> > profiles);
int main(int argc, char *argv[]) {
cout.setf(ios::scientific);
test_mgr t;
t.set_output_level(1);
nucleus_rmf rn;
rn.set_verbose(0);
rmf_load(rn.def_rmf,"NL3");
rn.run_nucleus(82,126,3,3);
bool neutron=false;
cout << "Proton orbitals: " << endl;
cout << "Index State Energy Nodes Deg." << endl;
for(int i=0;i<rn.nlevels;i++) {
if (rn.levels[i].isospin<0.0 && neutron==false) {
cout << endl;
cout << "Neutron orbitals: " << endl;
cout << "Index State Energy Nodes Deg." << endl;
neutron=true;
}
cout.width(2);
cout << i << " ";
cout.width(12);
cout.setf(ios::left);
cout << rn.levels[i].state << " ";
cout.unsetf(ios::left);
cout << rn.levels[i].eigen << " "
<< rn.levels[i].nodes << " "
<< rn.levels[i].twojp1 << endl;
}
cout << endl;
cout << "Proton orbitals: " << endl;
cout << "Index State Energy Nodes Deg." << endl;
neutron=false;
for(int i=0;i<rn.nuolevels;i++) {
if (rn.unocc_levels[i].isospin<0.0 && neutron==false) {
cout << endl;
cout << "Neutron orbitals: " << endl;
cout << "Index State Energy Nodes Deg." << endl;
neutron=true;
}
cout.width(2);
cout << i << " ";
cout.width(12);
cout.setf(ios::left);
cout << rn.unocc_levels[i].state << " ";
cout.unsetf(ios::left);
cout << rn.unocc_levels[i].eigen << " "
<< rn.unocc_levels[i].nodes << " "
<< rn.unocc_levels[i].twojp1 << endl;
}
cout << endl;
cout << "Proton RMS radius : " << rn.rprms << " fm." << endl;
cout << "Neutron RMS radius : " << rn.rnrms << " fm." << endl;
cout << "Total energy per nucleon: " << rn.etot << " MeV." << endl;
cout << endl;
t.test_rel(-7.842551,rn.etot,1.0e-5,"Lead binding energy");
std::shared_ptr<table_units<> > profiles=rn.get_profiles();
lead_chden_exp(profiles);
// Output the results to a file
hdf_file hf;
hdf_file hf2;
hf.open_or_create("ex_nucleus_rmf_prof.o2");
hf2.open_or_create("ex_nucleus_rmf_chden.o2");
std::shared_ptr<table_units<> > charge_dens=rn.get_chden();
hdf_output(hf,*profiles,"profiles");
hdf_output(hf2,*charge_dens,"charge_densities");
hf.close();
hf2.close();
cout << "Now with IUFSU:\n" << endl;
// IUFSU requires a different initial guess
rn.ig.sigma0=1.73;
rn.ig.omega0=1.38;
rn.ig.rho0=1.0e-13;
rmf_load(rn.def_rmf,"IUFSU");
rn.run_nucleus(82,126,0,0);
cout << "Proton RMS radius : " << rn.rprms << " fm." << endl;
cout << "Neutron RMS radius : " << rn.rnrms << " fm." << endl;
cout << "Total energy per nucleon: " << rn.etot << " MeV." << endl;
t.report();
return 0;
}
// End of example
int lead_chden_exp(std::shared_ptr<table_units<> > profiles) {
// The experimental charge density for Lead 208 from de Vries, et
// al. At. Data Nucl. Data Tables 36 (1987) 495 using the
// Sum-of-Gaussians method
double rr[12]={0.1,0.7,1.6,2.1,2.7,3.5,4.2,
5.1,6.0,6.6,7.6,8.7};
double qq[12]={0.003845,0.009724,0.033093,0.000120,
0.083107,0.080869,0.139957,0.260892,
0.336013,0.0033637,0.018729,0.000020};
double g=1.7/sqrt(1.5);
double a[12];
for(size_t i=0;i<12;i++) {
a[i]=82.0*qq[i]/2.0/pow(o2scl_const::pi,1.5)/
pow(g,3.0)/(1.0+2.0*rr[i]*rr[i]/g/g);
}
// Add experimental profile to table
profiles->new_column("chden_exp");
for(size_t i=0;i<profiles->get_nlines();i++) {
double val=0.0;
for(size_t j=0;j<12;j++) {
val+=a[j]*(exp(-pow((profiles->get("r",i)-rr[j])/g,2.0))+
exp(-pow((profiles->get("r",i)+rr[j])/g,2.0)));
}
profiles->set("chden_exp",i,val);
}
return 0;
}
Typical output:
Proton orbitals:
Index State Energy Nodes Deg.
0 ^{1}s_{1/2} -4.848748e+01 0 2
1 ^{1}p_{3/2} -4.288159e+01 0 4
2 ^{1}p_{1/2} -4.220057e+01 0 2
3 ^{1}d_{5/2} -3.586023e+01 0 6
4 ^{1}d_{3/2} -3.429919e+01 0 4
5 ^{2}s_{1/2} -3.054878e+01 1 2
6 ^{1}f_{7/2} -2.786126e+01 0 8
7 ^{1}f_{5/2} -2.508489e+01 0 6
8 ^{2}p_{3/2} -2.058053e+01 1 4
9 ^{2}p_{1/2} -1.951978e+01 1 2
10 ^{1}g_{9/2} -1.922203e+01 0 10
11 ^{1}g_{7/2} -1.501610e+01 0 8
12 ^{2}d_{5/2} -1.081390e+01 1 6
13 ^{1}h_{11/2} -1.020303e+01 0 12
14 ^{2}d_{3/2} -9.184272e+00 1 4
15 ^{3}s_{1/2} -8.054054e+00 2 2
Neutron orbitals:
Index State Energy Nodes Deg.
16 ^{1}s_{1/2} -5.903999e+01 0 2
17 ^{1}p_{3/2} -5.283626e+01 0 4
18 ^{1}p_{1/2} -5.224434e+01 0 2
19 ^{1}d_{5/2} -4.536117e+01 0 6
20 ^{1}d_{3/2} -4.394391e+01 0 4
21 ^{2}s_{1/2} -4.083662e+01 1 2
22 ^{1}f_{7/2} -3.699136e+01 0 8
23 ^{1}f_{5/2} -3.438846e+01 0 6
24 ^{2}p_{3/2} -3.046551e+01 1 4
25 ^{2}p_{1/2} -2.940207e+01 1 2
26 ^{1}g_{9/2} -2.806148e+01 0 10
27 ^{1}g_{7/2} -2.402928e+01 0 8
28 ^{2}d_{5/2} -2.046789e+01 1 6
29 ^{1}h_{11/2} -1.884691e+01 0 12
30 ^{2}d_{3/2} -1.881771e+01 1 4
31 ^{3}s_{1/2} -1.807312e+01 2 2
32 ^{1}h_{9/2} -1.335125e+01 0 10
33 ^{2}f_{7/2} -1.100552e+01 1 8
34 ^{1}i_{13/2} -9.576804e+00 0 14
35 ^{2}f_{5/2} -8.991761e+00 1 6
36 ^{3}p_{3/2} -8.290647e+00 2 4
37 ^{3}p_{1/2} -7.515362e+00 2 2
Proton orbitals:
Index State Energy Nodes Deg.
0 ^{3}s_{1/2} -8.057079e+00 2 2
1 ^{1}h_{9/2} -4.550103e+00 0 10
2 ^{2}f_{7/2} -1.365592e+00 1 8
Neutron orbitals:
Index State Energy Nodes Deg.
3 ^{1}i_{11/2} -2.912219e+00 0 12
4 ^{1}j_{15/2} -5.639563e-01 0 16
5 ^{3}d_{5/2} -1.097961e+00 2 6
Proton RMS radius : 5.460268e+00 fm.
Neutron RMS radius : 5.740404e+00 fm.
Total energy per nucleon: -7.842549e+00 MeV.
Now with IUFSU:
Proton RMS radius : 5.429696e+00 fm.
Neutron RMS radius : 5.590974e+00 fm.
Total energy per nucleon: -7.865092e+00 MeV.
1 tests performed.
All tests passed.