#include <potentials.h>

A class that stores parameters of potentials, calculates potentials and their gradients.

Potentials are responsible for long-range interactions and stand in the left part of Boltzmann equation. Short-range interactions are taken into account in the right part of it - in the collision term.

Definition at line 36 of file potentials.h.

Collaboration diagram for smash::Potentials:

[legend]

Public Member Functions
	Potentials (Configuration conf, const DensityParameters &parameters)
	Potentials constructor. More...

virtual	~Potentials ()
	Standard destructor. More...

ThreeVector	single_particle_energy_gradient (DensityLattice *jB_lattice, const ThreeVector &position, const ThreeVector &momentum, double mass, ParticleList &plist) const
	Calculates the gradient of the single-particle energy (including potentials) in the calculation frame in MeV/fm. More...

double	calculation_frame_energy (const ThreeVector &momentum, const FourVector &jmu_B, double mass) const
	Evaluates the single-particle energy (including the potential) of a particle at a given position and momentum in the calculation frame. More...

double	skyrme_pot (const double baryon_density) const
	Evaluates Skyrme potential given a baryon density. More...

double	symmetry_pot (const double baryon_isospin_density, const double baryon_density) const
	Evaluates symmetry potential given baryon isospin density. More...

double	symmetry_S (const double baryon_density) const
	Calculate the factor \(S(\rho)\) in the symmetry potential. More...

FourVector	vdf_pot (double rhoB, const FourVector jmuB_net) const
	Evaluates the FourVector potential in the VDF model given the rest frame density and the computational frame baryon current. More...

double	potential (const ThreeVector &r, const ParticleList &plist, const ParticleType &acts_on) const
	Evaluates potential (Skyrme with optional Symmetry or VDF) at point r. More...

std::pair< ThreeVector, ThreeVector >	skyrme_force (const double rhoB, const ThreeVector grad_j0B, const ThreeVector dvecjB_dt, const ThreeVector curl_vecjB) const
	Evaluates the electric and magnetic components of the skyrme force. More...

std::pair< ThreeVector, ThreeVector >	symmetry_force (const double rhoI3, const ThreeVector grad_j0I3, const ThreeVector dvecjI3_dt, const ThreeVector curl_vecjI3, const double rhoB, const ThreeVector grad_j0B, const ThreeVector dvecjB_dt, const ThreeVector curl_vecjB) const
	Evaluates the electric and magnetic components of the symmetry force. More...

std::pair< ThreeVector, ThreeVector >	vdf_force (double rhoB, const double drhoB_dt, const ThreeVector grad_rhoB, const ThreeVector gradrhoB_cross_vecjB, const double j0B, const ThreeVector grad_j0B, const ThreeVector vecjB, const ThreeVector dvecjB_dt, const ThreeVector curl_vecjB) const
	Evaluates the electric and magnetic components of force in the VDF model given the derivatives of the baryon current \(j^{\mu}\). More...

std::pair< ThreeVector, ThreeVector >	vdf_force (const ThreeVector grad_A_0, const ThreeVector dA_dt, const ThreeVector curl_vecA) const
	Evaluates the electric and magnetic components of force in the VDF force given the derivatives of the VDF mean-field \(A^\mu\). More...

virtual std::tuple< ThreeVector, ThreeVector, ThreeVector, ThreeVector >	all_forces (const ThreeVector &r, const ParticleList &plist) const
	Evaluates the electric and magnetic components of the forces at point r. More...

virtual bool	use_skyrme () const

virtual bool	use_symmetry () const

virtual bool	use_coulomb () const

virtual bool	use_momentum_dependence () const

double	skyrme_a () const

double	skyrme_b () const

double	skyrme_tau () const

double	symmetry_S_pot () const

virtual bool	use_vdf () const

double	saturation_density () const

const std::vector< double > &	coeffs () const

const std::vector< double > &	powers () const

int	number_of_terms () const

double	coulomb_r_cut () const

bool	use_potentials_outside_lattice () const

Static Public Member Functions
static std::pair< double, int >	force_scale (const ParticleType &data)
	Evaluates the scaling factor of the forces acting on the particles. More...

