energy_inline.hpp

Go to the documentation of this file.

/*
 * Copyright (C) 2010-2026 The ESPResSo project
 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
 *   Max-Planck-Institute for Polymer Research, Theory Group
 *
 * This file is part of ESPResSo.
 *
 * ESPResSo is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * ESPResSo is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#pragma once
 
/** \file
 *  Energy calculation.
 */
 
#include "config/config.hpp"
 
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "electrostatics/coulomb_inline.hpp"
#include "magnetostatics/dipoles_inline.hpp"
#include "nonbonded_interactions/bmhtf-nacl.hpp"
#include "nonbonded_interactions/buckingham.hpp"
#include "nonbonded_interactions/gaussian.hpp"
#include "nonbonded_interactions/gay_berne.hpp"
#include "nonbonded_interactions/hat.hpp"
#include "nonbonded_interactions/hertzian.hpp"
#include "nonbonded_interactions/lj.hpp"
#include "nonbonded_interactions/ljcos.hpp"
#include "nonbonded_interactions/ljcos2.hpp"
#include "nonbonded_interactions/ljgen.hpp"
#include "nonbonded_interactions/morse.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "nonbonded_interactions/nonbonded_tab.hpp"
#include "nonbonded_interactions/smooth_step.hpp"
#include "nonbonded_interactions/soft_sphere.hpp"
#include "nonbonded_interactions/thole.hpp"
#include "nonbonded_interactions/wca.hpp"
 
#include "BoxGeometry.hpp"
#include "Observable_stat.hpp"
#include "Particle.hpp"
#include "bond_error.hpp"
#include "errorhandling.hpp"
#include "exclusions.hpp"
 
#include <utils/Vector.hpp>
 
#include <optional>
#include <span>
#include <string>
#include <variant>
 
inline double calc_central_radial_energy(IA_parameters const &ia_params,
                                         double const dist) {
 
  [[maybe_unused]] auto const mask = ia_params.active_pair_mask;
  double ret = 0.;
 
#ifdef ESPRESSO_LENNARD_JONES
  if (mask & pair_potential_bit(PairPotential::LennardJones))
    ret += lj_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_WCA
  if (mask & pair_potential_bit(PairPotential::WCA))
    ret += wca_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_LENNARD_JONES_GENERIC
  if (mask & pair_potential_bit(PairPotential::LennardJonesGeneric))
    ret += ljgen_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_SMOOTH_STEP
  if (mask & pair_potential_bit(PairPotential::SmoothStep))
    ret += SmSt_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_HERTZIAN
  if (mask & pair_potential_bit(PairPotential::Hertzian))
    ret += hertzian_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_GAUSSIAN
  if (mask & pair_potential_bit(PairPotential::Gaussian))
    ret += gaussian_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_BMHTF_NACL
  if (mask & pair_potential_bit(PairPotential::BMHTF))
    ret += BMHTF_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_MORSE
  if (mask & pair_potential_bit(PairPotential::Morse))
    ret += morse_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_BUCKINGHAM
  if (mask & pair_potential_bit(PairPotential::Buckingham))
    ret += buck_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_SOFT_SPHERE
  if (mask & pair_potential_bit(PairPotential::SoftSphere))
    ret += soft_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_HAT
  if (mask & pair_potential_bit(PairPotential::Hat))
    ret += hat_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_LJCOS2
  if (mask & pair_potential_bit(PairPotential::LJCos2))
    ret += ljcos2_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_TABULATED
  if (mask & pair_potential_bit(PairPotential::Tabulated))
    ret += tabulated_pair_energy(ia_params, dist);
#endif
 
#ifdef ESPRESSO_LJCOS
  if (mask & pair_potential_bit(PairPotential::LJCos))
    ret += ljcos_pair_energy(ia_params, dist);
#endif
 
  return ret;
}
 
/** Calculate non-bonded energies between a pair of particles.
 *  @param p1         particle 1.
 *  @param p2         particle 2.
 *  @param ia_params  the interaction parameters between the two particles
 *  @param d          vector between p1 and p2.
 *  @param dist       distance between p1 and p2.
 *  @param bonded_ias       bonded interaction kernels.
 *  @param coulomb          Electrostatics solver.
 *  @param coulomb_kernel   Coulomb energy kernel.
 *  @return the short-range interaction energy between the two particles
 */
inline double calc_non_bonded_pair_energy(
    Particle const &p1, Particle const &p2, IA_parameters const &ia_params,
    Utils::Vector3d const &d, double const dist,
    [[maybe_unused]] BondedInteractionsMap const &bonded_ias,
    [[maybe_unused]] Coulomb::Solver const &coulomb,
    [[maybe_unused]] Coulomb::ShortRangeEnergyKernel::kernel_type const
        *coulomb_kernel) {
 
  double ret = 0.;
 
  ret += calc_central_radial_energy(ia_params, dist);
 
#ifdef ESPRESSO_THOLE
  /* Thole damping */
  ret += thole_pair_energy(p1, p2, ia_params, d, dist, bonded_ias, coulomb,
                           coulomb_kernel);
#endif
 
#ifdef ESPRESSO_GAY_BERNE
  /* Gay-Berne */
  ret += gb_pair_energy(p1.quat(), p2.quat(), ia_params, d, dist);
#endif
 
  return ret;
}
 
