Classes | Public Attributes | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Types | List of all members
o2scl::nucleus_rmf Class Reference

Spherical closed-shell nuclei with a relativistic mean-field model in the Hartree approximation. More...

#include <nucleus_rmf.h>


struct  initial_guess
 Initial guess structure. More...
struct  odparms
 A convenient struct for the solution of the Dirac equations. More...
struct  shell
 A shell of nucleons for nucleus_rmf. More...

Public Member Functions

Basic operation
int run_nucleus (int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
 Computes the structure of a nucleus with the specified number of levels. More...
void set_verbose (int v)
 Set output level.
Lower-level interface
void init_run (int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
 Initialize a run. More...
void iterate (int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N, int &iconverged)
 Perform an iteration.
int post_converge (int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
 After convergence, make CM corrections, etc.

Public Attributes

initial_guess ig
 Parameters for initial guess. More...
bool generic_ode
 If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta.

Protected Member Functions

int load_nl3 (eos_had_rmf &r)
 Load the default model NL3 into the given eos_had_rmf object.
void init_meson_density ()
 Initialize the meson and photon fields, the densities, etc.
void energies (double xpro, double xnu, double e)
 Calculate the energy profile.
void center_mass_corr (double atot)
 Compute the center of mass correction.

Protected Attributes

double a_proton
 The parameter for the charge density of the proton (default is about 4.27073)
 The base EOS.
std::shared_ptr< table_units<> > profiles
 The radial profiles.
std::shared_ptr< table_units<> > chden_table
 The final charge densities.
std::vector< shell > * levp
 A pointer to the current vector of levels (either levels or unocc_levels)
int verbose
 Control output (default 1)
shell neutron_shells [n_internal_levels]
 The starting neutron levels.
shell proton_shells [n_internal_levels]
 The starting proton levels.
double step_size
 The grid step size (default 0.04)
double mnuc
 The nucleon mass (automatically set in init_fun())
ubvector energy
 Energy integrand.
bool init_called
 True if init() has been called.
ubvector ode_y
 ODE functions.
ubvector ode_dydx
 ODE derivatives.
ubvector ode_yerr
 ODE errors.
Density information (protected)
ubmatrix xrho
 The densities times radius squared.
ubvector xrhosp
 The proton scalar density times radius squared.
ubvector xrhos
 The scalar field RHS.
ubvector xrhov
 The vector field RHS.
ubvector xrhor
 The isubvector field RHS.
ubvector chden1
 Charge density.
ubvector chdenc
 Charge density.
ubvector arho
 Baryon density.
Gauss-Legendre integration points and weights
double x12 [6]
double w12 [6]
double x100 [50]
double w100 [50]

Static Protected Attributes

static const int n_internal_levels =29
 The total number of shells stored internally.
static const int grid_size =300
 The grid size.

Private Types

typedef boost::numeric::ublas::vector< double > ubvector
typedef boost::numeric::ublas::matrix< double > ubmatrix


int nlevels
 The number of levels.
std::vector< shelllevels
 The levels (protons first, then neutrons) More...
int nuolevels
 The number of unoccupied levels (equal to unocc_Z + unocc_N)
std::vector< shellunocc_levels
 The unoccupied levels (protons first, then neutrons) More...
double stens
 Surface tension (in $ \mathrm{fm}^{-3} $ ) More...
double rnrp
 Skin thickness (in fm) More...
double rnrms
 Neutron RMS radius (in fm) More...
double rprms
 Proton RMS radius (in fm) More...
double etot
 Total energy (in MeV) More...
double r_charge
 Charge radius (in fm) More...
double r_charge_cm
 Charge radius corrected by the center of mass (in fm) More...
std::shared_ptr< table_units<> > get_profiles ()
 Get the radial profiles. More...
std::shared_ptr< table_units<> > get_chden ()
 The final charge densities.

Equation of state

eos_had_rmf def_rmf
 The default equation of state (default NL3) More...
thermo hb
 thermo object for the EOS More...
fermion n
 The neutron. More...
fermion p
 The proton. More...
int set_eos (eos_had_rmf &r)
 Set the base EOS to be used. More...

Numeric configuration

bool err_nonconv
 If true, call the error handler if the routine does not converge or reach the desired tolerance (default true)
int itmax
 Maximum number of total iterations (default 70)
int meson_itmax
 Maximum number of iterations for solving the meson field equations (default 10000)
int dirac_itmax
 Maximum number of iterations for solving the Dirac equations (default 100)
double dirac_tol
 Tolerance for Dirac equations (default $ 5 \times 10^{-3} $ ).
double dirac_tol2
 Second tolerance for Dirac equations (default $ 5 \times 10^{-4} $ ).
double meson_tol
 Tolerance for meson field equations (default $ 10^{-6} $ ).
void set_step (ode_step< ubvector, ubvector, ubvector, ode_funct11 > &step)
 Set the stepper for the Dirac differential equation.

The meson and photon fields and field equations (protected)

ubmatrix field0
 Values of the fields from the last iteration.
ubmatrix fields
 The values of the fields.
ubmatrix gin
 The Green's functions inside.
ubmatrix gout
 The Green's functions outside.
int surf_index
 The grid index corresponding to the nuclear surface (computed by init_run())
double sigma_rhs (double sig, double ome, double rho)
 Scalar density RHS.
double omega_rhs (double sig, double ome, double rho)
 Vector density RHS.
double rho_rhs (double sig, double ome, double rho)
 Isubvector density RHS.
void meson_init ()
 Calculate meson and photon Green's functions gin and gout.
void meson_iter (double ic)
 Calculate meson and photon fields. More...
void meson_solve ()
 Solve for the meson profiles.

Calculating the form factor, etc. (protected)

interp_vec< ubvector > * gi
 Interpolation object.
void pfold (double x, double &xrhof)
 Fold in proton form factor.
double xpform (double x, double xp, double a)
 Function representing proton form factor.
void gauss (double xmin, double xmax, double x, double &xi)
 Perform integrations for form factor.
double xrhop (double x1)

Solving the Dirac equations (protected)

ode_rkck_gsl< ubvector, ubvector, ubvector, ode_funct11def_step
 The default stepper.
ode_step< ubvector, ubvector, ubvector, ode_funct11 > * ostep
 The ODE stepper.
void dirac (int ilevel)
 Solve the Dirac equations. More...
void dirac_step (double &x, double h, double eigen, double kappa, ubvector &varr)
 Take a step in the Dirac equations.
int odefun (double x, size_t nv, const ubvector &y, ubvector &dydx, odparms &op)
 The form of the Dirac equations for the ODE stepper.
void field (double x, double &s, double &v, ubvector &varr)
 Compute the fields for the Dirac equations.
double dirac_rk4 (double x, double g1, double f1, double &funt, double eigen, double kappa, ubvector &varr)
 Integrate the Dirac equations using a simple inline 4th order Runge-Kutta.

Detailed Description

This code is very experimental.

This class is based on a code developed by C.J. Horowitz and B.D. Serot, and used in Horowitz81 which was then adapted by P.J. Ellis and used in Heide94 and Prakash94. Ellis and A.W. Steiner adapted it for the parameterization in in eos_had_rmf for Steiner05b, and then converted to C++ by Steiner afterwards.

The standard usage is something like:

nucleus_rmf rn;
cout << rn.rnrp << endl;

which computes the structure of $ ^{208}\mathrm{Pb} $ and outputs the neutron skin thickness using the model 'NL4'.

Potential exceptions are

The initial level pattern is

1 S 1/2
// 2 nucleons
1 P 3/2
1 P 1/2
// 8 nucleus
1 D 5/2
1 D 3/2
2 S 1/2
// 20 nucleons
1 F 7/2
// 28 nucleons
1 F 5/2
2 P 3/2
2 P 1/2
// 40 nucleons
1 G 9/2
// 50 nucleus
1 G 7/2
2 D 5/2
1 H 11/2
2 D 3/2
3 S 1/2
// 82 nucleons
1 H 9/2
2 F 7/2
1 I 13/2
2 F 5/2
3 P 3/2
3 P 1/2
// 126 nucleons
2 G 9/2
1 I 11/2
1 J 15/2
3 D 5/2
4 S 1/2
2 G 7/2
3 D 3/2
// 184 nucleons

Below, $ \alpha $ is a generic index for the isospin, the radial quantum number $ n $ and the angular quantum numbers $ \kappa $ and $ m $. The meson fields are $ \sigma(r), \omega(r) $ and $ \rho(r) $. The baryon density is $ n(r) $, the neutron and proton densities are $ n_n(r) $ and $ n_p(r) $, and the baryon scalar density is $ n_s(r) $. The nucleon field equations are

\begin{eqnarray*} F^{\prime}_{\alpha}(r)- \frac{\kappa}{r} F_{\alpha}(r) + \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r) - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right) e A(r) - M + g_{\sigma} \sigma(r) \right] G_{\alpha}(r) &=& 0 \\ G^{\prime}_{\alpha}(r)+ \frac{\kappa}{r} G_{\alpha}(r) - \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r) - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right) e A(r) + M - g_{\sigma} \sigma(r) \right] F_{\alpha}(r) &=& 0 \end{eqnarray*}

