Nuclear structure in the Hartree approximation

See class o2scl::nucleus_rmf .

Nucleus in the Hartree approximation example

This example uses the NL3 EOS as implemented in o2scl::eos_had_rmf and computes the structure of using o2scl::nucleus_rmf. The results are stored in o2scl::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>
#ifdef O2SCL_HDF
#include <o2scl/hdf_io.h>
#include <o2scl/hdf_eos_io.h>
#endif

using namespace std;
using namespace o2scl;
using namespace o2scl_const;
#ifdef O2SCL_HDF
using namespace o2scl_hdf;
#endif

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");

#ifdef O2SCL_HDF

  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<> > profiles=rn.get_profiles();
  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();

  //profiles.line_of_names(((string)"r rhop rhon sig ome rho ")+
  //"coul chden rhosp rhosn ");
  //chden_table.line_of_names("x chden1 chdenc");

#endif

  t.report();

  return 0;
}
// End of example

Nucleon density distributions in Pb-208

Typical output:

Proton orbitals: 
Index State      Energy       Nodes Deg.
 0 ^{1}s_{1/2}  -4.848749e+01 0     2
 1 ^{1}p_{3/2}  -4.288160e+01 0     4
 2 ^{1}p_{1/2}  -4.220058e+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.508490e+01 0     6
 8 ^{2}p_{3/2}  -2.058053e+01 1     4
 9 ^{2}p_{1/2}  -1.951979e+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.184274e+00 1     4
15 ^{3}s_{1/2}  -8.054056e+00 2     2

Neutron orbitals: 
Index State      Energy       Nodes Deg.
16 ^{1}s_{1/2}  -5.904000e+01 0     2
17 ^{1}p_{3/2}  -5.283627e+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.394392e+01 0     4
21 ^{2}s_{1/2}  -4.083663e+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.940208e+01 1     2
26 ^{1}g_{9/2}  -2.806149e+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.884692e+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.335126e+01 0     10
33 ^{2}f_{7/2}  -1.100552e+01 1     8
34 ^{1}i_{13/2} -9.576805e+00 0     14
35 ^{2}f_{5/2}  -8.991762e+00 1     6
36 ^{3}p_{3/2}  -8.290648e+00 2     4
37 ^{3}p_{1/2}  -7.515363e+00 2     2

Proton orbitals: 
Index State      Energy       Nodes Deg.
 0 ^{3}s_{1/2}  -8.057080e+00 2     2
 1 ^{1}h_{9/2}  -4.550104e+00 0     10
 2 ^{2}f_{7/2}  -1.365593e+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.639565e-01 0     16
 5 ^{3}d_{5/2}  -1.097961e+00 2     6

Proton RMS radius       :  5.460267e+00 fm.
Neutron RMS radius      :  5.740403e+00 fm.
Total energy per nucleon: -7.842550e+00 MeV.

1 tests performed.
All tests passed.