inline std::optional<double> calc_pair_bonded_energy(
    Bonded_IA_Parameters const &iaparams, Utils::Vector3d const &dx,
    Utils::Vector3d const &pos1, Utils::Vector3d const &pos2, double q1q2,
    Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel) {
 
  if (auto const *iap = std::get_if<FeneBond>(&iaparams)) {
    return iap->energy(dx);
  }
  if (auto const *iap = std::get_if<HarmonicBond>(&iaparams)) {
    return iap->energy(dx);
  }
  if (auto const *iap = std::get_if<QuarticBond>(&iaparams)) {
    return iap->energy(dx);
  }
#ifdef ESPRESSO_ELECTROSTATICS
  if (auto const *iap = std::get_if<BondedCoulomb>(&iaparams)) {
    return iap->energy(q1q2, dx);
  }
  if (auto const *iap = std::get_if<BondedCoulombSR>(&iaparams)) {
    return iap->energy(pos1, pos2, dx, *kernel);
  }
#endif
#ifdef ESPRESSO_BOND_CONSTRAINT
  if (std::get_if<RigidBond>(&iaparams)) {
    return {0.};
  }
#endif
#ifdef ESPRESSO_TABULATED
  if (auto const *iap = std::get_if<TabulatedDistanceBond>(&iaparams)) {
    return iap->energy(dx);
  }
#endif
  if (std::get_if<VirtualBond>(&iaparams)) {
    return {0.};
  }
  throw BondUnknownTypeError();
}
 
inline std::optional<double>
calc_angle_bonded_energy(Bonded_IA_Parameters const &iaparams,
                         Utils::Vector3d const &vec1,
                         Utils::Vector3d const &vec2) {
  if (auto const *iap = std::get_if<AngleHarmonicBond>(&iaparams)) {
    return iap->energy(vec1, vec2);
  }
  if (auto const *iap = std::get_if<AngleCosineBond>(&iaparams)) {
    return iap->energy(vec1, vec2);
  }
  if (auto const *iap = std::get_if<AngleCossquareBond>(&iaparams)) {
    return iap->energy(vec1, vec2);
  }
  if (auto const *iap = std::get_if<TabulatedAngleBond>(&iaparams)) {
    return iap->energy(vec1, vec2);
  }
  if (std::get_if<IBMTriel>(&iaparams)) {
    runtimeWarningMsg() << "Unsupported bond type " +
                               std::to_string(iaparams.index()) +
                               " in energy calculation.";
    return 0.;
  }
  throw BondUnknownTypeError();
}
 
inline std::optional<double> calc_dihedral_bonded_energy(
    Bonded_IA_Parameters const &iaparams, Utils::Vector3d const &v12,
    Utils::Vector3d const &v23, Utils::Vector3d const &v34) {
  if (auto const *iap = std::get_if<DihedralBond>(&iaparams)) {
    return iap->energy(v12, v23, v34);
  }
  if (auto const *iap = std::get_if<TabulatedDihedralBond>(&iaparams)) {
    return iap->energy(v12, v23, v34);
  }
  if (std::get_if<IBMTribend>(&iaparams)) {
    runtimeWarningMsg() << "Unsupported bond type " +
                               std::to_string(iaparams.index()) +
                               " in energy calculation.";
    return 0.;
  }
  throw BondUnknownTypeError();
}
 
inline std::optional<double>
calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1,
                   std::span<Particle *> partners, BoxGeometry const &box_geo,
                   Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel) {
  auto const n_partners = static_cast<int>(partners.size());
 
  auto p2 = (n_partners > 0) ? partners[0] : nullptr;
  auto p3 = (n_partners > 1) ? partners[1] : nullptr;
  auto p4 = (n_partners > 2) ? partners[2] : nullptr;
 
  if (n_partners == 1) {
    auto const dx = box_geo.get_mi_vector(p1.pos(), p2->pos());
    return calc_pair_bonded_energy(iaparams, dx, p1.pos(), p2->pos(),
#ifdef ESPRESSO_ELECTROSTATICS
                                   p1.q() * p2->q(), kernel
#else
                                   0.0, nullptr
#endif
    );
  } // 1 partner
  if (n_partners == 2) {
    auto const vec1 = box_geo.get_mi_vector(p2->pos(), p1.pos());
    auto const vec2 = box_geo.get_mi_vector(p3->pos(), p1.pos());
    return calc_angle_bonded_energy(iaparams, vec1, vec2);
  } // 2 partners
  if (n_partners == 3) {
    // note: particles in a dihedral bond are ordered as p2-p1-p3-p4
    auto const v12 = box_geo.get_mi_vector(p1.pos(), p2->pos());
    auto const v23 = box_geo.get_mi_vector(p3->pos(), p1.pos());
    auto const v34 = box_geo.get_mi_vector(p4->pos(), p3->pos());
    return calc_dihedral_bonded_energy(iaparams, v12, v23, v34);
  } // 3 partners
  if (n_partners == 0) {
    return 0.;
  }
 
  throw BondInvalidSizeError(n_partners);
}
 
/** Calculate kinetic energies from translation for one particle.
 *  @param p   particle for which to calculate energies
 */
inline double translational_kinetic_energy(Particle const &p) {
  return p.is_virtual() ? 0. : 0.5 * p.mass() * p.v().norm2();
}
 
/** Calculate kinetic energies from rotation for one particle.
 *  @param p   particle for which to calculate energies
 */
inline double rotational_kinetic_energy([[maybe_unused]] Particle const &p) {
#ifdef ESPRESSO_ROTATION
  return (p.can_rotate() and not p.is_virtual())
             ? 0.5 * (hadamard_product(p.omega(), p.omega()) * p.rinertia())
             : 0.0;
#else
  return 0.0;
#endif
}