static ThreeVector	E_field_integrand (ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
	Integrand for calculating the electric field. More...

static ThreeVector	B_field_integrand (ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
	Integrand for calculating the magnetic field using the Biot-Savart formula. More...

Private Member Functions
double	dVsym_drhoI3 (const double rhoB, const double rhoI3) const
	Calculate the derivative of the symmetry potential with respect to the isospin density in GeV * fm^3. More...

double	dVsym_drhoB (const double rhoB, const double rhoI3) const
	Calculate the derivative of the symmetry potential with respect to the net baryon density in GeV * fm^3. More...

Static Private Member Functions
static double	skyrme_pot (const double baryon_density, const double A, const double B, const double tau)
	Single particle Skyrme potential in MeV. More...

static double	root_eq_potentials (double energy_calc, const ThreeVector &momentum_calc, const FourVector &jmu, double m, double A, double B, double tau, double C, double Lambda)
	Root equation used to determine the energy in the calculation frame. More...

static double	momentum_dependent_part (double momentum, double rho, double C, double Lambda)
	Momentum dependent term of the potential. More...

Private Attributes
const DensityParameters	param_
	Struct that contains the gaussian smearing width \(\sigma\), the distance cutoff \(r_{\rm cut}\) and the testparticle number needed for the density calculation. More...

bool	use_skyrme_
	Skyrme potential on/off. More...

bool	use_symmetry_
	Symmetry potential on/off. More...

bool	use_coulomb_
	Coulomb potential on/Off. More...

bool	use_vdf_
	VDF potential on/off. More...

bool	use_momentum_dependence_
	Momentum-dependent part on/off. More...

double	skyrme_a_
	Parameter of skyrme potentials: the coefficient in front of \(\frac{\rho}{\rho_0}\) in GeV. More...

double	skyrme_b_
	Parameters of skyrme potentials: the coefficient in front of \((\frac{\rho}{\rho_0})^\tau\) in GeV. More...

double	skyrme_tau_
	Parameters of skyrme potentials: the power index. More...

double	mom_dependence_Lambda_
	Parameter Lambda of the momentum-dependent part of the potentials given in 1/fm. More...

double	mom_dependence_C_
	Parameter C of the momentum-dependent part of the potentials given in MeV. More...

double	symmetry_S_Pot_
	Parameter S_Pot in the symmetry potential in MeV. More...

bool	symmetry_is_rhoB_dependent_ = false
	Whether the baryon density dependence of the symmetry potential is included. More...

double	symmetry_gamma_
	Power \( \gamma \) in formula for \( S(\rho) \): More...

double	coulomb_r_cut_
	Cutoff in integration for coulomb potential. More...

bool	use_potentials_outside_lattice_
	Wether potentials should be included outside of the lattice. More...

double	saturation_density_
	Saturation density of nuclear matter used in the VDF potential; it may vary between different parameterizations. More...

std::vector< double >	coeffs_
	Parameters of the VDF potential: coefficients \(C_i\), in GeV. More...

std::vector< double >	powers_
	Parameters of the VDF potential: exponents \(b_i\). More...

Constructor & Destructor Documentation

◆ Potentials()

smash::Potentials::Potentials	(	Configuration	conf,
		const DensityParameters &	parameters
	)

Potentials constructor.

Parameters

[in]	conf	Configuration which contains the switches determining whether to turn on the Skyrme or the symmetry potentials, and the coefficients controlling how strong the potentials are.
[in]	parameters	Struct that contains the gaussian smearing factor \(\sigma\), the distance cutoff \(r_{\rm cut}\) and the testparticle number needed for the density calculation.

Definition at line 20 of file potentials.cc.

     : param_(param),
       use_skyrme_(conf.has_section(InputSections::p_skyrme)),
       use_symmetry_(conf.has_section(InputSections::p_symmetry)),
       use_coulomb_(conf.has_section(InputSections::p_coulomb)),
       use_vdf_(conf.has_section(InputSections::p_vdf)),
       use_momentum_dependence_(
           conf.has_section(InputSections::p_momentumDependence)),
       use_potentials_outside_lattice_(
           conf.take(InputKeys::potentials_use_potentials_outside_lattice)) {
   if (use_skyrme_) {
     skyrme_a_ = conf.take(InputKeys::potentials_skyrme_skyrmeA);
     skyrme_b_ = conf.take(InputKeys::potentials_skyrme_skyrmeB);
     skyrme_tau_ = conf.take(InputKeys::potentials_skyrme_skyrmeTau);
   }
   if (use_momentum_dependence_) {
     mom_dependence_Lambda_ =
         conf.take(InputKeys::potentials_momentum_dependence_Lambda);
     mom_dependence_C_ = conf.take(InputKeys::potentials_momentum_dependence_C);
   }
  
   if (use_symmetry_) {
     symmetry_S_Pot_ = conf.take(InputKeys::potentials_symmetry_sPot);
     if (conf.has_value(InputKeys::potentials_symmetry_gamma)) {
       symmetry_gamma_ = conf.take(InputKeys::potentials_symmetry_gamma);
       symmetry_is_rhoB_dependent_ = true;
     }
   }
   if (use_coulomb_) {
     coulomb_r_cut_ = conf.take(InputKeys::potentials_coulomb_rCut);
   }
   if (use_vdf_) {
     saturation_density_ = conf.take(InputKeys::potentials_vdf_satRhoB);
     coeffs_ = conf.take(InputKeys::potentials_vdf_coeffs);
     powers_ = conf.take(InputKeys::potentials_vdf_powers);
     if (coeffs_.size() != powers_.size()) {
       throw std::invalid_argument(
           "The number of coefficients should equal the number of powers.");
     }
     // coefficients are provided in MeV, but the code uses GeV
     std::transform(coeffs_.cbegin(), coeffs_.cend(), coeffs_.begin(),
                    [](double c) { return c * mev_to_gev; });
   }
 }

Here is the call graph for this function:

◆ ~Potentials()

smash::Potentials::~Potentials ( )

virtual

Standard destructor.

Definition at line 65 of file potentials.cc.

65 {}

Member Function Documentation

◆ single_particle_energy_gradient()

ThreeVector smash::Potentials::single_particle_energy_gradient	(	DensityLattice *	jB_lattice,
		const ThreeVector &	position,
		const ThreeVector &	momentum,
		double	mass,
		ParticleList &	plist
	)		const

inline

Calculates the gradient of the single-particle energy (including potentials) in the calculation frame in MeV/fm.

Parameters

jB_lattice	Pointer to the baryon density lattice
position	Position of the particle of interest in fm
momentum	Momentum of the particle of interest in GeV
mass	Mass of the particle of interest in GeV
plist	List of all particles

Returns: ThreeVector gradient of the single particle energy in the calculation frame in MeV/fm

Definition at line 65 of file potentials.h.

                                                                          {
     const std::array<double, 3> dr = (jB_lattice)
                                          ? jB_lattice->cell_sizes()
                                          : std::array<double, 3>{0.1, 0.1, 0.1};
     ThreeVector result, position_left, position_right;
     DensityOnLattice jmu_left, jmu_right;
     FourVector net_4current_left, net_4current_right;
     for (int i = 0; i < 3; i++) {
       position_left = position;
       position_left[i] -= dr[i];
       position_right = position;
       position_right[i] += dr[i];
  
       if (jB_lattice && jB_lattice->value_at(position_left, jmu_left)) {
         net_4current_left = jmu_left.jmu_net();
       } else if (use_potentials_outside_lattice_) {
         auto current = current_eckart(position_left, plist, param_,
                                       DensityType::Baryon, false, true);
         net_4current_left = std::get<1>(current);
       } else {
         return {0., 0., 0.};
       }
  
       if (jB_lattice && jB_lattice->value_at(position_right, jmu_right)) {
         net_4current_right = jmu_right.jmu_net();
       } else if (use_potentials_outside_lattice_) {
         auto current = current_eckart(position_right, plist, param_,
                                       DensityType::Baryon, false, true);
         net_4current_right = std::get<1>(current);
       } else {
         return {0., 0., 0.};
       }
       result[i] =
           (calculation_frame_energy(momentum, net_4current_right, mass) -
            calculation_frame_energy(momentum, net_4current_left, mass)) /
           (2 * dr[i]);
     }
     return result;
   }

Here is the call graph for this function:

◆ calculation_frame_energy()

double smash::Potentials::calculation_frame_energy	(	const ThreeVector &	momentum,
		const FourVector &	jmu_B,
		double	mass
	)		const

inline

Evaluates the single-particle energy (including the potential) of a particle at a given position and momentum in the calculation frame.

Parameters

[in]	momentum	Momentum of interest in GeV
[in]	jmu_B	Baryon current density at pos
[in]	mass	mass of the particle of interest

Returns: the energy of a particle in the calculation frame

Definition at line 118 of file potentials.h.

                                                                               {
     std::function<double(double)> root_equation = [momentum, jmu_B, mass,
                                                    this](double energy) {
       return root_eq_potentials(energy, momentum, jmu_B, mass, skyrme_a_,
                                 skyrme_b_, skyrme_tau_, mom_dependence_C_,
                                 mom_dependence_Lambda_);
     };
     RootSolver1D root_solver{root_equation};
     constexpr std::size_t max_number_root_solver_iterations = 100000;
     const double initial_guess = std::sqrt(mass * mass + momentum * momentum);
     const std::array<double, 4> interval_half_widths = {0.05, 0.5, 5.0, 50.0};
     for (const double half_width : interval_half_widths) {
       const std::pair<double, double> x_range = {initial_guess - half_width,
                                                  initial_guess + half_width};
       const auto calc_frame_energy = root_solver.try_find_root(
           x_range.first, x_range.second, max_number_root_solver_iterations);
       if (calc_frame_energy) {
         return *calc_frame_energy;
       } else {
         logg[LPotentials].debug()
             << "Did not find a root for potentials in the interval ["
             << x_range.first << " GeV ," << x_range.second << " GeV].";
       }
     }
     /* If the above attempts failed, try a scan in a larger interval with a fine
      * resolution. This is a compromise between trying to be sure to find a root
      * and not introduce a too large overhead. As plan B it should be fine. */
     constexpr double half_interval_width = 100;  // GeV
     constexpr double scanning_resolution = 0.1;  // GeV
     const double upper_bound = initial_guess + half_interval_width;
     const double lower_bound = initial_guess - half_interval_width;
     std::pair<double, double> x_range = {initial_guess - scanning_resolution,
                                          initial_guess + scanning_resolution};
     logg[LPotentials].debug()
         << "Trying to find a root for potentials around " << initial_guess
         << " GeV in the range [" << lower_bound << ", " << upper_bound
         << "] GeV\nExploring interval with a resolution of "
         << scanning_resolution << " GeV, starting with x_range: ["
         << x_range.first << ", " << x_range.second << "] GeV";
     while (true) {
       const auto calc_frame_energy = root_solver.try_find_root(
           x_range.first, x_range.second, max_number_root_solver_iterations);
       if (calc_frame_energy) {
         return *calc_frame_energy;
       } else {
         const bool is_possible_to_go_right = (x_range.second < upper_bound);
         const bool is_possible_to_go_left = (x_range.first > lower_bound);
         // The following logic might be compacted, but it would be much harder
         // to read and follow, hence accept here a trivial duplication of code
         if (!is_possible_to_go_left && !is_possible_to_go_right) {
           break;
         } else if (!is_possible_to_go_left) {
           x_range.second += scanning_resolution;
         } else if (!is_possible_to_go_right) {
           x_range.first -= scanning_resolution;
         } else {
           const bool go_right = random::uniform_int(0, 1);
           if (go_right) {
             x_range.second += scanning_resolution;
           } else {
             x_range.first -= scanning_resolution;
           }
         }
         logg[LPotentials].trace() << "x_range: [" << x_range.first << ", "
                                   << x_range.second << "] GeV";
       }
     }
     logg[LPotentials].debug()
         << "Did not find any sub-range with a root scanning the interval ["
         << x_range.first << ", " << x_range.second << "] GeV";
     logg[LPotentials].error(
         "Failed to find root for momentum-dependent potentials.");
     throw std::runtime_error("Unable to continue simulation.");
   }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ skyrme_pot() [1/2]

double smash::Potentials::skyrme_pot ( const double baryon_density ) const

inline

Evaluates Skyrme potential given a baryon density.

Parameters

[in] baryon_density Baryon density \(\rho\) evaluated in the local rest frame in fm \(^{-3}\).

Returns: Skyrme potential
\[U_B=10^{-3}\times\frac{\rho}{|\rho|} (A\frac{\rho}{\rho_0}+B(\frac{\rho}{\rho_0})^\tau)\]
in GeV

Definition at line 202 of file potentials.h.

                                                        {
     return skyrme_pot(baryon_density, skyrme_a_, skyrme_b_, skyrme_tau_);
   }

Here is the caller graph for this function:

◆ symmetry_pot()

double smash::Potentials::symmetry_pot	(	const double	baryon_isospin_density,
		const double	baryon_density
	)		const

Evaluates symmetry potential given baryon isospin density.

Note: The second term is neglected if \(\gamma\) is not specified in the config.

Parameters

[in]	baryon_isospin_density	The difference between the proton and the neutron density in the local rest frame in fm \(^{-3}\).
[in]	baryon_density

Returns: Symmetry potential
\[U_I=2\times 10^{-3}S_{\rm sym} \frac{\rho_{I_3}}{\rho_0} + \left[12.3\left(\frac{\rho_B}{\rho_0}\right)^{2/3} + 20\left(\frac{\rho_B}{\rho_0}\right)^\gamma\right] \left(\frac{\rho_{I_3}}{\rho_B}\right)^2\]
in GeV

Definition at line 86 of file potentials.cc.

                                                                    {
   double pot = mev_to_gev * 2. * symmetry_S_Pot_ * baryon_isospin_density /
                nuclear_density;
   if (symmetry_is_rhoB_dependent_) {
     pot += mev_to_gev * symmetry_S(baryon_density) * baryon_isospin_density *
            baryon_isospin_density / (baryon_density * baryon_density);
   }
   return pot;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ symmetry_S()

double smash::Potentials::symmetry_S ( const double baryon_density ) const

Calculate the factor \(S(\rho)\) in the symmetry potential.

Parameters

[in] baryon_density baryon density

Returns: factor S in symmetry potenial

Definition at line 78 of file potentials.cc.

                                                                {
   if (symmetry_is_rhoB_dependent_) {
     return 12.3 * std::pow(baryon_density / nuclear_density, 2. / 3.) +
            20.0 * std::pow(baryon_density / nuclear_density, symmetry_gamma_);
   } else {
     return 0.;
   }
 }

Here is the caller graph for this function:

◆ vdf_pot()

FourVector smash::Potentials::vdf_pot	(	double	rhoB,
		const FourVector	jmuB_net
	)		const

Evaluates the FourVector potential in the VDF model given the rest frame density and the computational frame baryon current.

Parameters

[in]	rhoB	rest frame baryon density, in fm \(^{-3}\)
[in]	jmuB_net	net baryon current in the computational frame, in fm \(^{-3}\)

Returns: VDF potential
\[A^{\mu} = 10^{-3}\times \sum_i C_i \left(\frac{\rho}{\rho_0}\right)^{b_i - 2} \frac{j^{\mu}}{\rho_0}\]
in GeV

Definition at line 97 of file potentials.cc.

                                                                            {
   // this needs to be used in order to prevent trying to calculate something
   // like
   // (-rho_B)^{3.4}
   const int sgn = rhoB > 0 ? 1 : -1;
   double abs_rhoB = std::abs(rhoB);
   // to prevent NAN expressions
   if (abs_rhoB < very_small_double) {
     abs_rhoB = very_small_double;
   }
   // F_2 is a multiplicative factor in front of the baryon current
   // in the VDF potential
   double F_2 = 0.0;
   for (int i = 0; i < number_of_terms(); i++) {
     F_2 += coeffs_[i] * std::pow(abs_rhoB, powers_[i] - 2.0) /
            std::pow(saturation_density_, powers_[i] - 1.0);
   }
   F_2 = F_2 * sgn;
   // Return in GeV
   return F_2 * jmuB_net;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ potential()

double smash::Potentials::potential	(	const ThreeVector &	r,
		const ParticleList &	plist,
		const ParticleType &	acts_on
	)		const

Evaluates potential (Skyrme with optional Symmetry or VDF) at point r.

For Skyrme and Symmetry options, potential is always taken in the local Eckart rest frame, but point r is in the computational frame.

Parameters

[in]	r	Arbitrary space point where potential is calculated
[in]	plist	List of all particles to be used in \(j^{\mu}\) calculation. If the distance between particle and calculation point r, \( \|r-r_i\| > r_{cut} \) then particle input to density will be ignored.
[in]	acts_on	Type of particle on which potential is going to act. It gives the charges (or more precisely, the scaling factors) of the particle moving in the potential field.

Returns

Total potential energy acting on the particle: for Skyrme and Symmetry potentials,

\[U_{\rm tot} =Q_BU_B+2I_3U_I\]

in GeV, while for the VDF potential

\[U_{\rm tot} =Q_B A^0\]

in GeV, where \(Q_B\) is the baryon charge scaled by the ratio of the light (u, d) quark to the total quark number and \(I_3\) is the third compnent of the isospin.

Definition at line 119 of file potentials.cc.

                                                                 {
   double total_potential = 0.0;
   const bool compute_gradient = false;
   const bool smearing = true;
   const auto scale = force_scale(acts_on);
  
   if (!(acts_on.is_baryon() || acts_on.is_nucleus())) {
     return total_potential;
   }
   const auto baryon_density_and_gradient = current_eckart(
       r, plist, param_, DensityType::Baryon, compute_gradient, smearing);
   const double rhoB = std::get<0>(baryon_density_and_gradient);
   if (use_skyrme_) {
     total_potential += scale.first * skyrme_pot(rhoB);
   }
   if (use_symmetry_) {
     const double rho_iso = std::get<0>(
         current_eckart(r, plist, param_, DensityType::BaryonicIsospin,
                        compute_gradient, smearing));
     const double sym_pot = symmetry_pot(rho_iso, rhoB) * acts_on.isospin3_rel();
     total_potential += scale.second * sym_pot;
   }
  
   if (use_vdf_) {
     const FourVector jmuB = std::get<1>(baryon_density_and_gradient);
     const FourVector VDF_potential = vdf_pot(rhoB, jmuB);
     total_potential += scale.first * VDF_potential.x0();
   }
  
   return total_potential;
 }

Here is the call graph for this function:

◆ force_scale()

std::pair< double, int > smash::Potentials::force_scale ( const ParticleType & data )

static

Evaluates the scaling factor of the forces acting on the particles.

The forces are equal to the product of the scaling factor and the gradient of the potential. We need these scaling factors to describe the motions of the hyperons as well as the anti-particles in the potentials. For Lambda and Sigma, since they carry 2 light (u or d) quarks, they are affected by 2/3 of the Skyrme force. Xi carries 1 light quark, it is affected by 1/3 of the Skyrme force. Omega carries no light quark, so it's not affected by the Skyrme force. Anti-baryons are affected by the force as large as the force acting on baryons but with an opposite direction.

Parameters

[in] data Type of particle on which potential is going to act.

Returns: ( \(Q_B(1-\frac{|Q_S|}{3}), Q_B\)) where \(Q_B\) is the baryon charge and \(Q_S\) is the strangeness.

Definition at line 152 of file potentials.cc.

                                                                      {
   const auto &pdg = data.pdgcode();
   const double skyrme_or_VDF_scale =
       (3 - std::abs(pdg.strangeness())) / 3. * pdg.baryon_number();
   const int symmetry_scale = pdg.baryon_number();
   return std::make_pair(skyrme_or_VDF_scale, symmetry_scale);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ skyrme_force()

std::pair< ThreeVector, ThreeVector > smash::Potentials::skyrme_force	(	const double	rhoB,
		const ThreeVector	grad_j0B,
		const ThreeVector	dvecjB_dt,
		const ThreeVector	curl_vecjB
	)		const

Evaluates the electric and magnetic components of the skyrme force.

Parameters

[in]	rhoB	Eckart baryon density [fm \(^{-3}\)].
[in]	grad_j0B	Gradient of baryon density [fm \(^{-4}\)]. This density is evaluated in the computational frame.
[in]	dvecjB_dt	Time derivative of the vector baryon current density [fm \(^{-4}\)
[in]	curl_vecjB	Curl of the baryon vector current density [fm \(^{-4}\)

Returns

( \(E_B, B_B\)), where

\[ E_B = - V_B^\prime(\rho^\ast)(\boldsymbol{\nabla}\rho_B + \partial_t\,\mathbf{j}_B) \]

is the electro component of Skyrme force and

\[ B_B = V_B^\prime(\rho^\ast) \boldsymbol{\nabla}\times\mathbf{j}_B \]

is the magnetic component of the Skyrme force with \(\rho^\ast\) being the Eckart baryon density.

Definition at line 160 of file potentials.cc.

                                         {
   ThreeVector E_component(0.0, 0.0, 0.0), B_component(0.0, 0.0, 0.0);
   if (use_skyrme_) {
     const int sgn = rhoB > 0 ? 1 : -1;
     const double abs_rhoB = std::abs(rhoB);
     const double dV_drho = sgn *
                            (skyrme_a_ + skyrme_b_ * skyrme_tau_ *
                                             std::pow(abs_rhoB / nuclear_density,
                                                      skyrme_tau_ - 1)) *
                            mev_to_gev / nuclear_density;
     E_component -= dV_drho * (grad_j0B + dvecjB_dt);
     B_component += dV_drho * curl_vecjB;
   }
   return std::make_pair(E_component, B_component);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ symmetry_force()

std::pair< ThreeVector, ThreeVector > smash::Potentials::symmetry_force	(	const double	rhoI3,
		const ThreeVector	grad_j0I3,
		const ThreeVector	dvecjI3_dt,
		const ThreeVector	curl_vecjI3,
		const double	rhoB,
		const ThreeVector	grad_j0B,
		const ThreeVector	dvecjB_dt,
		const ThreeVector	curl_vecjB
	)		const

Evaluates the electric and magnetic components of the symmetry force.

Parameters

[in]	rhoI3	Relative isospin 3 density.
[in]	grad_j0I3	Gradient of I3/I density [fm \(^{-4}\)]. This density is evaluated in the computational frame.
[in]	dvecjI3_dt	Time derivative of the I3/I vector current density [fm \(^{-4}\)]
[in]	curl_vecjI3	Curl of the I3/I vector current density [fm \(^{-4}\)]
[in]	rhoB	Net-baryon density in the rest frame
[in]	grad_j0B	Gradient of the net-baryon density in the computational frame
[in]	dvecjB_dt	Time derivative of the net-baryon vector current density
[in]	curl_vecjB	Curl of the net-baryon vector current density

Returns

( \(E_{I_3}, B_{I_3}\)) [GeV/fm], where

\[ \mathbf{E} = - \frac{\partial V^\ast}{\partial\rho_{I_3}^\ast} (\boldsymbol{\nabla}\rho_{I_3} + \partial_t \mathbf{j}_{I_3}) - \frac{\partial V^\ast}{\partial\rho_B^\ast} (\boldsymbol{\nabla}\rho_B + \partial_t\,\mathbf{j}_B) \]

is the electrical component of symmetry force and

\[ \mathbf{B} = \frac{\partial V^\ast}{\rho_{I_3}^\ast} \boldsymbol{\nabla}\times\mathbf{j}_{I_3} + \frac{\partial V^\ast}{\rho_B^\ast} \boldsymbol{\nabla}\times\mathbf{j}_B \]

is the magnetic component of the symmetry force with \(\rho^\ast\) being the respective Eckart density.

Definition at line 178 of file potentials.cc.

                                         {
   ThreeVector E_component(0.0, 0.0, 0.0), B_component(0.0, 0.0, 0.0);
   if (use_symmetry_) {
     E_component -= dVsym_drhoI3(rhoB, rhoI3) * (grad_j0I3 + dvecjI3_dt) +
                    dVsym_drhoB(rhoB, rhoI3) * (grad_j0B + dvecjB_dt);
     B_component += dVsym_drhoI3(rhoB, rhoI3) * curl_vecjI3 +
                    dVsym_drhoB(rhoB, rhoI3) * curl_vecjB;
   }
   return std::make_pair(E_component, B_component);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ E_field_integrand()

static ThreeVector smash::Potentials::E_field_integrand	(	ThreeVector	pos,
		DensityOnLattice &	charge_density,
		ThreeVector	point
	)

inlinestatic

Integrand for calculating the electric field.

The field is calculated via

\[ \mathbf{E}(\mathbf{r}) = \int\frac{(\mathbf{r}-\mathbf{r}^\prime)\rho(\mathbf{r}^\prime)} {|\mathbf{r}-\mathbf{r}^\prime|^3}d^3r^\prime \]

Parameters

[in]	pos	position vector to be integrated over
[in]	charge_density	electric charge density at position pos
[in]	point	position where to calculate the field

Definition at line 365 of file potentials.h.

                                                           {
     ThreeVector dr = point - pos;
     if (dr.abs() < really_small) {
       return {0., 0., 0.};
     }
     return elementary_charge * charge_density.rho() * dr /
            std::pow(dr.abs(), 3);
   }

Here is the call graph for this function:

◆ B_field_integrand()

static ThreeVector smash::Potentials::B_field_integrand	(	ThreeVector	pos,
		DensityOnLattice &	charge_density,
		ThreeVector	point
	)

inlinestatic

Integrand for calculating the magnetic field using the Biot-Savart formula.

Parameters

[in]	pos	position vector to be integrated over
[in]	charge_density	electric charge density and current
[in]	point	position where the magnetic field will be calculated

Definition at line 383 of file potentials.h.

                                                           {
     ThreeVector dr = point - pos;
     if (dr.abs() < really_small) {
       return {0., 0., 0.};
     }
     return elementary_charge *
            charge_density.jmu_net().threevec().cross_product(dr) /
            std::pow(dr.abs(), 3);
   }

Here is the call graph for this function:

◆ vdf_force() [1/2]

std::pair< ThreeVector, ThreeVector > smash::Potentials::vdf_force	(	double	rhoB,
		const double	drhoB_dt,
		const ThreeVector	grad_rhoB,
		const ThreeVector	gradrhoB_cross_vecjB,
		const double	j0B,
		const ThreeVector	grad_j0B,
		const ThreeVector	vecjB,
		const ThreeVector	dvecjB_dt,
		const ThreeVector	curl_vecjB
	)		const

Evaluates the electric and magnetic components of force in the VDF model given the derivatives of the baryon current \(j^{\mu}\).

Parameters

[in]	rhoB	rest frame baryon density in fm \(^{-3}\)
[in]	drhoB_dt	time derivative of the rest frame density
[in]	grad_rhoB	gradient of the rest frame density
[in]	gradrhoB_cross_vecjB	cross product of the gradient of the rest frame density and the 3-vector baryon current density
[in]	j0B	computational frame baryon density in fm \(^{-3}\)
[in]	grad_j0B	gradient of the computational frame baryon density
[in]	vecjB	3-vector baryon current
[in]	dvecjB_dt	time derivative of the computational frame 3-vector baryon current
[in]	curl_vecjB	curl of the 3-vector baryon current

Returns

( \(E_{VDF}, B_{VDF}\)) [GeV/fm], where

\[ \mathbf{E}_{VDF} = - F_1 \big[(\boldsymbol{\nabla} \rho)\,j^0 + (\partial_t \rho)\,\mathbf{j}\big] - F_2 (\boldsymbol{\nabla} j^0 + \partial_t\,\mathbf{j}) \]

is the electrical component of VDF force and

\[ \mathbf{B}_{VDF} = F_1 (\boldsymbol{\nabla} \rho) \times \mathbf{j} + F_2 \boldsymbol{\nabla} \times \mathbf{j} \]

is the magnetic component of the VDF force, with

\begin{aligned} F_1 &= \sum_i C_i (b_i - 2) \frac{\rho^{b_i - 3}}{\rho_0^{b_i - 1}}\\ F_2 &= \sum_i C_i \frac{\rho^{b_i - 2}}{\rho_0^{b_i - 1}} \;, \end{aligned}

where \(\rho_0\) is the saturation density.

Definition at line 193 of file potentials.cc.

                                                                      {
   // this needs to be used to prevent trying to calculate something like
   // (-rhoB)^{3.4}
   const int sgn = rhoB > 0 ? 1 : -1;
   ThreeVector E_component(0.0, 0.0, 0.0), B_component(0.0, 0.0, 0.0);
   // to prevent NAN expressions
   double abs_rhoB = std::abs(rhoB);
   if (abs_rhoB < very_small_double) {
     abs_rhoB = very_small_double;
   }
   if (use_vdf_) {
     // F_1 and F_2 are multiplicative factors in front of the baryon current
     // in the VDF potential
     double F_1 = 0.0;
     for (int i = 0; i < number_of_terms(); i++) {
       F_1 += coeffs_[i] * (powers_[i] - 2.0) *
              std::pow(abs_rhoB, powers_[i] - 3.0) /
              std::pow(saturation_density_, powers_[i] - 1.0);
     }
     F_1 = F_1 * sgn;
  
     double F_2 = 0.0;
     for (int i = 0; i < number_of_terms(); i++) {
       F_2 += coeffs_[i] * std::pow(abs_rhoB, powers_[i] - 2.0) /
              std::pow(saturation_density_, powers_[i] - 1.0);
     }
     F_2 = F_2 * sgn;
  
     E_component -= (F_1 * (grad_rhoB * j0B + drhoB_dt * vecjB) +
                     F_2 * (grad_j0B + dvecjB_dt));
     B_component += F_1 * gradrhoB_cross_vecjB + F_2 * curl_vecjB;
   }
   return std::make_pair(E_component, B_component);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ vdf_force() [2/2]

std::pair< ThreeVector, ThreeVector > smash::Potentials::vdf_force	(	const ThreeVector	grad_A_0,
		const ThreeVector	dA_dt,
		const ThreeVector	curl_vecA
	)		const

Evaluates the electric and magnetic components of force in the VDF force given the derivatives of the VDF mean-field \(A^\mu\).

Parameters

[in]	grad_A_0	gradient of the zeroth component of the field A^mu
[in]	dA_dt	time derivative of the field A^mu
[in]	curl_vecA	curl of the vector component of the field A^mu

Returns

( \(E_{VDF}, B_{VDF}\)) [GeV/fm], where

\[ \mathbf{E}_{VDF} = - \boldsymbol{\nabla} A^0 - \partial_t\mathbf{A} \]

is the electrical component of VDF force and

\[ \mathbf{B}_{VDF} = \boldsymbol{\nabla} \times \mathbf{A} \]

is the magnetic component of the VDF force.

Definition at line 233 of file potentials.cc.

                                     {
   ThreeVector E_component(0.0, 0.0, 0.0), B_component(0.0, 0.0, 0.0);
   if (use_vdf_) {
     E_component -= (grad_A_0 + dA_dt);
     B_component += curl_A;
   }
   return std::make_pair(E_component, B_component);
 }

◆ all_forces()

std::tuple< ThreeVector, ThreeVector, ThreeVector, ThreeVector > smash::Potentials::all_forces	(	const ThreeVector &	r,
		const ParticleList &	plist
	)		const

virtual

Evaluates the electric and magnetic components of the forces at point r.

Point r is in the computational frame.

Parameters

[in]	r	Arbitrary space point where potential gradient is calculated
[in]	plist	List of all particles to be used in \(j^{\mu}\) calculation. If the distance between particle and calculation point r, \( \|r-r_i\| > r_{cut} \) then particle input to density will be ignored.

Returns: ( \(E_B, B_B, E_{I_3}, B_{I_3}\)) [GeV/fm], where \(E_B\): the electric component of the Skyrme or VDF force, \(B_B\): the magnetic component of the Skyrme or VDF force, \(E_{I_3}\): the electric component of the symmetry force, \(B_{I_3}\): the magnetic component of the symmetry force

Definition at line 268 of file potentials.cc.

                                                                             {
   const bool compute_gradient = true;
   const bool smearing = true;
   auto F_skyrme_or_VDF =
       std::make_pair(ThreeVector(0., 0., 0.), ThreeVector(0., 0., 0.));
   auto F_symmetry =
       std::make_pair(ThreeVector(0., 0., 0.), ThreeVector(0., 0., 0.));
  
   const auto baryon_density_and_gradient = current_eckart(
       r, plist, param_, DensityType::Baryon, compute_gradient, smearing);
   double rhoB = std::get<0>(baryon_density_and_gradient);
   const ThreeVector grad_j0B = std::get<2>(baryon_density_and_gradient);
   const ThreeVector curl_vecjB = std::get<3>(baryon_density_and_gradient);
   const FourVector djmuB_dt = std::get<4>(baryon_density_and_gradient);
   if (use_skyrme_) {
     F_skyrme_or_VDF =
         skyrme_force(rhoB, grad_j0B, djmuB_dt.threevec(), curl_vecjB);
   }
  
   if (use_symmetry_) {
     const auto density_and_gradient =
         current_eckart(r, plist, param_, DensityType::BaryonicIsospin,
                        compute_gradient, smearing);
     const double rhoI3 = std::get<0>(density_and_gradient);
     const ThreeVector grad_j0I3 = std::get<2>(density_and_gradient);
     const ThreeVector curl_vecjI3 = std::get<3>(density_and_gradient);
     const FourVector dvecjI3_dt = std::get<4>(density_and_gradient);
     F_symmetry =
         symmetry_force(rhoI3, grad_j0I3, dvecjI3_dt.threevec(), curl_vecjI3,
                        rhoB, grad_j0B, djmuB_dt.threevec(), curl_vecjB);
   }
  
   if (use_vdf_) {
     const FourVector jmuB = std::get<1>(baryon_density_and_gradient);
     const FourVector djmuB_dx = std::get<5>(baryon_density_and_gradient);
     const FourVector djmuB_dy = std::get<6>(baryon_density_and_gradient);
     const FourVector djmuB_dz = std::get<7>(baryon_density_and_gradient);
  
     // safety check to not divide by zero
     const int sgn = rhoB > 0 ? 1 : -1;
     if (std::abs(rhoB) < very_small_double) {
       rhoB = sgn * very_small_double;
     }
  
     const double drhoB_dt =
         (1 / rhoB) * (jmuB.x0() * djmuB_dt.x0() - jmuB.x1() * djmuB_dt.x1() -
                       jmuB.x2() * djmuB_dt.x2() - jmuB.x3() * djmuB_dt.x3());
  
     const double drhoB_dx =
         (1 / rhoB) * (jmuB.x0() * djmuB_dx.x0() - jmuB.x1() * djmuB_dx.x1() -
                       jmuB.x2() * djmuB_dx.x2() - jmuB.x3() * djmuB_dx.x3());
  
     const double drhoB_dy =
         (1 / rhoB) * (jmuB.x0() * djmuB_dy.x0() - jmuB.x1() * djmuB_dy.x1() -
                       jmuB.x2() * djmuB_dy.x2() - jmuB.x3() * djmuB_dy.x3());
  
     const double drhoB_dz =
         (1 / rhoB) * (jmuB.x0() * djmuB_dz.x0() - jmuB.x1() * djmuB_dz.x1() -
                       jmuB.x2() * djmuB_dz.x2() - jmuB.x3() * djmuB_dz.x3());
  
     const FourVector drhoB_dxnu = {drhoB_dt, drhoB_dx, drhoB_dy, drhoB_dz};
  
     const ThreeVector grad_rhoB = drhoB_dxnu.threevec();
     const ThreeVector vecjB = jmuB.threevec();
     const ThreeVector Drho_cross_vecj = grad_rhoB.cross_product(vecjB);
  
     F_skyrme_or_VDF = vdf_force(
         rhoB, drhoB_dt, drhoB_dxnu.threevec(), Drho_cross_vecj, jmuB.x0(),
         grad_j0B, jmuB.threevec(), djmuB_dt.threevec(), curl_vecjB);
   }
  
   return std::make_tuple(F_skyrme_or_VDF.first, F_skyrme_or_VDF.second,
                          F_symmetry.first, F_symmetry.second);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ use_skyrme()

virtual bool smash::Potentials::use_skyrme ( ) const

inlinevirtual

Returns: Is Skyrme potential on?

Definition at line 478 of file potentials.h.

478 { return use_skyrme_; }

Here is the caller graph for this function:

◆ use_symmetry()

virtual bool smash::Potentials::use_symmetry ( ) const

inlinevirtual

Returns: Is symmetry potential on?

Definition at line 480 of file potentials.h.

480 { return use_symmetry_; }

Here is the caller graph for this function:

◆ use_coulomb()

virtual bool smash::Potentials::use_coulomb ( ) const

inlinevirtual

Returns: Is Coulomb potential on?

Definition at line 482 of file potentials.h.

482 { return use_coulomb_; }

Here is the caller graph for this function:

◆ use_momentum_dependence()

virtual bool smash::Potentials::use_momentum_dependence ( ) const

inlinevirtual

Returns: Use momentum-dependent part of the potential?

Definition at line 484 of file potentials.h.

                                                {
     return use_momentum_dependence_;
   }

Here is the caller graph for this function:

◆ skyrme_a()

double smash::Potentials::skyrme_a ( ) const

inline

Returns: Skyrme parameter skyrme_a, in MeV

Definition at line 489 of file potentials.h.

489 { return skyrme_a_; }

◆ skyrme_b()

double smash::Potentials::skyrme_b ( ) const

inline

Returns: Skyrme parameter skyrme_b, in MeV

Definition at line 491 of file potentials.h.

491 { return skyrme_b_; }

◆ skyrme_tau()

double smash::Potentials::skyrme_tau ( ) const

inline

Returns: Skyrme parameter skyrme_tau

Definition at line 493 of file potentials.h.

493 { return skyrme_tau_; }

◆ symmetry_S_pot()

double smash::Potentials::symmetry_S_pot ( ) const

inline

Returns: Skyrme parameter S_pot, in MeV

Definition at line 495 of file potentials.h.

495 { return symmetry_S_Pot_; }

◆ use_vdf()

virtual bool smash::Potentials::use_vdf ( ) const

inlinevirtual

Returns: Is VDF potential on?

Definition at line 498 of file potentials.h.

498 { return use_vdf_; }

Here is the caller graph for this function:

◆ saturation_density()

double smash::Potentials::saturation_density ( ) const

inline

Returns: Value of the saturation density used in the VDF potential

Definition at line 500 of file potentials.h.

500 { return saturation_density_; }

◆ coeffs()

const std::vector<double>& smash::Potentials::coeffs ( ) const

inline

Returns: Vector of the VDF coefficients \(C_i\), coefficients_

Definition at line 502 of file potentials.h.

502 { return coeffs_; }

◆ powers()

const std::vector<double>& smash::Potentials::powers ( ) const

inline

Returns: Vector of the VDF exponents \(b_i\), powers_

Definition at line 504 of file potentials.h.

504 { return powers_; }

◆ number_of_terms()

int smash::Potentials::number_of_terms ( ) const

inline

Returns: Number of terms in the VDF potential

Definition at line 506 of file potentials.h.

506 { return powers_.size(); }

Here is the caller graph for this function:

◆ coulomb_r_cut()

double smash::Potentials::coulomb_r_cut ( ) const

inline

Returns: cutoff radius in ntegration for coulomb potential in fm

Definition at line 509 of file potentials.h.

509 { return coulomb_r_cut_; }

◆ use_potentials_outside_lattice()

bool smash::Potentials::use_potentials_outside_lattice ( ) const

inline

Returns: Wether to take potentials into account for particles outside of the lattice

Definition at line 515 of file potentials.h.

                                               {
     return use_potentials_outside_lattice_;
   }

Here is the caller graph for this function:

◆ dVsym_drhoI3()

double smash::Potentials::dVsym_drhoI3	(	const double	rhoB,
		const double	rhoI3
	)		const

private

Calculate the derivative of the symmetry potential with respect to the isospin density in GeV * fm^3.

\[ \frac{\partial V_\mathrm{sym}}{\partial \rho_{I_3}} = 2\frac{S_\mathrm{Pot}}{\rho_0} + \frac{2\rho_{I_3}\left[12.3\left(\frac{\rho_B}{\rho_0}\right)^{2/3} + 20 \left(\frac{\rho_B}{\rho_0}\right)^\gamma\right]}{\rho_B^2} \]

Note: The isospin 3 density here is actually the density of I3 / I.

Parameters

[in]	rhoB	net baryon density
[in]	rhoI3	isospin density

Returns: partial derivative of the symmetry potenital with respect to the isospin density.

Definition at line 244 of file potentials.cc.

                                                                            {
   double term1 = 2. * symmetry_S_Pot_ / nuclear_density;
   if (symmetry_is_rhoB_dependent_) {
     double term2 = 2. * rhoI3 * symmetry_S(rhoB) / (rhoB * rhoB);
     return mev_to_gev * (term1 + term2);
   } else {
     return mev_to_gev * term1;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ dVsym_drhoB()

double smash::Potentials::dVsym_drhoB	(	const double	rhoB,
		const double	rhoI3
	)		const

private

Calculate the derivative of the symmetry potential with respect to the net baryon density in GeV * fm^3.

\[ \frac{\partial V_\mathrm{sym}}{\partial \rho_B} = \left(\frac{\rho_{I_3}}{\rho_B}\right)^2 \left[\frac{8.2}{\rho_0}\left(\frac{\rho_B}{\rho_0}\right)^{-1/3} + \frac{20\gamma}{\rho_B}\left(\frac{\rho_B}{\rho_0}\right)^\gamma\right] -2\frac{\rho_{I_3}^2}{\rho_B^3} \left[12.3\left(\frac{\rho_B}{\rho_0}\right)^{2/3} + 20\left(\frac{\rho_B}{\rho_0}\right)^\gamma\right]\]

Note: The isospin 3 density here is actually the density of I3 / I

Parameters

[in]	rhoB	net baryon density
[in]	rhoI3	isospin density

Returns: partial derivative of the symmetry potenital with respect to the net baryon density.

Definition at line 254 of file potentials.cc.

                                                                           {
   if (symmetry_is_rhoB_dependent_) {
     double rhoB_over_rho0 = rhoB / nuclear_density;
     double term1 = 8.2 * std::pow(rhoB_over_rho0, -1. / 3.) / nuclear_density +
                    20. * symmetry_gamma_ *
                        std::pow(rhoB_over_rho0, symmetry_gamma_) / rhoB;
     double term2 = -2. * symmetry_S(rhoB) / rhoB;
     return mev_to_gev * (term1 + term2) * rhoI3 * rhoI3 / (rhoB * rhoB);
   } else {
     return 0.;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ skyrme_pot() [2/2]

double smash::Potentials::skyrme_pot	(	const double	baryon_density,
		const double	A,
		const double	B,
		const double	tau
	)

staticprivate

Single particle Skyrme potential in MeV.

Parameters

baryon_density	net baryon density in the local rest-frame in 1/fm^3
A	Skyrme parameter A in MeV
B	Skyrme parameter B in MeV
tau	Skyrme parameter tau

Returns: Single particle Skyrme potential in MeV

Definition at line 67 of file potentials.cc.

                                                                 {
   const double tmp = baryon_density / nuclear_density;
   /* U = U(|rho|) * sgn , because the sign of the potential changes
    * under a charge reversal transformation. */
   const int sgn = tmp > 0 ? 1 : -1;
   // Return in GeV
   return mev_to_gev * sgn *
          (A * std::abs(tmp) + B * std::pow(std::abs(tmp), tau));
 }

Here is the call graph for this function:

◆ root_eq_potentials()

static double smash::Potentials::root_eq_potentials	(	double	energy_calc,
		const ThreeVector &	momentum_calc,
		const FourVector &	jmu,
		double	m,
		double	A,
		double	B,
		double	tau,
		double	C,
		double	Lambda
	)

inlinestaticprivate

Root equation used to determine the energy in the calculation frame.

It is the difference between energy (including potential) squared minus momentum squared in the calculation frame and in the local rest-frame. The equation should be zero due to Lorentz invariance but a root finder is required to determine the calculation frame energy such that this is indeed the case.

Parameters

[in]	energy_calc	Energy in the calculation frame at which the equation should be evaluated in GeV
[in]	momentum_calc	Momentum of the particle in calculation frame of interest in GeV
[in]	jmu	Baryon current fourvector at the position of the particle
[in]	m	Mass of the particle of interest in GeV
[in]	A	Skyrme parameter A in MeV
[in]	B	Skyrme parameter B in MeV
[in]	tau	Skyrme parameter tau
[in]	C	Parameter C of the momentum dependent part of the potential in MeV
[in]	Lambda	Parameter Lambda of the momentum-dependent term of the potential in 1/fm

Returns: effective mass squared in calculation frame minus effective mass in rest_frame in GeV^2

Definition at line 679 of file potentials.h.

                                                   {
     // get velocity for boost to the local rest frame
     double rho_LRF = jmu.abs();
     ThreeVector beta_LRF = jmu.x0() > really_small ? jmu.threevec() / jmu.x0()
                                                    : ThreeVector(0, 0, 0);
     // get momentum in the local rest frame
     FourVector pmu_calc = FourVector(energy_calc, momentum_calc);
     FourVector pmu_LRF = beta_LRF.abs() > really_small
                              ? pmu_calc.lorentz_boost(beta_LRF)
                              : pmu_calc;
     double p_LRF = pmu_LRF.threevec().abs();
     double energy_LRF = std::sqrt(m * m + p_LRF * p_LRF) +
                         skyrme_pot(rho_LRF, A, B, tau) +
                         momentum_dependent_part(p_LRF, rho_LRF, C, Lambda);
     const double result = energy_calc * energy_calc - momentum_calc.sqr() -
                           (energy_LRF * energy_LRF - p_LRF * p_LRF);
     logg[LPotentials].debug()
         << "root equation for potentials called with E_calc=" << energy_calc
         << " p_calc=" << momentum_calc << " jmu=" << jmu << " m=" << m
         << " tau=" << tau << " A=" << A << " B=" << B << " C=" << C
         << " Lambda=" << Lambda << " and the root equation is " << result;
     return result;
   }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ momentum_dependent_part()

static double smash::Potentials::momentum_dependent_part	(	double	momentum,
		double	rho,
		double	C,
		double	Lambda
	)

inlinestaticprivate

Momentum dependent term of the potential.

To be added to the momentum independent part.

Parameters

[in]	momentum	Absolute momentum of the particle of interest in the local rest-frame in GeV
[in]	rho	Baryon density in the Eckart frame in 1/fm^3
[in]	C	Parameter C of the momentum dependent part of the potential in MeV
[in]	Lambda	Parameter Lambda of the momentum-dependent part of the potential in 1/fm

Returns: momentum dependent part of the potential in GeV

Definition at line 721 of file potentials.h.

                                                        {
     /* We assume here that the distribution function is the one of cold
      * nuclear matter, which consists only of protons and neutrons.
      * That is, the degeneracy factor simply equals nucleon spin degeneracy
      * times nucleon isospin degeneracy, even though all baryons contribute
      * to the density and the potential is applied to all baryons. */
     int g = 4;
     const double fermi_momentum =
         std::cbrt(6. * M_PI * M_PI * rho / g);  // in 1/fm
     momentum = momentum / hbarc;                // convert to 1/fm
     if (unlikely(momentum < really_small)) {
       return mev_to_gev * g * C / (M_PI * M_PI * nuclear_density) *
              (Lambda * Lambda * fermi_momentum -
               std::pow(Lambda, 3) * std::atan(fermi_momentum / Lambda));
     }
     const std::array<double, 7> temp = {
         2 * g * C * M_PI * std::pow(Lambda, 3) /
             (std::pow(2 * M_PI, 3) * nuclear_density),
         (fermi_momentum * fermi_momentum + Lambda * Lambda -
          momentum * momentum) /
             (2 * momentum * Lambda),
         std::pow(momentum + fermi_momentum, 2) + Lambda * Lambda,
         std::pow(momentum - fermi_momentum, 2) + Lambda * Lambda,
         2 * fermi_momentum / Lambda,
         (momentum + fermi_momentum) / Lambda,
         (momentum - fermi_momentum) / Lambda};
     const double result =
         temp[0] * (temp[1] * std::log(temp[2] / temp[3]) + temp[4] -
                    2 * (std::atan(temp[5]) - std::atan(temp[6])));
     return mev_to_gev * result;
   }

Here is the caller graph for this function:

Member Data Documentation

◆ param_

const DensityParameters smash::Potentials::param_

private

Struct that contains the gaussian smearing width \(\sigma\), the distance cutoff \(r_{\rm cut}\) and the testparticle number needed for the density calculation.

Definition at line 525 of file potentials.h.

◆ use_skyrme_

bool smash::Potentials::use_skyrme_

private

Skyrme potential on/off.

Definition at line 528 of file potentials.h.

◆ use_symmetry_

bool smash::Potentials::use_symmetry_

private

Symmetry potential on/off.

Definition at line 531 of file potentials.h.

◆ use_coulomb_

bool smash::Potentials::use_coulomb_

private

Coulomb potential on/Off.

Definition at line 534 of file potentials.h.

◆ use_vdf_

bool smash::Potentials::use_vdf_

private

VDF potential on/off.

Definition at line 537 of file potentials.h.

◆ use_momentum_dependence_

bool smash::Potentials::use_momentum_dependence_

private

Momentum-dependent part on/off.

Definition at line 540 of file potentials.h.

◆ skyrme_a_

double smash::Potentials::skyrme_a_

private

Parameter of skyrme potentials: the coefficient in front of \(\frac{\rho}{\rho_0}\) in GeV.

Definition at line 546 of file potentials.h.

◆ skyrme_b_

double smash::Potentials::skyrme_b_

private

Parameters of skyrme potentials: the coefficient in front of \((\frac{\rho}{\rho_0})^\tau\) in GeV.

Definition at line 552 of file potentials.h.

◆ skyrme_tau_

double smash::Potentials::skyrme_tau_

private

Parameters of skyrme potentials: the power index.

Definition at line 558 of file potentials.h.

◆ mom_dependence_Lambda_

double smash::Potentials::mom_dependence_Lambda_

private

Parameter Lambda of the momentum-dependent part of the potentials given in 1/fm.

Definition at line 564 of file potentials.h.

◆ mom_dependence_C_

double smash::Potentials::mom_dependence_C_

private

Parameter C of the momentum-dependent part of the potentials given in MeV.

Definition at line 570 of file potentials.h.

◆ symmetry_S_Pot_

double smash::Potentials::symmetry_S_Pot_

private

Parameter S_Pot in the symmetry potential in MeV.

Definition at line 573 of file potentials.h.

◆ symmetry_is_rhoB_dependent_

bool smash::Potentials::symmetry_is_rhoB_dependent_ = false

private

Whether the baryon density dependence of the symmetry potential is included.

Definition at line 579 of file potentials.h.

◆ symmetry_gamma_

double smash::Potentials::symmetry_gamma_

private

Power \( \gamma \) in formula for \( S(\rho) \):

\[ S(\rho)=12.3\,\mathrm{MeV}\times \left(\frac{\rho}{\rho_0}\right)^{2/3}+20\,\mathrm{MeV}\times \left(\frac{\rho}{\rho_0}\right)^\gamma \]

Definition at line 586 of file potentials.h.

◆ coulomb_r_cut_

double smash::Potentials::coulomb_r_cut_

private

Cutoff in integration for coulomb potential.

Definition at line 589 of file potentials.h.

◆ use_potentials_outside_lattice_

bool smash::Potentials::use_potentials_outside_lattice_

private

Wether potentials should be included outside of the lattice.

Definition at line 592 of file potentials.h.

◆ saturation_density_

double smash::Potentials::saturation_density_

private

Saturation density of nuclear matter used in the VDF potential; it may vary between different parameterizations.

Definition at line 598 of file potentials.h.

◆ coeffs_

std::vector<double> smash::Potentials::coeffs_

private

Parameters of the VDF potential: coefficients \(C_i\), in GeV.

Definition at line 600 of file potentials.h.

◆ powers_

std::vector<double> smash::Potentials::powers_

private

Parameters of the VDF potential: exponents \(b_i\).

Definition at line 602 of file potentials.h.

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

src/include/smash/potentials.h
src/potentials.cc

Public Member Functions

Static Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ Potentials()

◆ ~Potentials()

Member Function Documentation

◆ single_particle_energy_gradient()

◆ calculation_frame_energy()

◆ skyrme_pot() [1/2]

◆ symmetry_pot()

◆ symmetry_S()

◆ vdf_pot()

◆ potential()

◆ force_scale()

◆ skyrme_force()

◆ symmetry_force()

◆ E_field_integrand()

◆ B_field_integrand()

◆ vdf_force() [1/2]

◆ vdf_force() [2/2]

◆ all_forces()

◆ use_skyrme()

◆ use_symmetry()

◆ use_coulomb()

◆ use_momentum_dependence()

◆ skyrme_a()

◆ skyrme_b()

◆ skyrme_tau()

◆ symmetry_S_pot()

◆ use_vdf()

◆ saturation_density()

◆ coeffs()

◆ powers()

◆ number_of_terms()

◆ coulomb_r_cut()

◆ use_potentials_outside_lattice()

◆ dVsym_drhoI3()

◆ dVsym_drhoB()

◆ skyrme_pot() [2/2]

◆ root_eq_potentials()

◆ momentum_dependent_part()

Member Data Documentation

◆ param_

◆ use_skyrme_

◆ use_symmetry_

◆ use_coulomb_

◆ use_vdf_

◆ use_momentum_dependence_

◆ skyrme_a_

◆ skyrme_b_

◆ skyrme_tau_

◆ mom_dependence_Lambda_

◆ mom_dependence_C_

◆ symmetry_S_Pot_

◆ symmetry_is_rhoB_dependent_

◆ symmetry_gamma_

◆ coulomb_r_cut_

◆ use_potentials_outside_lattice_

◆ saturation_density_

◆ coeffs_

◆ powers_