where the isospin number, $ t_{\alpha} $ is $ 1/2 $ for protons and $ -1/2 $ for neutrons. The meson field equations are

\begin{eqnarray*} \sigma^{\prime \prime}(r) + \frac{2}{r} \sigma^{\prime}(r) - m_{\sigma}^2 \sigma &=& - g_{\sigma} n_s(r) + b M g_{\sigma}^3 \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \sigma} \\ \omega^{\prime \prime}(r) + \frac{2}{r} \omega^{\prime}(r) - m_{\omega}^2 \omega &=& - g_{\omega} n(r) + \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \omega} \\ \rho^{\prime \prime}(r) + \frac{2}{r} \rho^{\prime}(r) - m_{\rho}^2 \rho &=& - \frac{g_{\rho}}{2} \left[n_n(r)-n_p(r)\right] + 2 g_{\rho}^2 \rho f + \frac{\xi}{6} g_{\rho}^4 \rho^3 \end{eqnarray*}

and the Coulomb field equation is

\[ A^{\prime \prime}(r) + \frac{2}{r} A^{\prime}(r) = - e n_p(r) \]

The meson field equations plus a pair of Dirac-like nucleon field equations for each index $ \alpha $ must be solved to compute the structure of a given nucleus.

The densities (scalar, baryon, isospin, and charge) are

