nucleus_rmf.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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 I believe currently the center of mass correction
207  for the binding energy per nucleon is not done and has
208  to be added in after the fact
209 
210  \future Sort energy levels at the end by energy
211  \future Improve the numerical methods
212  \future Make the neutron and proton orbitals more configurable
213  \future Generalize to \f$ m_n \neq m_p \f$ .
214  \future Allow more freedom in the integrations
215  \future Consider converting everything to inverse fermis.
216  \future Convert to zero-indexed arrays (mostly done)
217  \future Warn when the level ordering is wrong, and unoccupied levels
218  are lower energy than occupied levels
219  \future Connect with \ref o2scl::nucmass ?
220  */
221  class nucleus_rmf {
222 
225 
226 #ifndef DOXYGEN_INTERNAL
227 
228  protected:
229 
230  /// A convenient struct for the solution of the Dirac equations
231  typedef struct {
232  /// Eigenvalue
233  double eigen;
234  /// Quantum number \f$ \kappa \f$ .
235  double kappa;
236  /// The meson fields
237  ubmatrix *fields;
238  /// Desc
239  ubvector *varr;
240  } odparms;
241 
242  /// The total number of shells stored internally
243  static const int n_internal_levels=29;
244 
245 #endif
246 
247  public:
248 
249  nucleus_rmf();
250 
251  ~nucleus_rmf();
252 
253  /** \name Basic operation
254  */
255  //@{
256  /// A shell of nucleons for \ref nucleus_rmf
257  typedef struct {
258  /// Degeneracy \f$ 2 j+1 \f$ .
259  int twojp1;
260  /// \f$ \kappa \f$
261  int kappa;
262  /// Energy eigenvalue
263  double energy;
264  /// Isospin ( \f$ +1/2 \f$ or \f$ -1/2 \f$ .
265  double isospin;
266  /// Angular momentum-spin state \f$ ^{2s+1} \ell_{j} \f$
267  std::string state;
268  /// Matching radius (in fm)
269  double match_point;
270  /// Desc.
271  double eigen;
272  /// Desc.
273  double eigenc;
274  /// Number of nodes in the wave function
275  int nodes;
276  } shell;
277 
278  /** \brief Computes the structure of a nucleus with the specified
279  number of levels
280 
281  Note that \ref rmf must be set before run_nucleus() is called.
282 
283  This calls init_run(), and then iterate() until \c iconverged is
284  1, and then post_converge().
285  */
286  int run_nucleus(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
287 
288  /// Set output level
289  void set_verbose(int v) { verbose=v; };
290  //@}
291 
292  /** \name Lower-level interface
293  */
294  //@{
295  /** \brief Initialize a run
296 
297  Note that \ref rmf must be set before run_nucleus() is called.
298  */
299  void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
300 
301  /// Perform an iteration
302  int iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N,
303  int &iconverged, int &dirac_converged, int &meson_converged);
304 
305  /// After convergence, make CM corrections, etc.
306  int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N);
307 
308  //@}
309 
310  /// \name Results
311  //@{
312  /** \brief Get the radial profiles
313 
314  The profiles are calculated each iteration by iterate().
315  */
316  std::shared_ptr<table_units<> > get_profiles() { return profiles; };
317 
318  /** \brief The final charge densities
319  */
320  std::shared_ptr<table_units<> > get_chden() { return chden_table; };
321 
322  /// The number of levels
323  int nlevels;
324 
325  /** \brief The levels (protons first, then neutrons)
326 
327  An array of size \ref nlevels
328  */
329  std::vector<shell> levels;
330 
331  /// The number of unoccupied levels (equal to \c unocc_Z + \c unocc_N)
333 
334  /** \brief The unoccupied levels (protons first, then neutrons)
335 
336  An array of size \ref nuolevels
337  */
338  std::vector<shell> unocc_levels;
339 
340  /** \brief Surface tension (in \f$ \mathrm{fm}^{-3} \f$ )
341 
342  Computed in post_converge() or automatically in run_nucleus()
343  */
344  double stens;
345 
346  /** \brief Skin thickness (in fm)
347 
348  Computed every iteration in iterate()
349  or automatically in run_nucleus()
350  */
351  double rnrp;
352 
353  /** \brief Neutron RMS radius (in fm)
354 
355  Computed every iteration in iterate() or automatically in
356  run_nucleus()
357  */
358  double rnrms;
359 
360  /** \brief Proton RMS radius (in fm)
361 
362  Computed every iteration in iterate() or automatically in
363  run_nucleus()
364  */
365  double rprms;
366 
367  /** \brief Total energy (in MeV)
368 
369  Computed every iteration in iterate() or automatically in
370  run_nucleus()
371  */
372  double etot;
373 
374  /** \brief Charge radius (in fm)
375 
376  Computed in post_converge() or automatically in run_nucleus()
377  */
378  double r_charge;
379 
380  /** \brief Charge radius corrected by the center of mass (in fm)
381 
382  Computed in post_converge() or automatically in run_nucleus()
383  */
384  double r_charge_cm;
385  //@}
386 
387  /** \name Equation of state
388  */
389  //@{
390  /** \brief The default equation of state (default NL3)
391 
392  This is set in the constructor to be the default
393  model, NL3, using the function \ref load_nl3().
394  */
396 
397  /** \brief Set the base EOS to be used
398 
399  The equation of state must be set before run_nucleus() or
400  init_fun() are called, including the value of eos_had_rmf::mnuc.
401  */
403  rmf=&r;
404  return 0;
405  }
406 
407  /** \brief \ref thermo object for the EOS
408 
409  This is just used as temporary storage.
410  */
412 
413  /** \brief The neutron
414 
415  The mass of the neutron is ignored and set by init_run()
416  to be eos_had_rmf::mnuc from \ref rmf.
417  */
419 
420  /** \brief The proton
421 
422  The mass of the proton is ignored and set by init_run()
423  to be eos_had_rmf::mnuc from \ref rmf.
424  */
426  //@}
427 
428  /** \name Numeric configuration
429  */
430  //@{
431  /** \brief If true, call the error handler if the routine does not
432  converge or reach the desired tolerance (default true)
433  */
435 
436  /// Set the stepper for the Dirac differential equation
437  void set_step(ode_step<ubvector,ubvector,ubvector,
438  ode_funct> &step) {
439  ostep=&step;
440  return;
441  }
442 
443  /// Maximum number of total iterations (default 70)
444  int itmax;
445 
446  /** \brief Maximum number of iterations for solving the meson
447  field equations (default 10000)
448  */
450 
451  /** \brief Maximum number of iterations for solving the Dirac
452  equations (default 100)
453  */
455 
456  /// Tolerance for Dirac equations (default \f$ 5 \times 10^{-3} \f$ ).
457  double dirac_tol;
458 
459  /** \brief Second tolerance for Dirac equations
460  (default \f$ 5 \times 10^{-4} \f$ ).
461  */
462  double dirac_tol2;
463 
464  /// Tolerance for meson field equations (default \f$ 10^{-6} \f$ ).
465  double meson_tol;
466  //@}
467 
468  /** \brief Initial guess structure
469 
470  The initial guess for the meson field profiles is
471  a set of Fermi-Dirac functions, i.e.
472  \f[
473  \sigma(r)=\mathrm{sigma0}/
474  [1+\exp((r-\mathrm{fermi\_radius})/\mathrm{fermi\_width})]
475  \f]
476  */
477  typedef struct {
478  /// Scalar field at r=0
479  double sigma0;
480  /// Vector field at r=0
481  double omega0;
482  /// Isubvector field at r=0
483  double rho0;
484  /// Coulomb field at r=0
485  double A0;
486  /// The radius for which the fields are half their central value
487  double fermi_radius;
488  /// The "width" of the Fermi-Dirac function
489  double fermi_width;
490  } initial_guess;
491 
492  /** \brief Parameters for initial guess
493 
494  Default is {310,240,-6,25.9,6.85,0.6}
495  */
497 
498  /** \brief If true, use the generic ODE solver instead of the
499  internal 4th order Runge-Kutta
500  */
502 
503 #ifndef DOXYGEN_INTERNAL
504 
505  protected:
506 
507  /** \brief The parameter for the charge density of the proton
508  (default is about 4.27073)
509  */
510  double a_proton;
511 
512  /// Load the default model NL3 into the given \ref eos_had_rmf object
513  int load_nl3(eos_had_rmf &r);
514 
515  /// The base EOS
517 
518  /** \brief The radial profiles
519  */
520  std::shared_ptr<table_units<> > profiles;
521 
522  /** \brief The final charge densities
523  */
524  std::shared_ptr<table_units<> > chden_table;
525 
526  /** \brief A pointer to the current vector of levels
527  (either \ref levels or \ref unocc_levels)
528  */
529  std::vector<shell> *levp;
530 
531  /** \brief Control output (default 1)
532  */
533  int verbose;
534 
535  /// The starting neutron levels
537 
538  /// The starting proton levels
540 
541  /// The grid size
542  static const int grid_size=300;
543 
544  /// The grid step size (default 0.04)
545  double step_size;
546 
547  /// The nucleon mass (automatically set in init_fun())
548  double mnuc;
549 
550  /** \name The meson and photon fields and field equations (protected)
551  */
552  //@{
553  /// Values of the fields from the last iteration
554  ubmatrix field0;
555 
556  /// The values of the fields
557  ubmatrix fields;
558 
559  /// The Green's functions inside
560  ubmatrix gin;
561 
562  /// The Green's functions outside
563  ubmatrix gout;
564 
565  /// Scalar density RHS
566  double sigma_rhs(double sig, double ome, double rho);
567 
568  /// Vector density RHS
569  double omega_rhs(double sig, double ome, double rho);
570 
571  /// Isubvector density RHS
572  double rho_rhs(double sig, double ome, double rho);
573 
574  /** \brief Calculate meson and photon
575  Green's functions \ref gin and \ref gout
576  */
577  void meson_init();
578 
579  /** \brief Calculate meson and photon fields
580 
581  The meson and photon field equations are of the
582  Sturm-Liouville form, e.g.
583  \f[
584  \left[r^2 \sigma^{\prime}(r) \right]^{\prime} -
585  r^2 m_{\sigma}^2 \sigma(r) = f(r)
586  \f]
587  where \f$ \sigma(0) = \sigma_0 \f$ and \f$ \sigma(+\infty) = 0
588  \f$. A solution of the homogeneous equation with \f$ f(r) =0
589  \f$ is \f$ \sigma(r) = \sigma_0 \sinh( m_{\sigma} r)/
590  (m_{\sigma} r) \f$. The associated Green's function is
591  \f[
592  D(r,r^{\prime},m_{\sigma}) = \frac{-1}{m_{\sigma} r r^{\prime}}
593  \sinh (m_{\sigma} r_{<}) \exp (-m_{\sigma} r_{>})
594  \f]
595  where \f$ r_{<} \f$ is the smaller of \f$ r \f$ and
596  \f$ r^{\prime} \f$ and \f$ r_{>} \f$ is the larger.
597  One can then write the solution of the meson field
598  equation given the density
599  \f[
600  \sigma(r) = \int_0^{\infty} r^{\prime 2}~d r^{\prime}
601  \left[ - g_{\sigma} n_{s}(r) \right]
602  D\left(r,r^{\prime},m_{\sigma}\right)
603  \f]
604 
605  Since radii are positive, \f$ \sinh (m r) \approx
606  e^{m r}/2 \f$ and
607  \f[
608  D(r,r^{\prime},m_{\sigma}) \approx \left[
609  \frac{-1}{2 m_{\sigma} r_{<}}
610  \exp (m_{\sigma} r_{<})
611  \right] \left[
612  \frac{1}{r_{>}}
613  \exp (-m_{\sigma} r_{>})
614  \right]
615  \f]
616  The two terms in the Green's function are separated into
617  \f[
618  \mathrm{gin}(r) = \frac{e^{m_{\sigma} r}}{2 m_{\sigma} r}
619  \f]
620  and
621  \f[
622  \mathrm{gout}(r) = \frac{e^{-m_{\sigma} r}}{r} \, .
623  \f]
624  These functions are computed in \ref meson_init() . Then the
625  field is given by
626  \f[
627  \sigma(r)= \left(\frac{e^{-m_{\sigma} r}}{r}\right)
628  \int_0^{r} r^{\prime 2} g_{\sigma} n_{s}
629  \left(\frac{e^{m_{\sigma} r^{\prime}}}{2 m_{\sigma}
630  r^{\prime}} \right)~d r^{\prime} +
631  \left(\frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \right)
632  \int_r^{\infty} r^{\prime 2} g_{\sigma} n_{s}
633  \left(\frac{e^{-m_{\sigma} r^{\prime}}}
634  {r^{\prime}}\right)~d r^{\prime}
635  \f]
636  where the first integral is stored in <tt>xi2</tt> and
637  the second is in <tt>xi1</tt> in the function \ref meson_iter() .
638  The part of <tt>xi2</tt> at \f$ r^{\prime}=0 \f$ is
639  stored in <tt>xi20</tt>.
640  */
641  void meson_iter(double ic);
642 
643  /** \brief Solve for the meson profiles
644  */
645  int meson_solve();
646 
647  /** \brief The grid index corresponding to the nuclear surface
648  (computed by init_run())
649  */
651  //@}
652 
653  /** \name Density information (protected)
654  */
655  //@{
656  /// The densities times radius squared
657  ubmatrix xrho;
658 
659  /// The proton scalar density times radius squared
660  ubvector xrhosp;
661 
662  /// The scalar field RHS
663  ubvector xrhos;
664 
665  /// The vector field RHS
666  ubvector xrhov;
667 
668  /// The isubvector field RHS
669  ubvector xrhor;
670 
671  /// Charge density
672  ubvector chden1;
673 
674  /// Charge density
675  ubvector chdenc;
676 
677  /// Baryon density
678  ubvector arho;
679  //@}
680 
681  /// Energy integrand
682  ubvector energy;
683 
684  /// Initialize the meson and photon fields, the densities, etc.
685  void init_meson_density();
686 
687  /// Calculate the energy profile
688  int energy_radii(double xpro, double xnu, double e);
689 
690  /// Compute the center of mass correction
691  void center_mass_corr(double atot);
692 
693  /** \name Calculating the form factor, etc. (protected)
694  */
695  //@{
696  /// Fold in proton form factor
697  void pfold(double x, double &xrhof);
698 
699  /// Function representing proton form factor
700  double xpform(double x, double xp, double a);
701 
702  /// Perform integrations for form factor
703  void gauss(double xmin, double xmax, double x, double &xi);
704 
705  /// Desc.
706  double xrhop(double x1);
707 
708  /// Interpolation object
710  //@}
711 
712  /** \name Solving the Dirac equations (protected)
713  */
714  //@{
715  /** \brief Solve the Dirac equations
716 
717  Solves the Dirac equation in from 12 fm to the match point and
718  then out from .04 fm and adjusts eigenvalue with
719  \f[
720  \Delta \varepsilon = -g(r=\mathrm{match\_point})
721  \times (f^{+}-f^{-})
722  \f]
723  */
724  int dirac(int ilevel);
725 
726  /// Take a step in the Dirac equations
727  void dirac_step(double &x, double h, double eigen,
728  double kappa, ubvector &varr);
729 
730  /// The form of the Dirac equations for the ODE stepper
731  int odefun(double x, size_t nv, const ubvector &y,
732  ubvector &dydx, odparms &op);
733 
734  /// Compute the fields for the Dirac equations
735  void field(double x, double &s, double &v, ubvector &varr);
736 
737  /// The default stepper
738  ode_rkck_gsl<ubvector,ubvector,ubvector,
740 
741  /// The ODE stepper
742  ode_step<ubvector,ubvector,ubvector,
744 
745  /** \brief Integrate the Dirac equations using a simple
746  inline 4th order Runge-Kutta
747  */
748  double dirac_rk4(double x, double g1, double f1, double &funt,
749  double eigen, double kappa, ubvector &varr);
750  //@}
751 
752  /// True if init() has been called
754 
755  /// ODE functions
756  ubvector ode_y;
757 
758  /// ODE derivatives
759  ubvector ode_dydx;
760 
761  /// ODE errors
762  ubvector ode_yerr;
763 
764  /// \name Gauss-Legendre integration points and weights
765  //@{
766  double x12[6], w12[6];
767  double x100[50], w100[50];
768  //@}
769 
770 #endif
771 
772  };
773 
774 #ifndef DOXYGEN_NO_O2NS
775 }
776 #endif
777 
778 #endif
void pfold(double x, double &xrhof)
Fold in proton form factor.
double dirac_tol
Tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:457
eos_had_rmf * rmf
The base EOS.
Definition: nucleus_rmf.h:516
double stens
Surface tension (in )
Definition: nucleus_rmf.h:344
int twojp1
Degeneracy .
Definition: nucleus_rmf.h:259
bool generic_ode
If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta.
Definition: nucleus_rmf.h:501
A convenient struct for the solution of the Dirac equations.
Definition: nucleus_rmf.h:231
int energy_radii(double xpro, double xnu, double e)
Calculate the energy profile.
double energy
Energy eigenvalue.
Definition: nucleus_rmf.h:263
double omega0
Vector field at r=0.
Definition: nucleus_rmf.h:481
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:365
shell proton_shells[n_internal_levels]
The starting proton levels.
Definition: nucleus_rmf.h:539
ubvector xrhor
The isubvector field RHS.
Definition: nucleus_rmf.h:669
ubmatrix field0
Values of the fields from the last iteration.
Definition: nucleus_rmf.h:554
Spherical closed-shell nuclei with a relativistic mean-field model in the Hartree approximation...
Definition: nucleus_rmf.h:221
int iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N, int &iconverged, int &dirac_converged, int &meson_converged)
Perform an iteration.
double etot
Total energy (in MeV)
Definition: nucleus_rmf.h:372
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:434
std::shared_ptr< table_units<> > profiles
The radial profiles.
Definition: nucleus_rmf.h:520
std::shared_ptr< table_units<> > get_profiles()
Get the radial profiles.
Definition: nucleus_rmf.h:316
static const int grid_size
The grid size.
Definition: nucleus_rmf.h:542
ubvector chdenc
Charge density.
Definition: nucleus_rmf.h:675
double A0
Coulomb field at r=0.
Definition: nucleus_rmf.h:485
double step_size
The grid step size (default 0.04)
Definition: nucleus_rmf.h:545
double r_charge_cm
Charge radius corrected by the center of mass (in fm)
Definition: nucleus_rmf.h:384
int nlevels
The number of levels.
Definition: nucleus_rmf.h:320
double sigma0
Scalar field at r=0.
Definition: nucleus_rmf.h:479
A shell of nucleons for nucleus_rmf.
Definition: nucleus_rmf.h:257
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:560
double dirac_tol2
Second tolerance for Dirac equations (default ).
Definition: nucleus_rmf.h:462
ubvector xrhos
The scalar field RHS.
Definition: nucleus_rmf.h:663
double match_point
Matching radius (in fm)
Definition: nucleus_rmf.h:269
double fermi_radius
The radius for which the fields are half their central value.
Definition: nucleus_rmf.h:487
std::vector< shell > levels
The levels (protons first, then neutrons)
Definition: nucleus_rmf.h:329
void meson_iter(double ic)
Calculate meson and photon fields.
Relativistic mean field theory EOS.
Definition: eos_had_rmf.h:298
int set_eos(eos_had_rmf &r)
Set the base EOS to be used.
Definition: nucleus_rmf.h:402
interp_vec< ubvector > * gi
Interpolation object.
Definition: nucleus_rmf.h:709
double rnrp
Skin thickness (in fm)
Definition: nucleus_rmf.h:351
shell neutron_shells[n_internal_levels]
The starting neutron levels.
Definition: nucleus_rmf.h:536
double meson_tol
Tolerance for meson field equations (default ).
Definition: nucleus_rmf.h:465
ubmatrix fields
The values of the fields.
Definition: nucleus_rmf.h:557
ubvector energy
Energy integrand.
Definition: nucleus_rmf.h:682
bool init_called
True if init() has been called.
Definition: nucleus_rmf.h:753
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:444
ode_rkck_gsl< ubvector, ubvector, ubvector, ode_funct > def_step
The default stepper.
Definition: nucleus_rmf.h:739
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
eos_had_rmf def_rmf
The default equation of state (default NL3)
Definition: nucleus_rmf.h:395
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:660
void set_step(ode_step< ubvector, ubvector, ubvector, ode_funct > &step)
Set the stepper for the Dirac differential equation.
Definition: nucleus_rmf.h:437
double rho0
Isubvector field at r=0.
Definition: nucleus_rmf.h:483
void center_mass_corr(double atot)
Compute the center of mass correction.
int meson_solve()
Solve for the meson profiles.
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:418
initial_guess ig
Parameters for initial guess.
Definition: nucleus_rmf.h:496
int dirac(int ilevel)
Solve the Dirac equations.
double fermi_width
The "width" of the Fermi-Dirac function.
Definition: nucleus_rmf.h:489
double eigen
Eigenvalue.
Definition: nucleus_rmf.h:233
std::vector< shell > * levp
A pointer to the current vector of levels (either levels or unocc_levels)
Definition: nucleus_rmf.h:529
double a_proton
The parameter for the charge density of the proton (default is about 4.27073)
Definition: nucleus_rmf.h:510
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:678
int surf_index
The grid index corresponding to the nuclear surface (computed by init_run())
Definition: nucleus_rmf.h:650
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:449
ubvector xrhov
The vector field RHS.
Definition: nucleus_rmf.h:666
ubvector ode_dydx
ODE derivatives.
Definition: nucleus_rmf.h:759
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:411
double isospin
Isospin ( or .
Definition: nucleus_rmf.h:265
Initial guess structure.
Definition: nucleus_rmf.h:477
ubvector chden1
Charge density.
Definition: nucleus_rmf.h:672
ubvector ode_yerr
ODE errors.
Definition: nucleus_rmf.h:762
int dirac_itmax
Maximum number of iterations for solving the Dirac equations (default 100)
Definition: nucleus_rmf.h:454
ode_step< ubvector, ubvector, ubvector, ode_funct > * ostep
The ODE stepper.
Definition: nucleus_rmf.h:743
void set_verbose(int v)
Set output level.
Definition: nucleus_rmf.h:289
double rnrms
Neutron RMS radius (in fm)
Definition: nucleus_rmf.h:358
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:756
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:267
double r_charge
Charge radius (in fm)
Definition: nucleus_rmf.h:378
double xrhop(double x1)
Desc.
int nodes
Number of nodes in the wave function.
Definition: nucleus_rmf.h:275
double mnuc
The nucleon mass (automatically set in init_fun())
Definition: nucleus_rmf.h:548
fermion p
The proton.
Definition: nucleus_rmf.h:425
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:657
std::shared_ptr< table_units<> > chden_table
The final charge densities.
Definition: nucleus_rmf.h:524
int nuolevels
The number of unoccupied levels (equal to unocc_Z + unocc_N)
Definition: nucleus_rmf.h:332
int verbose
Control output (default 1)
Definition: nucleus_rmf.h:533
std::vector< shell > unocc_levels
The unoccupied levels (protons first, then neutrons)
Definition: nucleus_rmf.h:338
double kappa
Quantum number .
Definition: nucleus_rmf.h:235
static const int n_internal_levels
The total number of shells stored internally.
Definition: nucleus_rmf.h:243
std::shared_ptr< table_units<> > get_chden()
The final charge densities.
Definition: nucleus_rmf.h:320
ubmatrix * fields
The meson fields.
Definition: nucleus_rmf.h:237
ubmatrix gout
The Green&#39;s functions outside.
Definition: nucleus_rmf.h:563

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