diff --git a/src/config/features.def b/src/config/features.def index 2045216378..4450b318d2 100644 --- a/src/config/features.def +++ b/src/config/features.def @@ -78,6 +78,7 @@ BUCKINGHAM SOFT_SPHERE HAT GAY_BERNE implies ROTATION +GAY_BERNE_WIDTH implies ROTATION THOLE requires ELECTROSTATICS /* ScaFaCoS */ diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 2995150aef..43276f1f1c 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -34,6 +34,7 @@ #include "nonbonded_interactions/buckingham.hpp" #include "nonbonded_interactions/gaussian.hpp" #include "nonbonded_interactions/gay_berne.hpp" +#include "nonbonded_interactions/gay_berne_width.hpp" #include "nonbonded_interactions/hat.hpp" #include "nonbonded_interactions/hertzian.hpp" #include "nonbonded_interactions/lj.hpp" @@ -175,6 +176,10 @@ inline double calc_non_bonded_pair_energy( ret += gb_pair_energy(p1.quat(), p2.quat(), ia_params, d, dist); #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + ret += gb_width_pair_energy(p1.quat(), p2.quat(), ia_params, d, dist); +#endif + return ret; } diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 19cf232faf..d73e6f9b96 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -40,6 +40,7 @@ #include "nonbonded_interactions/buckingham.hpp" #include "nonbonded_interactions/gaussian.hpp" #include "nonbonded_interactions/gay_berne.hpp" +#include "nonbonded_interactions/gay_berne_width.hpp" #include "nonbonded_interactions/hat.hpp" #include "nonbonded_interactions/hertzian.hpp" #include "nonbonded_interactions/lj.hpp" @@ -152,6 +153,11 @@ inline ParticleForce calc_non_central_force(Particle const &p1, #ifdef ESPRESSO_GAY_BERNE pf += gb_pair_force(p1.quat(), p2.quat(), ia_params, d, dist); #endif + +#ifdef ESPRESSO_GAY_BERNE_WIDTH + pf += gb_width_pair_force(p1.quat(), p2.quat(), ia_params, d, dist); +#endif + return pf; } @@ -166,6 +172,11 @@ inline Utils::Vector3d calc_non_central_force(Utils::Vector3d const &dir1, #ifdef ESPRESSO_GAY_BERNE f += gb_pair_force(dir1, dir2, ia_params, d, dist).f; #endif + +#ifdef ESPRESSO_GAY_BERNE_WIDTH + f += gb_width_pair_force(dir1, dir2, ia_params, d, dist).f; +#endif + return f; } diff --git a/src/core/nonbonded_interactions/CMakeLists.txt b/src/core/nonbonded_interactions/CMakeLists.txt index d008c69292..6ce2d24727 100644 --- a/src/core/nonbonded_interactions/CMakeLists.txt +++ b/src/core/nonbonded_interactions/CMakeLists.txt @@ -23,6 +23,7 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/buckingham.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gaussian.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gay_berne.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/gay_berne_width.cpp ${CMAKE_CURRENT_SOURCE_DIR}/hat.cpp ${CMAKE_CURRENT_SOURCE_DIR}/hertzian.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ljcos2.cpp diff --git a/src/core/nonbonded_interactions/gay_berne_width.cpp b/src/core/nonbonded_interactions/gay_berne_width.cpp new file mode 100644 index 0000000000..a0892ce93d --- /dev/null +++ b/src/core/nonbonded_interactions/gay_berne_width.cpp @@ -0,0 +1,41 @@ +/* + * Copyright (C) 2010-2026 The ESPResSo project + * + * 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 . + */ + +/** \file + * + * Implementation of \ref gay_berne_width.hpp + */ + +#include "gay_berne_width.hpp" + +#ifdef ESPRESSO_GAY_BERNE_WIDTH + +#include "nonbonded_interaction_data.hpp" + +#include + +GayBerneWidth_Parameters::GayBerneWidth_Parameters(double eps, double sig, + double wid, double cut, + double k1, double k2, + double mu, double nu) + : eps{eps}, sig{sig}, wid{wid}, cut{cut}, k1{k1}, k2{k2}, mu{mu}, nu{nu}, + chi1{((k1 * k1) - 1.) / ((k1 * k1) + 1.)}, + chi2{(std::pow(k2, 1. / mu) - 1.) / (std::pow(k2, 1. / mu) + 1.)} {} + +#endif // ESPRESSO_GAY_BERNE_WIDTH diff --git a/src/core/nonbonded_interactions/gay_berne_width.hpp b/src/core/nonbonded_interactions/gay_berne_width.hpp new file mode 100644 index 0000000000..fd293b1c24 --- /dev/null +++ b/src/core/nonbonded_interactions/gay_berne_width.hpp @@ -0,0 +1,384 @@ +/* + * Copyright (C) 2010-2026 The ESPResSo project + * + * 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 . + */ + +#pragma once + +/** \file + * Routines to calculate a modified Gay--Berne potential with an + * independently adjustable radial width. + * + * The original ESPResSo Gay--Berne interaction uses the radial + * coordinate + * + * \f[ + * x = \frac{r-s+\sigma}{\sigma}, + * \f] + * + * where \f$s\f$ is the orientation-dependent contact distance. + * + * This interaction replaces it with + * + * \f[ + * x_w = + * 2^{1/6} + \frac{r-r_{\min}}{w}, + * \f] + * + * where + * + * \f[ + * r_{\min} + * = + * s + \left(2^{1/6}-1\right)\sigma. + * \f] + * + * Therefore, \f$\sigma\f$ retains control of the position of the + * potential minimum, while \f$w\f$ controls the radial width and + * local stiffness of the well. + * + * Setting \f$w=\sigma\f$ recovers the radial form of the original + * ESPResSo Gay--Berne interaction. + */ + +#include "config/config.hpp" + +#ifdef ESPRESSO_GAY_BERNE_WIDTH + +#include "Particle.hpp" +#include "nonbonded_interaction_data.hpp" + +#include +#include +#include +#include + +#include + +/** + * Calculate the force and torque for the modified Gay--Berne interaction. + * + * This is the main implementation and operates directly on the particle + * directors. The director-based form is also required by pressure and + * virial calculations. + * + * @param ui Director of particle i. + * @param uj Director of particle j. + * @param ia_params Non-bonded interaction parameters. + * @param d Distance vector between the particles. + * @param dist Magnitude of the distance vector. + * + * @return Force acting on the particle and its torque. + */ +inline ParticleForce gb_width_pair_force(Utils::Vector3d const &ui, + Utils::Vector3d const &uj, + IA_parameters const &ia_params, + Utils::Vector3d const &d, + double dist) { + using Utils::int_pow; + using Utils::sqr; + + auto const ¶ms = ia_params.gay_berne_width; + + if (dist >= params.cut) { + return {}; + } + + auto const e0 = params.eps; + + /* + * sig controls the reference minimum position. + */ + auto const s0 = params.sig; + + /* + * wid controls the radial width and stiffness of the well. + */ + auto const w0 = params.wid; + + auto const chi1 = params.chi1; + auto const chi2 = params.chi2; + auto const mu = params.mu; + auto const nu = params.nu; + + /* + * Unit vector along the particle-particle separation. + */ + auto const r = Utils::Vector3d(d).normalize(); + + auto const dui = d * ui; + auto const duj = d * uj; + + auto const rui = r * ui; + auto const ruj = r * uj; + auto const uij = ui * uj; + + auto const oo1 = (dui + duj) / (1. + chi1 * uij); + auto const oo2 = (dui - duj) / (1. - chi1 * uij); + + auto const tt1 = (dui + duj) / (1. + chi2 * uij); + auto const tt2 = (dui - duj) / (1. - chi2 * uij); + + auto const o1 = sqr(rui + ruj) / (1. + chi1 * uij); + auto const o2 = sqr(rui - ruj) / (1. - chi1 * uij); + + auto const t1 = sqr(rui + ruj) / (1. + chi2 * uij); + auto const t2 = sqr(rui - ruj) / (1. - chi2 * uij); + + auto const Brhi1 = chi1 * (o1 + o2); + auto const Brhi2 = chi2 * (t1 + t2); + + auto const e1 = 1. / (1. - sqr(chi1 * uij)); + auto const e2 = 1. - 0.5 * Brhi2; + + /* + * As in ESPResSo's original Gay--Berne force implementation, + * this energy prefactor already contains the Lennard-Jones factor 4. + */ + auto const e = 4. * e0 * std::pow(e1, 0.5 * nu) * std::pow(e2, mu); + + /* + * Orientation-dependent contact distance: + * + * s = sig / sqrt(1 - Brhi1 / 2). + */ + auto const s1 = 1. / std::sqrt(1. - 0.5 * Brhi1); + auto const s = s0 * s1; + + /* + * Dimensionless position of the Lennard-Jones minimum. + */ + auto const alpha = std::pow(2., 1. / 6.); + + /* + * Preserve the original Gay--Berne minimum position. + */ + auto const r_min = s + (alpha - 1.) * s0; + + /* + * Modified inverse radial coordinate: + * + * X = 1 / x_w + * + * where + * + * x_w = alpha + (r - r_min) / wid. + * + * When wid == sig: + * + * x_w = (r - s + sig) / sig, + * + * which is the original ESPResSo Gay--Berne radial coordinate. + */ + auto const X = 1. / (alpha + (dist - r_min) / w0); + + /* + * Value of the modified inverse radial coordinate at the cutoff. + * This is used to shift the energy continuously to zero at the cutoff. + */ + auto const Xcut = 1. / (alpha + (params.cut - r_min) / w0); + + auto const X6 = int_pow<6>(X); + auto const Xcut6 = int_pow<6>(Xcut); + + /* + * Lennard-Jones-like radial terms: + * + * Brack = X^12 - X^6 + * + * and + * + * Bra12 = -d(Brack)/dx_w + * = 6 X^7 (2 X^6 - 1). + */ + auto const Bra12 = 6. * X6 * X * (2. * X6 - 1.); + auto const Bra12Cut = 6. * Xcut6 * Xcut * (2. * Xcut6 - 1.); + + auto const Brack = X6 * (X6 - 1.); + auto const BrackCut = Xcut6 * (Xcut6 - 1.); + + /* + * Angular energy-prefactor derivatives are unchanged from the original + * Gay--Berne interaction. + */ + auto Koef1 = mu / e2; + auto Koef2 = int_pow<3>(s1) * 0.5; + + /* + * The contact-distance derivative terms acquire sig / wid because + * + * dx_w / ds = -1 / wid + * + * instead of + * + * dx / ds = -1 / sig. + */ + auto const width_scale = s0 / w0; + + /* Radial derivative. */ + auto const dU_dr = + e * + (Koef1 * Brhi2 * (Brack - BrackCut) - + width_scale * Koef2 * Brhi1 * (Bra12 - Bra12Cut) - Bra12 * dist / w0) / + sqr(dist); + + Koef1 *= chi2 / sqr(dist); + Koef2 *= chi1 / sqr(dist); + + /* + * Derivatives with respect to the orientational scalar products. + */ + auto const dU_da = + e * (Koef1 * (tt1 + tt2) * (BrackCut - Brack) + + width_scale * Koef2 * (oo1 + oo2) * (Bra12 - Bra12Cut)); + + auto const dU_db = + e * (Koef1 * (tt2 - tt1) * (Brack - BrackCut) + + width_scale * Koef2 * (oo1 - oo2) * (Bra12 - Bra12Cut)); + + auto const dU_dc = + e * ((Brack - BrackCut) * (nu * e1 * sqr(chi1) * uij + + 0.5 * Koef1 * chi2 * (sqr(tt1) - sqr(tt2))) - + width_scale * (Bra12 - Bra12Cut) * 0.5 * Koef2 * chi1 * + (sqr(oo1) - sqr(oo2))); + + /* + * Force and torque structure follows the original ESPResSo + * Gay--Berne implementation. + * + * The torque is + * + * tau_i = ui x G2. + */ + auto const G2 = -dU_da * d - dU_dc * uj; + + ParticleForce pf{ + -dU_dr * d - dU_da * ui - dU_db * uj, + vector_product(ui, G2), + }; + + return pf; +} + +/** + * Calculate the energy of the modified Gay--Berne interaction. + * + * This is the main energy implementation and operates directly on + * particle directors. + * + * @param ui Director of particle i. + * @param uj Director of particle j. + * @param ia_params Non-bonded interaction parameters. + * @param d Distance vector between the particles. + * @param dist Magnitude of the distance vector. + * + * @return Interaction energy. + */ +inline double gb_width_pair_energy(Utils::Vector3d const &ui, + Utils::Vector3d const &uj, + IA_parameters const &ia_params, + Utils::Vector3d const &d, double dist) { + using Utils::int_pow; + using Utils::sqr; + + auto const ¶ms = ia_params.gay_berne_width; + + if (dist >= params.cut) { + return {}; + } + + auto const e0 = params.eps; + auto const s0 = params.sig; + auto const w0 = params.wid; + + auto const chi1 = params.chi1; + auto const chi2 = params.chi2; + auto const mu = params.mu; + auto const nu = params.nu; + + auto const r = Utils::Vector3d(d).normalize(); + + auto const uij = ui * uj; + auto const rui = r * ui; + auto const ruj = r * uj; + + auto const o1 = sqr(rui + ruj) / (1. + chi1 * uij); + auto const o2 = sqr(rui - ruj) / (1. - chi1 * uij); + + auto const t1 = sqr(rui + ruj) / (1. + chi2 * uij); + auto const t2 = sqr(rui - ruj) / (1. - chi2 * uij); + + /* + * Orientation-dependent energy prefactor. + */ + auto const e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu); + + auto const e2 = std::pow(1. - 0.5 * chi2 * (t1 + t2), mu); + + auto const e = e0 * e1 * e2; + + /* + * Orientation-dependent contact distance. + */ + auto const s1 = 1. / std::sqrt(1. - 0.5 * chi1 * (o1 + o2)); + + auto const s = s0 * s1; + + auto const alpha = std::pow(2., 1. / 6.); + + /* + * Original Gay--Berne minimum position. + */ + auto const r_min = s + (alpha - 1.) * s0; + + auto r_eff = [=](double r) { return alpha + (r - r_min) / w0; }; + + auto E = [=](double x) { + return 4. * e * (int_pow<12>(1. / x) - int_pow<6>(1. / x)); + }; + + return E(r_eff(dist)) - E(r_eff(params.cut)); +} + +/** + * Quaternion wrapper for the modified Gay--Berne force and torque. + */ +inline ParticleForce gb_width_pair_force(Utils::Quaternion const &qi, + Utils::Quaternion const &qj, + IA_parameters const &ia_params, + Utils::Vector3d const &d, + double dist) { + auto const ui = Utils::convert_quaternion_to_director(qi); + auto const uj = Utils::convert_quaternion_to_director(qj); + + return gb_width_pair_force(ui, uj, ia_params, d, dist); +} + +/** + * Quaternion wrapper for the modified Gay--Berne energy. + */ +inline double gb_width_pair_energy(Utils::Quaternion const &qi, + Utils::Quaternion const &qj, + IA_parameters const &ia_params, + Utils::Vector3d const &d, double dist) { + auto const ui = Utils::convert_quaternion_to_director(qi); + auto const uj = Utils::convert_quaternion_to_director(qj); + + return gb_width_pair_energy(ui, uj, ia_params, d, dist); +} + +#endif // ESPRESSO_GAY_BERNE_WIDTH diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index d98cf0263e..900f97d3a6 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -128,6 +128,11 @@ recalc_maximal_cutoff(IA_parameters const &data, data.gay_berne.max_cutoff()); #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + consider(max_cut_current, mask, PairPotential::GayBerneWidth, + data.gay_berne_width.max_cutoff()); +#endif + #ifdef ESPRESSO_TABULATED consider(max_cut_current, mask, PairPotential::Tabulated, data.tab.cutoff()); #endif diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index d09379aee4..7503d89c65 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -242,6 +242,27 @@ struct GayBerne_Parameters { double max_cutoff() const { return cut; } }; +/** Gay-Berne potential with independently adjustable radial width */ +struct GayBerneWidth_Parameters { + double eps = 0.0; + double sig = 0.0; + double wid = 0.0; + double cut = inactive_cutoff; + double k1 = 0.0; + double k2 = 0.0; + double mu = 0.0; + double nu = 0.0; + double chi1 = 0.0; + double chi2 = 0.0; + + GayBerneWidth_Parameters() = default; + + GayBerneWidth_Parameters(double eps, double sig, double wid, double cut, + double k1, double k2, double mu, double nu); + + double max_cutoff() const { return cut; } +}; + /** Thole potential */ struct Thole_Parameters { double scaling_coeff = 0.; // inactive cutoff is 0 @@ -298,6 +319,7 @@ enum class PairPotential : unsigned { LJCos2, Tabulated, GayBerne, + GayBerneWidth, DPD, }; @@ -375,6 +397,10 @@ struct IA_parameters { GayBerne_Parameters gay_berne; #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + GayBerneWidth_Parameters gay_berne_width; +#endif + #ifdef ESPRESSO_TABULATED TabulatedPotential tab; #endif diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index fd341a1f85..03d149d6a5 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -320,6 +320,46 @@ def default_params(self): return {} +@script_interface_register +class GayBerneWidthInteraction(NonBondedInteraction): + """Gay--Berne interaction with independently adjustable radial width. + + Methods + ------- + set_params() + Set new parameters for the interaction. + + Parameters + ---------- + eps : :obj:`float` + Potential well-depth scale. + sig : :obj:`float` + Gay--Berne reference length scale. Together with the particle + orientations and ``k1``, it determines the position of the + potential minimum. + wid : :obj:`float` + Radial width of the Lennard-Jones-like well. Setting ``wid`` + equal to ``sig`` recovers the radial form of the original + ESPResSo Gay--Berne interaction. + cut : :obj:`float` + Cutoff distance. + k1 : :obj:`float` + Molecular elongation parameter. + k2 : :obj:`float` + Ratio controlling the orientation dependence of the well depth. + mu : :obj:`float` + Gay--Berne energy-anisotropy exponent. + nu : :obj:`float` + Gay--Berne orientation exponent. + """ + + _so_name = "Interactions::InteractionGayBerneWidth" + _so_feature = "GAY_BERNE_WIDTH" + + def default_params(self): + return {} + + @script_interface_register class TabulatedNonBonded(NonBondedInteraction): """Tabulated interaction. diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index 70585b9214..87a6616ca2 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/src/script_interface/interactions/NonBondedInteraction.hpp @@ -567,6 +567,39 @@ class InteractionGayBerne }; #endif // ESPRESSO_GAY_BERNE +#ifdef ESPRESSO_GAY_BERNE_WIDTH +class InteractionGayBerneWidth + : public InteractionPotentialInterface<::GayBerneWidth_Parameters> { +protected: + CoreInteraction IA_parameters::*get_ptr_offset() const override { + return &::IA_parameters::gay_berne_width; + } + +public: + InteractionGayBerneWidth() { + add_parameters({ + make_autoparameter(&CoreInteraction::eps, "eps"), + make_autoparameter(&CoreInteraction::sig, "sig"), + make_autoparameter(&CoreInteraction::wid, "wid"), + make_autoparameter(&CoreInteraction::cut, "cut"), + make_autoparameter(&CoreInteraction::k1, "k1"), + make_autoparameter(&CoreInteraction::k2, "k2"), + make_autoparameter(&CoreInteraction::mu, "mu"), + make_autoparameter(&CoreInteraction::nu, "nu"), + }); + } + +private: + std::string inactive_parameter() const override { return "cut"; } + + void make_new_instance(VariantMap const ¶ms) override { + m_handle = make_shared_from_args( + params, "eps", "sig", "wid", "cut", "k1", "k2", "mu", "nu"); + } +}; +#endif // ESPRESSO_GAY_BERNE_WIDTH + #ifdef ESPRESSO_TABULATED class InteractionTabulated : public InteractionPotentialInterface<::TabulatedPotential> { @@ -754,6 +787,9 @@ class NonBondedInteractionHandle #ifdef ESPRESSO_GAY_BERNE std::shared_ptr m_gay_berne; #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + std::shared_ptr m_gay_berne_width; +#endif #ifdef ESPRESSO_TABULATED std::shared_ptr m_tabulated; #endif @@ -873,6 +909,10 @@ class NonBondedInteractionHandle #ifdef ESPRESSO_GAY_BERNE fun(m_gay_berne, "gay_berne", "Interactions::InteractionGayBerne"); #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + fun(m_gay_berne_width, "gay_berne_width", + "Interactions::InteractionGayBerneWidth"); +#endif #ifdef ESPRESSO_TABULATED fun(m_tabulated, "tabulated", "Interactions::InteractionTabulated"); #endif diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index 5d6e4a22de..d2c110c7a6 100644 --- a/src/script_interface/interactions/initialize.cpp +++ b/src/script_interface/interactions/initialize.cpp @@ -97,6 +97,10 @@ void initialize(Utils::Factory *om) { #ifdef ESPRESSO_GAY_BERNE om->register_new("Interactions::InteractionGayBerne"); #endif +#ifdef ESPRESSO_GAY_BERNE_WIDTH + om->register_new( + "Interactions::InteractionGayBerneWidth"); +#endif #ifdef ESPRESSO_TABULATED om->register_new("Interactions::InteractionTabulated"); #endif