nucleus_rmf.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file nucleus_rmf.h
24  \brief File defining \ref o2scl::nucleus_rmf
25 */
26 #ifndef RMF_NUCLEUS_H
27 #define RMF_NUCLEUS_H
28 
29 #include <iostream>
30 #include <string>
31 #include <vector>
32 #include <o2scl/interp.h>
33 #include <o2scl/constants.h>
34 #include <o2scl/part.h>
35 #include <o2scl/eos_had_rmf.h>
36 #include <o2scl/table_units.h>
37 #include <o2scl/ode_rkck_gsl.h>
38 #include <o2scl/ode_funct.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Spherical closed-shell nuclei with a relativistic
45  mean-field model in the Hartree approximation
46 
47  This code is very experimental.
48 
49  This class is based on a code developed by C.J. Horowitz and
50  B.D. Serot, and used in \ref Horowitz81 which was then adapted
51  by P.J. Ellis and used in \ref Heide94 and \ref Prakash94. Ellis
52  and A.W. Steiner adapted it for the parameterization in in \ref
53  eos_had_rmf for \ref Steiner05b, and then converted to C++ by
54  Steiner afterwards.
55 
56  The standard usage is something like:
57  \code
58  nucleus_rmf rn;
59  o2scl_hdf::rmf_load(rn.rmf,"NL4");
60  rn.run_nucleus(82,208,0,0);
61  cout << rn.rnrp << endl;
62  \endcode
63  which computes the structure of \f$ ^{208}\mathrm{Pb} \f$ and
64  outputs the neutron skin thickness using the model \c 'NL4'.
65 
66  Potential exceptions are
67  - Failed to converge
68  - Failed to solve meson field equations
69  - Energy not finite (usually a problem in the equation of
70  state)
71  - Energy not finite in final calculation
72  - Function \ref iterate() called before \ref init_run()
73  - Not a closed-shell nucleus
74 
75  The initial level pattern is
76  \verbatim
77  1 S 1/2
78  // 2 nucleons
79  1 P 3/2
80  1 P 1/2
81  // 8 nucleus
82  1 D 5/2
83  1 D 3/2
84  2 S 1/2
85  // 20 nucleons
86  1 F 7/2
87  // 28 nucleons
88  1 F 5/2
89  2 P 3/2
90  2 P 1/2
91  // 40 nucleons
92  1 G 9/2
93  // 50 nucleus
94  1 G 7/2
95  2 D 5/2
96  1 H 11/2
97  2 D 3/2
98  3 S 1/2
99  // 82 nucleons
100  1 H 9/2
101  2 F 7/2
102  1 I 13/2
103  2 F 5/2
104  3 P 3/2
105  3 P 1/2
106  // 126 nucleons
107  2 G 9/2
108  1 I 11/2
109  1 J 15/2
110  3 D 5/2
111  4 S 1/2
112  2 G 7/2
113  3 D 3/2
114  // 184 nucleons
115  \endverbatim
116 
117  Below, \f$ \alpha \f$ is a generic index for the isospin, the
118  radial quantum number \f$ n \f$ and the angular quantum numbers
119  \f$ \kappa \f$ and \f$ m \f$. The meson fields are \f$
120  \sigma(r), \omega(r) \f$ and \f$ \rho(r) \f$. The baryon density
121  is \f$ n(r) \f$, the neutron and proton densities are \f$ n_n(r)
122  \f$ and \f$ n_p(r) \f$, and the baryon scalar density is \f$
123  n_s(r) \f$.
124  The nucleon field equations are
125  \f{eqnarray*}
126  F^{\prime}_{\alpha}(r)- \frac{\kappa}{r} F_{\alpha}(r)
127  + \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r)
128  - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right)
129  e A(r) - M + g_{\sigma} \sigma(r) \right] G_{\alpha}(r) &=& 0 \\
130  G^{\prime}_{\alpha}(r)+ \frac{\kappa}{r} G_{\alpha}(r)
131  - \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r)
132  - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2}
133  \right) e A(r) + M - g_{\sigma} \sigma(r) \right] F_{\alpha}(r) &=& 0
134  \f}
135  where the isospin number, \f$ t_{\alpha} \f$ is \f$ 1/2 \f$ for
136  protons and \f$ -1/2 \f$ for neutrons.
137  The meson field equations are
138  \f{eqnarray*}
139  \sigma^{\prime \prime}(r) + \frac{2}{r} \sigma^{\prime}(r)
140  - m_{\sigma}^2 \sigma &=& - g_{\sigma} n_s(r) + b M g_{\sigma}^3
141  \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2
142  \frac{\partial f}{\partial \sigma} \\
143  \omega^{\prime \prime}(r) + \frac{2}{r} \omega^{\prime}(r)
144  - m_{\omega}^2 \omega &=& - g_{\omega} n(r) +
145  \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2
146  \frac{\partial f}{\partial \omega} \\
147  \rho^{\prime \prime}(r) + \frac{2}{r} \rho^{\prime}(r)
148  - m_{\rho}^2 \rho &=& - \frac{g_{\rho}}{2}
149  \left[n_n(r)-n_p(r)\right] + 2 g_{\rho}^2 \rho f + \frac{\xi}{6}
150  g_{\rho}^4 \rho^3
151  \f}
152  and the Coulomb field equation is
153  \f[
154  A^{\prime \prime}(r) + \frac{2}{r} A^{\prime}(r) = - e n_p(r)
155  \f]
156  The meson field equations plus a pair of Dirac-like nucleon
157  field equations for each index \f$ \alpha \f$ must be solved to
158  compute the structure of a given nucleus.
159 
160  The densities (scalar, baryon, isospin, and charge) are
161  \f{eqnarray*}
162  n_s(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2-F(r)^2
163  \right] \right\} \\
164  n(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2+F(r)^2
165  \right] \right\} \\
166  n_i(r) &=& \sum_{\alpha} \left\{ t_{\alpha} \int d^3 r
167  \left[ G(r)^2-F(r)^2 \right] \right\} \\
168  n_c(r) &=& \sum_{\alpha} \left\{ \left[t_{\alpha}+\frac{1}{2}\right]
169  \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\}
170  \f}
171 
172  The total energy is
173  \f[
174  E = \sum_{\alpha} \varepsilon_{\alpha} (2 j_{\alpha}+1)
175  - \frac{1}{2} \int d^{3} r
176  \left[- g_{\sigma} \sigma(r) \rho_s(r) + g_{\omega} \omega(r)
177  \rho(r) + \frac{1}{2} g_{\rho} \rho(r) + e A(r) n_p(r) \right]
178  \f]
179 
180  The charge density is the proton density folded with the
181  charge density of the proton, i.e.
182  \f[
183  \rho_{\mathrm{ch}}(r) = \int d^{3} r^{\prime}
184  \rho_{\mathrm{prot}}(r-r^{\prime}) \rho_p(r)
185  \f]
186  where the proton charge density is assumed to be of the form
187  \f[
188  \rho_{\mathrm{prot}}(r) = \frac{\mu^3}{8 \pi} \exp \left(
189  - \mu |r|\right)
190  \f]
191  and the parameter \f$ \mu = (0.71)^{1/2}~\mathrm{GeV} \f$ (see
192  Eq. 20b in \ref Horowitz81). The default value of \ref a_proton
193  is the value of \f$ \mu \f$ converted into \f$ \mathrm{fm}^{-1}
194  \f$.
195 
196  Generally, the first array index associated with a function
197  is not the value at \f$ r=0 \f$, but at \f$ r=\Delta r \f$
198  (stored in \ref step_size) and so the \f$ r=0 \f$ part of the
199  algorithm is handled separately.
200 
201  \todo Better documentation
202  \todo Convert energies() to use EOS and possibly
203  replace sigma_rhs() and related functions by the associated
204  field equation method of eos_had_rmf.
205 
206  \todo Document hw=3.923+23.265/cbrt(atot);
207 
208  \comment
209  the hw=blah blah term for the CM correction is discussed
210  a bit in Negele, PRC 1 (1970) 1260, esp. see Eq. 2.30 but
211  the numerical coefficients are different here.
212  \endcomment
213 
214  \todo I believe currently the center of mass correction
215  for the binding energy per nucleon is not done and has
216  to be added in after the fact
217 
218  \future Sort energy levels at the end by energy
219  \future Improve the numerical methods
220  \future Make the neutron and proton orbitals more configurable
221  \future Generalize to \f$ m_n \neq m_p \f$ .
222  \future Allow more freedom in the integrations
223  \future Consider converting everything to inverse fermis.
224  \future Convert to zero-indexed arrays (mostly done)
225  \future Warn when the level ordering is wrong, and unoccupied levels
226  are lower energy than occupied levels
227  \future Connect with \ref o2scl::nucmass ?
228  */
229  class nucleus_rmf {
230 
233 
234 #ifndef DOXYGEN_INTERNAL
235 
236  protected:
237 
238  /// A convenient struct for the solution of the Dirac equations
239  typedef struct {
240  /// Eigenvalue
241  double eigen;
242  /// Quantum number \f$ \kappa \f$ .
243  double kappa;
244  /// The meson fields
245  ubmatrix *fields;
246  /// Desc
247  ubvector *varr;
248  } odparms;
249 
250  /// The total number of shells stored internally
251  static const int n_internal_levels=29;
252 
253 #endif
254 
255  public:
256 
257  nucleus_rmf();
258 
259  ~nucleus_rmf();
260 
261  /** \name Basic operation
262  */
263  //@{
264  /// A shell of nucleons for \ref nucleus_rmf
265  typedef struct {
266  /// Degeneracy \f$ 2 j+1 \f$ .
267  int twojp1;
268  /// \f$ \kappa \f$
269  int kappa;
270  /// Energy eigenvalue
271  double energy;
272  /// Isospin ( \f$ +1/2 \f$ or \f$ -1/2 \f$ .
273  double isospin;
274  /// Angular momentum-spin state \f$ ^{2s+1} \ell_{j} \f$
275  std::string state;
276  /// Matching radius (in fm)
277  double match_point;
278  /// Desc.
279  double eigen;
280  /// Desc.
281  double eigenc;
282  /// Number of nodes in the wave function
283  int nodes;
284  } shell;
285 
286  /** \brief Computes the structure of a nucleus with the specified
287  number of levels
288 
289  Note that \ref rmf must be set before run_nucleus() is called.
290 
291  This calls init_run(), and then iterate() until \c iconverged is
292  1, and then post_converge().
293  */
294  int run_nucleus(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
295 
296  /// Set output level
297  void set_verbose(int v) { verbose=v; };
298  //@}
299 
300  /** \name Lower-level interface
301  */
302  //@{
303  /** \brief Initialize a run
304 
305  Note that \ref rmf must be set before run_nucleus() is called.
306  */
307  void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
308 
309  /// Perform an iteration
310  void iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N,
311  int &iconverged);
312 
313  /// After convergence, make CM corrections, etc.
314  int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
315 
316  //@}
317 
318  /// \name Results
319  //@{
320  /** \brief Get the radial profiles
321 
322  The profiles are calculated each iteration by iterate().
323  */
324  std::shared_ptr<table_units<> > get_profiles() { return profiles; };
325 
326  /** \brief The final charge densities
327  */
328  std::shared_ptr<table_units<> > get_chden() { return chden_table; };
329 
330  /// The number of levels
331  int nlevels;
332 
333  /** \brief The levels (protons first, then neutrons)
334 
335  An array of size \ref nlevels
336  */
337  std::vector<shell> levels;
338 
339  /// The number of unoccupied levels (equal to \c unocc_Z + \c unocc_N)
341 
342  /** \brief The unoccupied levels (protons first, then neutrons)
343 
344  An array of size \ref nuolevels
345  */
346  std::vector<shell> unocc_levels;
347 
348  /** \brief Surface tension (in \f$ \mathrm{fm}^{-3} \f$ )
349 
350  Computed in post_converge() or automatically in run_nucleus()
351  */
352  double stens;
353 
354  /** \brief Skin thickness (in fm)
355 
356  Computed every iteration in iterate()
357  or automatically in run_nucleus()
358  */
359  double rnrp;
360 
361  /** \brief Neutron RMS radius (in fm)
362 
363  Computed every iteration in iterate() or automatically in
364  run_nucleus()
365  */
366  double rnrms;
367 
368  /** \brief Proton RMS radius (in fm)
369 
370  Computed every iteration in iterate() or automatically in
371  run_nucleus()
372  */
373  double rprms;
374 
375  /** \brief Total energy (in MeV)
376 
377  Computed every iteration in iterate() or automatically in
378  run_nucleus()
379  */
380  double etot;
381 
382  /** \brief Charge radius (in fm)
383 
384  Computed in post_converge() or automatically in run_nucleus()
385  */
386  double r_charge;
387 
388  /** \brief Charge radius corrected by the center of mass (in fm)
389 
390  Computed in post_converge() or automatically in run_nucleus()
391  */
392  double r_charge_cm;
393  //@}
394 
395  /** \name Equation of state
396  */
397  //@{
398  /** \brief The default equation of state (default NL3)
399 
400  This is set in the constructor to be the default
401  model, NL3, using the function \ref load_nl3().
402  */
404 
405  /** \brief Set the base EOS to be used
406 
407  The equation of state must be set before run_nucleus() or
408  init_fun() are called, including the value of eos_had_rmf::mnuc.
409  */
411  rmf=&r;
412  return 0;
413  }
414 
415  /** \brief \ref thermo object for the EOS
416 
417  This is just used as temporary storage.
418  */
420 
421  /** \brief The neutron
422 
423  The mass of the neutron is ignored and set by init_run()
424  to be eos_had_rmf::mnuc from \ref rmf.
425  */
427 
428  /** \brief The proton
429 
430  The mass of the proton is ignored and set by init_run()
431  to be eos_had_rmf::mnuc from \ref rmf.
432  */
434  //@}
435 
436  /** \name Numeric configuration
437  */
438  //@{
439  /** \brief If true, call the error handler if the routine does not
440  converge or reach the desired tolerance (default true)
441  */
443 
444  /// Set the stepper for the Dirac differential equation
445  void set_step(ode_step<ubvector,ubvector,ubvector,
446  ode_funct11> &step) {
447  ostep=&step;
448  return;
449  }
450 
451  /// Maximum number of total iterations (default 70)
452  int itmax;
453 
454  /** \brief Maximum number of iterations for solving the meson
455  field equations (default 10000)
456  */
458 
459  /** \brief Maximum number of iterations for solving the Dirac
460  equations (default 100)
461  */
463 
464  /// Tolerance for Dirac equations (default \f$ 5 \times 10^{-3} \f$ ).
465  double dirac_tol;
466 
467  /** \brief Second tolerance for Dirac equations
468  (default \f$ 5 \times 10^{-4} \f$ ).
469  */
470  double dirac_tol2;
471 
472  /// Tolerance for meson field equations (default \f$ 10^{-6} \f$ ).
473  double meson_tol;
474  //@}
475 
476  /** \brief Initial guess structure
477 
478  The initial guess for the meson field profiles is
479  a set of Fermi-Dirac functions, i.e.
480  \f[
481  \sigma(r)=\mathrm{sigma0}/
482  [1+\exp((r-\mathrm{fermi\_radius})/\mathrm{fermi\_width})]
483  \f]
484  */
485  typedef struct {
486  /// Scalar field at r=0
487  double sigma0;
488  /// Vector field at r=0
489  double omega0;
490  /// Isubvector field at r=0
491  double rho0;
492  /// Coulomb field at r=0
493  double A0;
494  /// The radius for which the fields are half their central value
495  double fermi_radius;
496  /// The "width" of the Fermi-Dirac function
497  double fermi_width;
498  } initial_guess;
499 
500  /** \brief Parameters for initial guess
501 
502  Default is {310,240,-6,25.9,6.85,0.6}
503  */
505 
506  /** \brief If true, use the generic ODE solver instead of the
507  internal 4th order Runge-Kutta
508  */
510 
511 #ifndef DOXYGEN_INTERNAL
512 
513  protected:
514 
515  /** \brief The parameter for the charge density of the proton
516  (default is about 4.27073)
517  */
518  double a_proton;
519 
520  /// Load the default model NL3 into the given \ref eos_had_rmf object
521  int load_nl3(eos_had_rmf &r);
522 
523  /// The base EOS
525 
526  /** \brief The radial profiles
527  */
528  std::shared_ptr<table_units<> > profiles;
529 
530  /** \brief The final charge densities
531  */
532  std::shared_ptr<table_units<> > chden_table;
533 
534  /** \brief A pointer to the current vector of levels
535  (either \ref levels or \ref unocc_levels)
536  */
537  std::vector<shell> *levp;
538 
539  /** \brief Control output (default 1)
540  */
541  int verbose;
542 
543  /// The starting neutron levels
545 
546  /// The starting proton levels
548 
549  /// The grid size
550  static const int grid_size=300;
551 
552  /// The grid step size (default 0.04)
553  double step_size;
554 
555  /// The nucleon mass (automatically set in init_fun())
556  double mnuc;
557 
558  /** \name The meson and photon fields and field equations (protected)
559  */
560  //@{
561  /// Values of the fields from the last iteration
562  ubmatrix field0;
563 
564  /// The values of the fields
565  ubmatrix fields;
566 
567  /// The Green's functions inside
568  ubmatrix gin;
569 
570  /// The Green's functions outside
571  ubmatrix gout;
572 
573  /// Scalar density RHS
574  double sigma_rhs(double sig, double ome, double rho);
575 
576  /// Vector density RHS
577  double omega_rhs(double sig, double ome, double rho);
578 
579  /// Isubvector density RHS
580  double rho_rhs(double sig, double ome, double rho);
581 
582  /** \brief Calculate meson and photon
583  Green's functions \ref gin and \ref gout
584  */
585  void meson_init();
586 
587  /** \brief Calculate meson and photon fields
588 
589  The meson and photon field equations are of the
590  Sturm-Liouville form, e.g.
591  \f[
592  \left[r^2 \sigma^{\prime}(r) \right]^{\prime} -
593  r^2 m_{\sigma}^2 \sigma(r) = f(r)
594  \f]
595  where \f$ \sigma(0) = \sigma_0 \f$ and \f$ \sigma(+\infty) = 0
596  \f$. A solution of the homogeneous equation with \f$ f(r) =0
597  \f$ is \f$ \sigma(r) = \sigma_0 \sinh( m_{\sigma} r)/
598  (m_{\sigma} r) \f$. The associated Green's function is
599  \f[
600  D(r,r^{\prime},m_{\sigma}) = \frac{-1}{m_{\sigma} r r^{\prime}}
601  \sinh (m_{\sigma} r_{<}) \exp (-m_{\sigma} r_{>})
602  \f]
603  where \f$ r_{<} \f$ is the smaller of \f$ r \f$ and
604  \f$ r^{\prime} \f$ and \f$ r_{>} \f$ is the larger.
605  One can then write the solution of the meson field
606  equation given the density
607  \f[
608  \sigma(r) = \int_0^{\infty} r^{\prime 2}~d r^{\prime}
609  \left[ - g_{\sigma} n_{s}(r) \right]
610  D\left(r,r^{\prime},m_{\sigma}\right)
611  \f]
612 
613  Since radii are positive, \f$ \sinh (m r) \approx
614  e^{m r}/2 \f$ and
615  \f[
616  D(r,r^{\prime},m_{\sigma}) \approx \left[
617  \frac{-1}{2 m_{\sigma} r_{<}}
618  \exp (m_{\sigma} r_{<})
619  \right] \left[
620  \frac{1}{r_{>}}
621  \exp (-m_{\sigma} r_{>})
622  \right]
623  \f]
624  The two terms in the Green's function are separated into
625  \f[
626  \mathrm{gin}(r) = \frac{e^{m_{\sigma} r}}{2 m_{\sigma} r}
627  \f]
628  and
629  \f[
630  \mathrm{gout}(r) = \frac{e^{-m_{\sigma} r}}{r} \, .
631  \f]
632  These functions are computed in \ref meson_init() . Then the
633  field is given by
634  \f[
635  \sigma(r)= \left(\frac{e^{-m_{\sigma} r}}{r}\right)
636  \int_0^{r} r^{\prime 2} g_{\sigma} n_{s}
637  \left(\frac{e^{m_{\sigma} r^{\prime}}}{2 m_{\sigma}
638  r^{\prime}} \right)~d r^{\prime} +
639  \left(\frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \right)
640  \int_r^{\infty} r^{\prime 2} g_{\sigma} n_{s}
641  \left(\frac{e^{-m_{\sigma} r^{\prime}}}
642  {r^{\prime}}\right)~d r^{\prime}
643  \f]
644  where the first integral is stored in <tt>xi2</tt> and
645  the second is in <tt>xi1</tt> in the function \ref meson_iter() .
646  The part of <tt>xi2</tt> at \f$ r^{\prime}=0 \f$ is
647  stored in <tt>xi20</tt>.
648  */
649  void meson_iter(double ic);
650 
651  /** \brief Solve for the meson profiles
652  */
653  void meson_solve();
654 
655  /** \brief The grid index corresponding to the nuclear surface
656  (computed by init_run())
657  */
659  //@}
660 
661  /** \name Density information (protected)
662  */
663  //@{
664  /// The densities times radius squared
665  ubmatrix xrho;
666 
667  /// The proton scalar density times radius squared
668  ubvector xrhosp;
669 
670  /// The scalar field RHS
671  ubvector xrhos;
672 
673  /// The vector field RHS
674  ubvector xrhov;
675 
676  /// The isubvector field RHS
677  ubvector xrhor;
678 
679  /// Charge density
680  ubvector chden1;
681 
682  /// Charge density
683  ubvector chdenc;
684 
685  /// Baryon density
686  ubvector arho;
687  //@}
688 
689  /// Energy integrand
690  ubvector energy;
691 
692  /// Initialize the meson and photon fields, the densities, etc.
693  void init_meson_density();
694 
695  /// Calculate the energy profile
696  void energies(double xpro, double xnu, double e);
697 
698  /// Compute the center of mass correction
699  void center_mass_corr(double atot);
700 
701  /** \name Calculating the form factor, etc. (protected)
702  */
703  //@{
704  /// Fold in proton form factor
705  void pfold(double x, double &xrhof);
706 
707  /// Function representing proton form factor
708  double xpform(double x, double xp, double a);
709 
710  /// Perform integrations for form factor
711  void gauss(double xmin, double xmax, double x, double &xi);
712 
713  /// Desc.
714  double xrhop(double x1);
715 
716  /// Interpolation object
718  //@}
719 
720  /** \name Solving the Dirac equations (protected)
721  */
722  //@{
723  /** \brief Solve the Dirac equations
724 
725  Solves the Dirac equation in from 12 fm to the match point and
726  then out from .04 fm and adjusts eigenvalue with
727  \f[
728  \Delta \varepsilon = -g(r=\mathrm{match\_point})
729  \times (f^{+}-f^{-})
730  \f]
731  */
732  void dirac(int ilevel);
733 
734  /// Take a step in the Dirac equations
735  void dirac_step(double &x, double h, double eigen,
736  double kappa, ubvector &varr);
737 
738  /// The form of the Dirac equations for the ODE stepper
739  int odefun(double x, size_t nv, const ubvector &y,
740  ubvector &dydx, odparms &op);
741 
742  /// Compute the fields for the Dirac equations
743  void field(double x, double &s, double &v, ubvector &varr);
744 
745  /// The default stepper
746  ode_rkck_gsl<ubvector,ubvector,ubvector,
748 
749  /// The ODE stepper
750  ode_step<ubvector,ubvector,ubvector,
752 
753  /** \brief Integrate the Dirac equations using a simple
754  inline 4th order Runge-Kutta
755  */
756  double dirac_rk4(double x, double g1, double f1, double &funt,
757  double eigen, double kappa, ubvector &varr);
758  //@}
759 
760  /// True if init() has been called
762 
763  /// ODE functions
764  ubvector ode_y;
765 
766  /// ODE derivatives
767  ubvector ode_dydx;
768 
769  /// ODE errors
770  ubvector ode_yerr;
771 
772  /// \name Gauss-Legendre integration points and weights
773  //@{
774  double x12[6], w12[6];
775  double x100[50], w100[50];
776  //@}
777 
778 #endif
779 
780  };
781 
782 #ifndef DOXYGEN_NO_O2NS
783 }
784 #endif
785 
786 #endif
void pfold(double x, double &xrhof)
Fold in proton form factor.
void dirac(int ilevel)
Solve the Dirac equations.
double dirac_tol
Tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:465
eos_had_rmf * rmf
The base EOS.
Definition: nucleus_rmf.h:524
double stens
Surface tension (in )
Definition: nucleus_rmf.h:352
int twojp1
Degeneracy .
Definition: nucleus_rmf.h:267
bool generic_ode
If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta.
Definition: nucleus_rmf.h:509
A convenient struct for the solution of the Dirac equations.
Definition: nucleus_rmf.h:239
double energy
Energy eigenvalue.
Definition: nucleus_rmf.h:271
double omega0
Vector field at r=0.
Definition: nucleus_rmf.h:489
double rho_rhs(double sig, double ome, double rho)
Isubvector density RHS.
void meson_init()
Calculate meson and photon Green&#39;s functions gin and gout.
double rprms
Proton RMS radius (in fm)
Definition: nucleus_rmf.h:373
shell proton_shells[n_internal_levels]
The starting proton levels.
Definition: nucleus_rmf.h:547
ubvector xrhor
The isubvector field RHS.
Definition: nucleus_rmf.h:677
ubmatrix field0
Values of the fields from the last iteration.
Definition: nucleus_rmf.h:562
Spherical closed-shell nuclei with a relativistic mean-field model in the Hartree approximation...
Definition: nucleus_rmf.h:229
double etot
Total energy (in MeV)
Definition: nucleus_rmf.h:380
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: nucleus_rmf.h:442
std::shared_ptr< table_units<> > profiles
The radial profiles.
Definition: nucleus_rmf.h:528
std::shared_ptr< table_units<> > get_profiles()
Get the radial profiles.
Definition: nucleus_rmf.h:324
static const int grid_size
The grid size.
Definition: nucleus_rmf.h:550
ubvector chdenc
Charge density.
Definition: nucleus_rmf.h:683
double A0
Coulomb field at r=0.
Definition: nucleus_rmf.h:493
double step_size
The grid step size (default 0.04)
Definition: nucleus_rmf.h:553
double r_charge_cm
Charge radius corrected by the center of mass (in fm)
Definition: nucleus_rmf.h:392
ode_rkck_gsl< ubvector, ubvector, ubvector, ode_funct11 > def_step
The default stepper.
Definition: nucleus_rmf.h:747
int nlevels
The number of levels.
Definition: nucleus_rmf.h:328
double sigma0
Scalar field at r=0.
Definition: nucleus_rmf.h:487
void meson_solve()
Solve for the meson profiles.
A shell of nucleons for nucleus_rmf.
Definition: nucleus_rmf.h:265
void init_meson_density()
Initialize the meson and photon fields, the densities, etc.
ubmatrix gin
The Green&#39;s functions inside.
Definition: nucleus_rmf.h:568
ode_step< ubvector, ubvector, ubvector, ode_funct11 > * ostep
The ODE stepper.
Definition: nucleus_rmf.h:751
double dirac_tol2
Second tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:470
ubvector xrhos
The scalar field RHS.
Definition: nucleus_rmf.h:671
double match_point
Matching radius (in fm)
Definition: nucleus_rmf.h:277
double fermi_radius
The radius for which the fields are half their central value.
Definition: nucleus_rmf.h:495
std::vector< shell > levels
The levels (protons first, then neutrons)
Definition: nucleus_rmf.h:337
void meson_iter(double ic)
Calculate meson and photon fields.
Relativistic mean field theory EOS.
Definition: eos_had_rmf.h:295
int set_eos(eos_had_rmf &r)
Set the base EOS to be used.
Definition: nucleus_rmf.h:410
interp_vec< ubvector > * gi
Interpolation object.
Definition: nucleus_rmf.h:717
double rnrp
Skin thickness (in fm)
Definition: nucleus_rmf.h:359
shell neutron_shells[n_internal_levels]
The starting neutron levels.
Definition: nucleus_rmf.h:544
double meson_tol
Tolerance for meson field equations (default ).
Definition: nucleus_rmf.h:473
ubmatrix fields
The values of the fields.
Definition: nucleus_rmf.h:565
ubvector energy
Energy integrand.
Definition: nucleus_rmf.h:690
bool init_called
True if init() has been called.
Definition: nucleus_rmf.h:761
double xpform(double x, double xp, double a)
Function representing proton form factor.
int itmax
Maximum number of total iterations (default 70)
Definition: nucleus_rmf.h:452
eos_had_rmf def_rmf
The default equation of state (default NL3)
Definition: nucleus_rmf.h:403
int load_nl3(eos_had_rmf &r)
Load the default model NL3 into the given eos_had_rmf object.
void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
Initialize a run.
ubvector xrhosp
The proton scalar density times radius squared.
Definition: nucleus_rmf.h:668
void set_step(ode_step< ubvector, ubvector, ubvector, ode_funct11 > &step)
Set the stepper for the Dirac differential equation.
Definition: nucleus_rmf.h:445
double rho0
Isubvector field at r=0.
Definition: nucleus_rmf.h:491
void center_mass_corr(double atot)
Compute the center of mass correction.
int odefun(double x, size_t nv, const ubvector &y, ubvector &dydx, odparms &op)
The form of the Dirac equations for the ODE stepper.
fermion n
The neutron.
Definition: nucleus_rmf.h:426
initial_guess ig
Parameters for initial guess.
Definition: nucleus_rmf.h:504
double fermi_width
The "width" of the Fermi-Dirac function.
Definition: nucleus_rmf.h:497
double eigen
Eigenvalue.
Definition: nucleus_rmf.h:241
std::vector< shell > * levp
A pointer to the current vector of levels (either levels or unocc_levels)
Definition: nucleus_rmf.h:537
double a_proton
The parameter for the charge density of the proton (default is about 4.27073)
Definition: nucleus_rmf.h:518
void gauss(double xmin, double xmax, double x, double &xi)
Perform integrations for form factor.
double sigma_rhs(double sig, double ome, double rho)
Scalar density RHS.
double omega_rhs(double sig, double ome, double rho)
Vector density RHS.
ubvector arho
Baryon density.
Definition: nucleus_rmf.h:686
int surf_index
The grid index corresponding to the nuclear surface (computed by init_run())
Definition: nucleus_rmf.h:658
int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)
After convergence, make CM corrections, etc.
int meson_itmax
Maximum number of iterations for solving the meson field equations (default 10000) ...
Definition: nucleus_rmf.h:457
ubvector xrhov
The vector field RHS.
Definition: nucleus_rmf.h:674
ubvector ode_dydx
ODE derivatives.
Definition: nucleus_rmf.h:767
void field(double x, double &s, double &v, ubvector &varr)
Compute the fields for the Dirac equations.
thermo hb
thermo object for the EOS
Definition: nucleus_rmf.h:419
double isospin
Isospin ( or .
Definition: nucleus_rmf.h:273
Initial guess structure.
Definition: nucleus_rmf.h:485
ubvector chden1
Charge density.
Definition: nucleus_rmf.h:680
ubvector ode_yerr
ODE errors.
Definition: nucleus_rmf.h:770
int dirac_itmax
Maximum number of iterations for solving the Dirac equations (default 100)
Definition: nucleus_rmf.h:462
void set_verbose(int v)
Set output level.
Definition: nucleus_rmf.h:297
double rnrms
Neutron RMS radius (in fm)
Definition: nucleus_rmf.h:366
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.
ubvector ode_y
ODE functions.
Definition: nucleus_rmf.h:764
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.
std::string state
Angular momentum-spin state .
Definition: nucleus_rmf.h:275
void energies(double xpro, double xnu, double e)
Calculate the energy profile.
double r_charge
Charge radius (in fm)
Definition: nucleus_rmf.h:386
double xrhop(double x1)
Desc.
int nodes
Number of nodes in the wave function.
Definition: nucleus_rmf.h:283
double mnuc
The nucleon mass (automatically set in init_fun())
Definition: nucleus_rmf.h:556
fermion p
The proton.
Definition: nucleus_rmf.h:433
void dirac_step(double &x, double h, double eigen, double kappa, ubvector &varr)
Take a step in the Dirac equations.
ubmatrix xrho
The densities times radius squared.
Definition: nucleus_rmf.h:665
std::shared_ptr< table_units<> > chden_table
The final charge densities.
Definition: nucleus_rmf.h:532
int nuolevels
The number of unoccupied levels (equal to unocc_Z + unocc_N)
Definition: nucleus_rmf.h:340
int verbose
Control output (default 1)
Definition: nucleus_rmf.h:541
std::vector< shell > unocc_levels
The unoccupied levels (protons first, then neutrons)
Definition: nucleus_rmf.h:346
double kappa
Quantum number .
Definition: nucleus_rmf.h:243
static const int n_internal_levels
The total number of shells stored internally.
Definition: nucleus_rmf.h:251
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct11
void iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N, int &iconverged)
Perform an iteration.
std::shared_ptr< table_units<> > get_chden()
The final charge densities.
Definition: nucleus_rmf.h:328
ubmatrix * fields
The meson fields.
Definition: nucleus_rmf.h:245
ubmatrix gout
The Green&#39;s functions outside.
Definition: nucleus_rmf.h:571

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