\begin{eqnarray*} n_s(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \\ n(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2+F(r)^2 \right] \right\} \\ n_i(r) &=& \sum_{\alpha} \left\{ t_{\alpha} \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \\ n_c(r) &=& \sum_{\alpha} \left\{ \left[t_{\alpha}+\frac{1}{2}\right] \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \end{eqnarray*}

The total energy is

\[ E = \sum_{\alpha} \varepsilon_{\alpha} (2 j_{\alpha}+1) - \frac{1}{2} \int d^{3} r \left[- g_{\sigma} \sigma(r) \rho_s(r) + g_{\omega} \omega(r) \rho(r) + \frac{1}{2} g_{\rho} \rho(r) + e A(r) n_p(r) \right] \]

The charge density is the proton density folded with the charge density of the proton, i.e.

\[ \rho_{\mathrm{ch}}(r) = \int d^{3} r^{\prime} \rho_{\mathrm{prot}}(r-r^{\prime}) \rho_p(r) \]

where the proton charge density is assumed to be of the form

\[ \rho_{\mathrm{prot}}(r) = \frac{\mu^3}{8 \pi} \exp \left( - \mu |r|\right) \]

and the parameter $ \mu = (0.71)^{1/2}~\mathrm{GeV} $ (see Eq. 20b in Horowitz81). The default value of a_proton is the value of $ \mu $ converted into $ \mathrm{fm}^{-1} $.

Generally, the first array index associated with a function is not the value at $ r=0 $, but at $ r=\Delta r $ (stored in step_size) and so the $ r=0 $ part of the algorithm is handled separately.


Better documentation

Convert energies() to use EOS and possibly replace sigma_rhs() and related functions by the associated field equation method of eos_had_rmf.

Document hw=3.923+23.265/cbrt(atot);
I believe currently the center of mass correction for the binding energy per nucleon is not done and has to be added in after the fact
Idea for Future:

Sort energy levels at the end by energy

Improve the numerical methods

Make the neutron and proton orbitals more configurable

Generalize to $ m_n \neq m_p $ .

Allow more freedom in the integrations

Consider converting everything to inverse fermis.

Convert to zero-indexed arrays (mostly done)

Warn when the level ordering is wrong, and unoccupied levels are lower energy than occupied levels

Connect with o2scl::nucmass ?

Definition at line 229 of file nucleus_rmf.h.

Member Function Documentation

◆ dirac()

void o2scl::nucleus_rmf::dirac ( int  ilevel)

Solves the Dirac equation in from 12 fm to the match point and then out from .04 fm and adjusts eigenvalue with

\[ \Delta \varepsilon = -g(r=\mathrm{match\_point}) \times (f^{+}-f^{-}) \]

◆ get_profiles()

std::shared_ptr<table_units<> > o2scl::nucleus_rmf::get_profiles ( )

The profiles are calculated each iteration by iterate().

Definition at line 324 of file nucleus_rmf.h.

◆ init_run()

void o2scl::nucleus_rmf::init_run ( int  nucleus_Z,
int  nucleus_N,
int  unocc_Z,
int  unocc_N 

Note that rmf must be set before run_nucleus() is called.

◆ meson_iter()

void o2scl::nucleus_rmf::meson_iter ( double  ic)

The meson and photon field equations are of the Sturm-Liouville form, e.g.

\[ \left[r^2 \sigma^{\prime}(r) \right]^{\prime} - r^2 m_{\sigma}^2 \sigma(r) = f(r) \]

where $ \sigma(0) = \sigma_0 $ and $ \sigma(+\infty) = 0 $. A solution of the homogeneous equation with $ f(r) =0 $ is $ \sigma(r) = \sigma_0 \sinh( m_{\sigma} r)/ (m_{\sigma} r) $. The associated Green's function is

\[ D(r,r^{\prime},m_{\sigma}) = \frac{-1}{m_{\sigma} r r^{\prime}} \sinh (m_{\sigma} r_{<}) \exp (-m_{\sigma} r_{>}) \]

where $ r_{<} $ is the smaller of $ r $ and $ r^{\prime} $ and $ r_{>} $ is the larger. One can then write the solution of the meson field equation given the density

\[ \sigma(r) = \int_0^{\infty} r^{\prime 2}~d r^{\prime} \left[ - g_{\sigma} n_{s}(r) \right] D\left(r,r^{\prime},m_{\sigma}\right) \]

Since radii are positive, $ \sinh (m r) \approx e^{m r}/2 $ and

\[ D(r,r^{\prime},m_{\sigma}) \approx \left[ \frac{-1}{2 m_{\sigma} r_{<}} \exp (m_{\sigma} r_{<}) \right] \left[ \frac{1}{r_{>}} \exp (-m_{\sigma} r_{>}) \right] \]

The two terms in the Green's function are separated into

\[ \mathrm{gin}(r) = \frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \]


\[ \mathrm{gout}(r) = \frac{e^{-m_{\sigma} r}}{r} \, . \]

These functions are computed in meson_init() . Then the field is given by

\[ \sigma(r)= \left(\frac{e^{-m_{\sigma} r}}{r}\right) \int_0^{r} r^{\prime 2} g_{\sigma} n_{s} \left(\frac{e^{m_{\sigma} r^{\prime}}}{2 m_{\sigma} r^{\prime}} \right)~d r^{\prime} + \left(\frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \right) \int_r^{\infty} r^{\prime 2} g_{\sigma} n_{s} \left(\frac{e^{-m_{\sigma} r^{\prime}}} {r^{\prime}}\right)~d r^{\prime} \]

where the first integral is stored in xi2 and the second is in xi1 in the function meson_iter() . The part of xi2 at $ r^{\prime}=0 $ is stored in xi20.

◆ run_nucleus()

int o2scl::nucleus_rmf::run_nucleus ( int  nucleus_Z,
int  nucleus_N,
int  unocc_Z,
int  unocc_N 

Note that rmf must be set before run_nucleus() is called.

This calls init_run(), and then iterate() until iconverged is 1, and then post_converge().

◆ set_eos()

int o2scl::nucleus_rmf::set_eos ( eos_had_rmf r)

The equation of state must be set before run_nucleus() or init_fun() are called, including the value of eos_had_rmf::mnuc.

Definition at line 410 of file nucleus_rmf.h.

Member Data Documentation

◆ def_rmf

eos_had_rmf o2scl::nucleus_rmf::def_rmf

This is set in the constructor to be the default model, NL3, using the function load_nl3().

Definition at line 403 of file nucleus_rmf.h.

◆ etot

double o2scl::nucleus_rmf::etot

Computed every iteration in iterate() or automatically in run_nucleus()

Definition at line 380 of file nucleus_rmf.h.

◆ hb

thermo o2scl::nucleus_rmf::hb

This is just used as temporary storage.

Definition at line 419 of file nucleus_rmf.h.

◆ ig

initial_guess o2scl::nucleus_rmf::ig

Default is {310,240,-6,25.9,6.85,0.6}

Definition at line 504 of file nucleus_rmf.h.

◆ levels

std::vector<shell> o2scl::nucleus_rmf::levels

An array of size nlevels

Definition at line 337 of file nucleus_rmf.h.

◆ n

fermion o2scl::nucleus_rmf::n

The mass of the neutron is ignored and set by init_run() to be eos_had_rmf::mnuc from rmf.

Definition at line 426 of file nucleus_rmf.h.

◆ p

fermion o2scl::nucleus_rmf::p

The mass of the proton is ignored and set by init_run() to be eos_had_rmf::mnuc from rmf.

Definition at line 433 of file nucleus_rmf.h.

◆ r_charge

double o2scl::nucleus_rmf::r_charge

Computed in post_converge() or automatically in run_nucleus()

Definition at line 386 of file nucleus_rmf.h.

◆ r_charge_cm

double o2scl::nucleus_rmf::r_charge_cm

Computed in post_converge() or automatically in run_nucleus()

Definition at line 392 of file nucleus_rmf.h.

◆ rnrms

double o2scl::nucleus_rmf::rnrms

Computed every iteration in iterate() or automatically in run_nucleus()

Definition at line 366 of file nucleus_rmf.h.

◆ rnrp

double o2scl::nucleus_rmf::rnrp

Computed every iteration in iterate() or automatically in run_nucleus()

Definition at line 359 of file nucleus_rmf.h.

◆ rprms

double o2scl::nucleus_rmf::rprms

Computed every iteration in iterate() or automatically in run_nucleus()

Definition at line 373 of file nucleus_rmf.h.

◆ stens

double o2scl::nucleus_rmf::stens

Computed in post_converge() or automatically in run_nucleus()

Definition at line 352 of file nucleus_rmf.h.

◆ unocc_levels

std::vector<shell> o2scl::nucleus_rmf::unocc_levels

An array of size nuolevels

Definition at line 346 of file nucleus_rmf.h.

The documentation for this class was generated from the following file:

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).