potentials.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2015,2017-2025
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_POTENTIALS_H_
#define SRC_INCLUDE_SMASH_POTENTIALS_H_
 
#include <tuple>
#include <utility>
#include <vector>
 
#include "gsl/gsl_multiroots.h"
#include "gsl/gsl_vector.h"
 
#include "configuration.h"
#include "density.h"
#include "forwarddeclarations.h"
#include "particledata.h"
#include "rootsolver.h"
#include "threevector.h"
 
namespace smash {
static constexpr int LPotentials = LogArea::Potentials::id;
class Potentials {
 public:
  Potentials(Configuration conf, const DensityParameters &parameters);
  virtual ~Potentials();
 
  ThreeVector single_particle_energy_gradient(DensityLattice *jB_lattice,
                                              const ThreeVector &position,
                                              const ThreeVector &momentum,
                                              double mass,
                                              ParticleList &plist) const {
    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;
  }
 
  double calculation_frame_energy(const ThreeVector &momentum,
                                  const FourVector &jmu_B, double mass) const {
    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.");
  }
 
  double skyrme_pot(const double baryon_density) const {
    return skyrme_pot(baryon_density, skyrme_a_, skyrme_b_, skyrme_tau_);
  }
 
  double symmetry_pot(const double baryon_isospin_density,
                      const double baryon_density) const;
 
  double symmetry_S(const double baryon_density) const;
 
  FourVector vdf_pot(double rhoB, const FourVector jmuB_net) const;
 
  double potential(const ThreeVector &r, const ParticleList &plist,
                   const ParticleType &acts_on) const;
 
  static std::pair<double, int> force_scale(const ParticleType &data);
 
  std::pair<ThreeVector, ThreeVector> skyrme_force(
      const double rhoB, const ThreeVector grad_j0B,
      const ThreeVector dvecjB_dt, const ThreeVector curl_vecjB) const;
 
  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;
 
  static ThreeVector E_field_integrand(ThreeVector pos,
                                       DensityOnLattice &charge_density,
                                       ThreeVector point) {
    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);
  }
 
  static ThreeVector B_field_integrand(ThreeVector pos,
                                       DensityOnLattice &charge_density,
                                       ThreeVector point) {
    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);
  }
  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;
 
  std::pair<ThreeVector, ThreeVector> vdf_force(
      const ThreeVector grad_A_0, const ThreeVector dA_dt,
      const ThreeVector curl_vecA) const;
 
  virtual std::tuple<ThreeVector, ThreeVector, ThreeVector, ThreeVector>
  all_forces(const ThreeVector &r, const ParticleList &plist) const;
 
  virtual bool use_skyrme() const { return use_skyrme_; }
  virtual bool use_symmetry() const { return use_symmetry_; }
  virtual bool use_coulomb() const { return use_coulomb_; }
  virtual bool use_momentum_dependence() const {
    return use_momentum_dependence_;
  }
 
  double skyrme_a() const { return skyrme_a_; }
  double skyrme_b() const { return skyrme_b_; }
  double skyrme_tau() const { return skyrme_tau_; }
  double symmetry_S_pot() const { return symmetry_S_Pot_; }
 
  virtual bool use_vdf() const { return use_vdf_; }
  double saturation_density() const { return saturation_density_; }
  const std::vector<double> &coeffs() const { return coeffs_; }
  const std::vector<double> &powers() const { return powers_; }
  int number_of_terms() const { return powers_.size(); }
 
  double coulomb_r_cut() const { return coulomb_r_cut_; }
 
  bool use_potentials_outside_lattice() const {
    return use_potentials_outside_lattice_;
  }
 
 private:
  const DensityParameters param_;
 
  bool use_skyrme_;
 
  bool use_symmetry_;
 
  bool use_coulomb_;
 
  bool use_vdf_;
 
  bool use_momentum_dependence_;
 
  double skyrme_a_;
 
  double skyrme_b_;
 
  double skyrme_tau_;
 
  double mom_dependence_Lambda_;
 
  double mom_dependence_C_;
 
  double symmetry_S_Pot_;
 
  bool symmetry_is_rhoB_dependent_ = false;
  double symmetry_gamma_;
 
  double coulomb_r_cut_;
 
  bool use_potentials_outside_lattice_;
 
  double saturation_density_;
  std::vector<double> coeffs_;
  std::vector<double> powers_;
 
  double dVsym_drhoI3(const double rhoB, const double rhoI3) const;
 
  double dVsym_drhoB(const double rhoB, const double rhoI3) const;
 
  static double skyrme_pot(const double baryon_density, const double A,
                           const double B, const double tau);
 
  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) {
    // 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;
  }
 
  static double momentum_dependent_part(double momentum, double rho, double C,
                                        double Lambda) {
    /* 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;
  }
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_POTENTIALS_H_