#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 32 of file potentials.h.

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

virtual	~Potentials ()
	Standard destructor. 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

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

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...

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...

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	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...

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 coefficents 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 17 of file potentials.cc.

     : param_(param),
       use_skyrme_(conf.has_value({"Skyrme"})),
       use_symmetry_(conf.has_value({"Symmetry"})),
       use_coulomb_(conf.has_value({"Coulomb"})),
       use_vdf_(conf.has_value({"VDF"})) {
   if (use_skyrme_) {
     skyrme_a_ = conf.take({"Skyrme", "Skyrme_A"});
     skyrme_b_ = conf.take({"Skyrme", "Skyrme_B"});
     skyrme_tau_ = conf.take({"Skyrme", "Skyrme_Tau"});
   }
  
   if (use_symmetry_) {
     symmetry_S_Pot_ = conf.take({"Symmetry", "S_Pot"});
     if (conf.has_value({"Symmetry", "gamma"})) {
       symmetry_gamma_ = conf.take({"Symmetry", "gamma"});
       symmetry_is_rhoB_dependent_ = true;
     }
   }
   if (use_coulomb_) {
     coulomb_r_cut_ = conf.take({"Coulomb", "R_Cut"});
   }
   if (use_vdf_) {
     saturation_density_ = conf.take({"VDF", "Sat_rhoB"});
     std::vector<double> aux_coeffs = conf.take({"VDF", "Coeffs"});
     std::vector<double> aux_powers = conf.take({"VDF", "Powers"});
     if (aux_coeffs.size() != aux_powers.size()) {
       throw std::invalid_argument(
           "The number of coefficients should equal the number of powers.");
     }
     const int n_terms = aux_powers.size();
     for (int i = 0; i < n_terms; i++) {
       if (aux_powers[i] < 0.0) {
         throw std::invalid_argument("Powers need to be positive real numbers.");
       }
       // coefficients are provided in MeV, but the code uses GeV
       coeffs_.push_back(aux_coeffs[i] * mev_to_gev);
       powers_.push_back(aux_powers[i]);
     }
   }
 }

◆ ~Potentials()

smash::Potentials::~Potentials ( )

virtual

Standard destructor.

Definition at line 59 of file potentials.cc.

59 {}

Member Function Documentation

◆ skyrme_pot()

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

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 61 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 *
          (skyrme_a_ * std::abs(tmp) +
           skyrme_b_ * std::pow(std::abs(tmp), skyrme_tau_));
 }

◆ 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 79 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;
 }

◆ 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 71 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.;
   }
 }

◆ 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 90 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;
 }

◆ 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 112 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;
 }

◆ 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 145 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);
 }

◆ 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 153 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);
 }

◆ 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 171 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);
 }

◆ 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 217 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);
   }

◆ 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 235 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);
   }

◆ 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 186 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);
 }

◆ 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 226 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 261 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);
 }

◆ use_skyrme()

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

inlinevirtual

Returns: Is Skyrme potential on?

Definition at line 330 of file potentials.h.

330 { return use_skyrme_; }

◆ use_symmetry()

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

inlinevirtual

Returns: Is symmetry potential on?

Definition at line 332 of file potentials.h.

332 { return use_symmetry_; }

◆ use_coulomb()

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

inlinevirtual

Returns: Is Coulomb potential on?

Definition at line 334 of file potentials.h.

334 { return use_coulomb_; }

◆ skyrme_a()

double smash::Potentials::skyrme_a ( ) const

inline

Returns: Skyrme parameter skyrme_a, in MeV

Definition at line 337 of file potentials.h.

337 { return skyrme_a_; }

◆ skyrme_b()

double smash::Potentials::skyrme_b ( ) const

inline

Returns: Skyrme parameter skyrme_b, in MeV

Definition at line 339 of file potentials.h.

339 { return skyrme_b_; }

◆ skyrme_tau()

double smash::Potentials::skyrme_tau ( ) const

inline

Returns: Skyrme parameter skyrme_tau

Definition at line 341 of file potentials.h.

341 { return skyrme_tau_; }

◆ symmetry_S_pot()

double smash::Potentials::symmetry_S_pot ( ) const

inline

Returns: Skyrme parameter S_pot, in MeV

Definition at line 343 of file potentials.h.

343 { return symmetry_S_Pot_; }

◆ use_vdf()

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

inlinevirtual

Returns: Is VDF potential on?

Definition at line 346 of file potentials.h.

346 { return use_vdf_; }

◆ saturation_density()

double smash::Potentials::saturation_density ( ) const

inline

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

Definition at line 348 of file potentials.h.

348 { 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 350 of file potentials.h.

350 { return coeffs_; }

◆ powers()

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

inline

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

Definition at line 352 of file potentials.h.

352 { return powers_; }

◆ number_of_terms()

int smash::Potentials::number_of_terms ( ) const

inline

Returns: Number of terms in the VDF potential

Definition at line 354 of file potentials.h.

354 { return powers_.size(); }

◆ coulomb_r_cut()

double smash::Potentials::coulomb_r_cut ( ) const

inline

Returns: cutoff radius in ntegration for coulomb potential in fm

Definition at line 357 of file potentials.h.

357 { return coulomb_r_cut_; }

◆ 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 237 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;
   }
 }

◆ 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 247 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.;
   }
 }

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 365 of file potentials.h.

◆ use_skyrme_

bool smash::Potentials::use_skyrme_

private

Skyrme potential on/off.

Definition at line 368 of file potentials.h.

◆ use_symmetry_

bool smash::Potentials::use_symmetry_

private

Symmetry potential on/off.

Definition at line 371 of file potentials.h.

◆ use_coulomb_

bool smash::Potentials::use_coulomb_

private

Coulomb potential on/Off.

Definition at line 374 of file potentials.h.

◆ use_vdf_

bool smash::Potentials::use_vdf_

private

VDF potential on/off.

Definition at line 377 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 383 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 389 of file potentials.h.

◆ skyrme_tau_

double smash::Potentials::skyrme_tau_

private

Parameters of skyrme potentials: the power index.

Definition at line 395 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 398 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 404 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 411 of file potentials.h.

◆ coulomb_r_cut_

double smash::Potentials::coulomb_r_cut_

private

Cutoff in integration for coulomb potential.

Definition at line 414 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 420 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 422 of file potentials.h.

◆ powers_

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

private

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

Definition at line 424 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

Private Attributes

Constructor & Destructor Documentation

◆ Potentials()

◆ ~Potentials()

Member Function Documentation

◆ skyrme_pot()

◆ 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()

◆ skyrme_a()

◆ skyrme_b()

◆ skyrme_tau()

◆ symmetry_S_pot()

◆ use_vdf()

◆ saturation_density()

◆ coeffs()

◆ powers()

◆ number_of_terms()

◆ coulomb_r_cut()

◆ dVsym_drhoI3()

◆ dVsym_drhoB()

Member Data Documentation

◆ param_

◆ use_skyrme_

◆ use_symmetry_

◆ use_coulomb_

◆ use_vdf_

◆ skyrme_a_

◆ skyrme_b_

◆ skyrme_tau_

◆ symmetry_S_Pot_

◆ symmetry_is_rhoB_dependent_

◆ symmetry_gamma_

◆ coulomb_r_cut_

◆ saturation_density_

◆ coeffs_

◆ powers_