diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index 58e4d965ee..59b960f391 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/doc/sphinx/inter_non-bonded.rst @@ -545,6 +545,49 @@ caution when performing energy calculations. However, you can often choose the cutoff such that the energy difference at the cutoff is less than a desired accuracy, since the potential decays very rapidly. +.. _Anisotropic Gaussian interaction: + +Anisotropic Gaussian interaction +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. note:: + Feature ``GAUSSIAN_ANISO`` required. + +The interface for the anisotropic Gaussian interaction is implemented in +:class:`espressomd.interactions.GaussianAnisoInteraction`. The anisotropic +Gaussian interaction parameters can be set via:: + + system.non_bonded_inter[type1, type2].gaussian_aniso.set_params( + eps=1.0, sig_x=0.1, sig_y=0.2, sig_z=0.3, cutoff=1.0) + +This defines an anisotropic Gaussian interaction between particles of the +types ``type1`` and ``type2``. The potential is defined by + +.. math:: + + V(\Delta x, \Delta y, \Delta z) = + \begin{cases} + \epsilon \exp\left[ + -\frac{1}{2}\left( + \frac{\Delta x^2}{\sigma_x^2} + + \frac{\Delta y^2}{\sigma_y^2} + + \frac{\Delta z^2}{\sigma_z^2} + \right)\right] + & r < r_\mathrm{cut}\\ + 0 & r \ge r_\mathrm{cut} + \end{cases} + +where :math:`r = \sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}`. The +parameters :math:`\sigma_x`, :math:`\sigma_y` and :math:`\sigma_z` set the +width of the Gaussian interaction along the Cartesian coordinate axes, while +``cutoff`` defines a radial cutoff based on the particle separation distance. + +For :math:`\sigma_x = \sigma_y = \sigma_z`, this reduces to an isotropic +three-dimensional Gaussian interaction with a radial cutoff. Currently, there +is no shift implemented, which means that the potential is discontinuous at +:math:`r=r_\mathrm{cut}`. Therefore use caution when performing energy +calculations. + .. _DPD interaction: DPD interaction diff --git a/myconfig.hpp b/myconfig.hpp new file mode 100644 index 0000000000..c14e424f30 --- /dev/null +++ b/myconfig.hpp @@ -0,0 +1,63 @@ + + +/* Generic features */ +#define EXTERNAL_FORCES +#define MASS +#define EXCLUSIONS +#define BOND_CONSTRAINT +#define THERMOSTAT_PER_PARTICLE +#define COLLISION_DETECTION +#define NPT +#define ENGINE +#define PARTICLE_ANISOTROPY +#define H5MD + +/* Rotation */ +#define ROTATION +#define ROTATIONAL_INERTIA + +/* Electrostatics */ +#define ELECTROSTATICS + +/* Magnetostatics */ +#define DIPOLES +#define DIPOLE_FIELD_TRACKING +#define THERMAL_STONER_WOHLFARTH + +/* Virtual sites features */ +#define VIRTUAL_SITES +#define VIRTUAL_SITES_RELATIVE +#define VIRTUAL_SITES_CENTER_OF_MASS +#define VIRTUAL_SITES_INERTIALESS_TRACERS + +/* DPD features */ +#define DPD + +/* Lattice-Boltzmann features */ +#define LB_ELECTROHYDRODYNAMICS + +/* Interaction features */ +#define TABULATED +#define LENNARD_JONES +#define WCA +#define LENNARD_JONES_GENERIC +#define LJCOS +#define LJCOS2 +#define LJGEN_SOFTCORE +#define SMOOTH_STEP +#define HERTZIAN +#define GAUSSIAN +#define GAUSSIAN_ANISO +#define BMHTF_NACL +#define MORSE +#define BUCKINGHAM +#define SOFT_SPHERE +#define HAT +#define GAY_BERNE +#define THOLE + +/* ScaFaCoS */ +#define SCAFACOS_DIPOLES + +/* Debugging */ +#define ADDITIONAL_CHECKS diff --git a/src/config/features.def b/src/config/features.def index 2045216378..de7f22bda4 100644 --- a/src/config/features.def +++ b/src/config/features.def @@ -72,6 +72,7 @@ LJGEN_SOFTCORE implies LENNARD_JONES_GENERIC SMOOTH_STEP HERTZIAN GAUSSIAN +GAUSSIAN_ANISO BMHTF_NACL MORSE BUCKINGHAM diff --git a/src/core/energy_cabana.hpp b/src/core/energy_cabana.hpp index a67d4c9ee0..494fc42545 100644 --- a/src/core/energy_cabana.hpp +++ b/src/core/energy_cabana.hpp @@ -154,6 +154,13 @@ struct EnergyKernel { { e_nb += calc_central_radial_energy(ia_params, dist); + // Only call GAUSSIAN_ANISO energy kernel if active +#ifdef ESPRESSO_GAUSSIAN_ANISO + if (ia_params.active_pair_mask & + pair_potential_bit(PairPotential::GaussianAniso)) { + e_nb += gaussian_aniso_pair_energy(ia_params, d); + } +#endif // Only call Thole energy kernel if active #ifdef ESPRESSO_THOLE if (thole_active(ia_params, coulomb_u_kernel != nullptr)) { diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 2995150aef..ea9a0196c0 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -33,6 +33,7 @@ #include "nonbonded_interactions/bmhtf-nacl.hpp" #include "nonbonded_interactions/buckingham.hpp" #include "nonbonded_interactions/gaussian.hpp" +#include "nonbonded_interactions/gaussian_aniso.hpp" #include "nonbonded_interactions/gay_berne.hpp" #include "nonbonded_interactions/hat.hpp" #include "nonbonded_interactions/hertzian.hpp" @@ -164,6 +165,13 @@ inline double calc_non_bonded_pair_energy( ret += calc_central_radial_energy(ia_params, dist); +#ifdef ESPRESSO_GAUSSIAN_ANISO + if (ia_params.active_pair_mask & + pair_potential_bit(PairPotential::GaussianAniso)) { + ret += gaussian_aniso_pair_energy(ia_params, d); + } +#endif + #ifdef ESPRESSO_THOLE /* Thole damping */ ret += thole_pair_energy(p1, p2, ia_params, d, dist, bonded_ias, coulomb, diff --git a/src/core/forces_cabana.hpp b/src/core/forces_cabana.hpp index c75cbba78b..f19e9880ac 100644 --- a/src/core/forces_cabana.hpp +++ b/src/core/forces_cabana.hpp @@ -163,6 +163,14 @@ struct ForcesKernel { if (not skip_non_bonded) { pf.f += calc_central_radial_force(ia_params, d, dist); + // Only call ESPRESSO_ANISO force kernel if active +#ifdef ESPRESSO_GAUSSIAN_ANISO + if (ia_params.active_pair_mask & + pair_potential_bit(PairPotential::GaussianAniso)) { + pf.f += gaussian_aniso_pair_force(ia_params, d); + } +#endif + // Only call Thole force kernel if active #ifdef ESPRESSO_THOLE if (thole_active(ia_params, coulomb_kernel != nullptr)) { diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 19cf232faf..74ba4fccc4 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -39,6 +39,7 @@ #include "nonbonded_interactions/bmhtf-nacl.hpp" #include "nonbonded_interactions/buckingham.hpp" #include "nonbonded_interactions/gaussian.hpp" +#include "nonbonded_interactions/gaussian_aniso.hpp" #include "nonbonded_interactions/gay_berne.hpp" #include "nonbonded_interactions/hat.hpp" #include "nonbonded_interactions/hertzian.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_GAUSSIAN_ANISO + pf.f += gaussian_aniso_pair_force(ia_params, d); +#endif + return pf; } @@ -163,9 +169,15 @@ inline Utils::Vector3d calc_non_central_force(Utils::Vector3d const &dir1, Utils::Vector3d const &d, double const dist) { Utils::Vector3d f{}; + #ifdef ESPRESSO_GAY_BERNE f += gb_pair_force(dir1, dir2, ia_params, d, dist).f; #endif + +#ifdef ESPRESSO_GAUSSIAN_ANISO + f += gaussian_aniso_pair_force(ia_params, d); +#endif + return f; } diff --git a/src/core/nonbonded_interactions/CMakeLists.txt b/src/core/nonbonded_interactions/CMakeLists.txt index d008c69292..c4ce443ef7 100644 --- a/src/core/nonbonded_interactions/CMakeLists.txt +++ b/src/core/nonbonded_interactions/CMakeLists.txt @@ -22,6 +22,7 @@ target_sources( PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/bmhtf-nacl.cpp ${CMAKE_CURRENT_SOURCE_DIR}/buckingham.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gaussian.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/gaussian_aniso.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gay_berne.cpp ${CMAKE_CURRENT_SOURCE_DIR}/hat.cpp ${CMAKE_CURRENT_SOURCE_DIR}/hertzian.cpp diff --git a/src/core/nonbonded_interactions/gaussian_aniso.cpp b/src/core/nonbonded_interactions/gaussian_aniso.cpp new file mode 100644 index 0000000000..0385c16b98 --- /dev/null +++ b/src/core/nonbonded_interactions/gaussian_aniso.cpp @@ -0,0 +1,49 @@ +/* + * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 gaussian_aniso.hpp + */ + +#include "gaussian_aniso.hpp" + +#ifdef ESPRESSO_GAUSSIAN_ANISO +#include "nonbonded_interaction_data.hpp" + +#include + +GaussianAniso_Parameters::GaussianAniso_Parameters(double eps, double sig_x, + double sig_y, double sig_z, + double cutoff) + : eps{eps}, sig_x{sig_x}, sig_y{sig_y}, sig_z{sig_z}, cut{cutoff} { + if (sig_x <= 0.) { + throw std::domain_error("GaussianAniso parameter 'sig_x' has to be >= 0"); + } + if (sig_y <= 0.) { + throw std::domain_error("GaussianAniso parameter 'sig_y' has to be >= 0"); + } + if (sig_z <= 0.) { + throw std::domain_error("GaussianAniso parameter 'sig_z' has to be >= 0"); + } + if (cutoff <= 0.) { + throw std::domain_error("GaussianAniso parameter 'cutoff' has to be >= 0"); + } +} + +#endif // ESPRESSO_GAUSSIAN_ANISO diff --git a/src/core/nonbonded_interactions/gaussian_aniso.hpp b/src/core/nonbonded_interactions/gaussian_aniso.hpp new file mode 100644 index 0000000000..11fff1c8d5 --- /dev/null +++ b/src/core/nonbonded_interactions/gaussian_aniso.hpp @@ -0,0 +1,128 @@ +/* + * Copyright (C) 2010-2022 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 . + */ +#ifndef GAUSSIAN_ANISO_H +#define GAUSSIAN_ANISO_H + +/** \file + * Routines to calculate an anisotropic 3D Gaussian potential between + * particle pairs with independent widths sigma_x, sigma_y, sigma_z. + * + * Potential: + * U(dx, dy, dz) = eps * exp( -0.5 * [ (dx/sig_x)^2 + * + (dy/sig_y)^2 + * + (dz/sig_z)^2 ] ) + * + * If sig_x == sig_y == sig_z, this reduces to the ordinary isotropic + * 3D Gaussian (up to using a radial cutoff on |r|, same semantics + * as the built-in GAUSSIAN potential). + * + * Implementation in \ref gaussian_aniso.cpp. + */ + +#include "config/config.hpp" + +#include "nonbonded_interaction_data.hpp" + +#include +#include +#include + +#ifdef ESPRESSO_GAUSSIAN_ANISO + +/** Calculate anisotropic Gaussian energy. + * + * Parameters from ia_params.gaussian_aniso: + * eps : amplitude + * sig_x : width along x + * sig_y : width along y + * sig_z : width along z + * cut : radial cutoff on |r| (same semantics as GAUSSIAN) + * + * Input: + * dx, dy, dz : displacement components (r_j - r_i) + * + * Returns: + * U_ij + */ +inline double gaussian_aniso_pair_energy(IA_parameters const &ia_params, + Utils::Vector3d const &d) { + double dx = d[0]; + double dy = d[1]; + double dz = d[2]; + + auto const &p = ia_params.gaussian_aniso; + + const double r2 = Utils::sqr(dx) + Utils::sqr(dy) + Utils::sqr(dz); + + if (r2 >= Utils::sqr(p.cut)) { + return 0.0; + } + + const double sig_x2 = Utils::sqr(p.sig_x); + const double sig_y2 = Utils::sqr(p.sig_y); + const double sig_z2 = Utils::sqr(p.sig_z); + + const double A = 0.5 * (Utils::sqr(dx) / sig_x2 + Utils::sqr(dy) / sig_y2 + + Utils::sqr(dz) / sig_z2); + + return p.eps * std::exp(-A); +} + +/** Calculate anisotropic Gaussian force. + * + * This is the analogue of gaussian_pair_force_factor, but because + * the potential is anisotropic, the force is not simply parallel to r, + * so we return the components directly. + * + * Input: + * dx, dy, dz : displacement components (r_j - r_i) + * + * Output: + * Force vector on particle i due to particle j + */ +inline Utils::Vector3d gaussian_aniso_pair_force(IA_parameters const &ia_params, + Utils::Vector3d const &d) { + auto const &p = ia_params.gaussian_aniso; + + const double dx = d[0]; + const double dy = d[1]; + const double dz = d[2]; + + const double r2 = Utils::sqr(dx) + Utils::sqr(dy) + Utils::sqr(dz); + + if (r2 >= Utils::sqr(p.cut)) { + return {}; + } + + const double sig_x2 = Utils::sqr(p.sig_x); + const double sig_y2 = Utils::sqr(p.sig_y); + const double sig_z2 = Utils::sqr(p.sig_z); + + const double A = 0.5 * (Utils::sqr(dx) / sig_x2 + Utils::sqr(dy) / sig_y2 + + Utils::sqr(dz) / sig_z2); + + const double U = p.eps * std::exp(-A); + + // F = -∇U + // With d = r_j - r_i, check sign against ESPResSo's force convention. + return {U * dx / sig_x2, U * dy / sig_y2, U * dz / sig_z2}; +} + +#endif /* ifdef ESPRESSO_GAUSSIAN_ANISO */ +#endif /* GAUSSIAN_ANISO_H */ diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index d98cf0263e..6af3556f7d 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -89,6 +89,11 @@ recalc_maximal_cutoff(IA_parameters const &data, data.gaussian.max_cutoff()); #endif +#ifdef ESPRESSO_GAUSSIAN_ANISO + consider(max_cut_current, mask, PairPotential::GaussianAniso, + data.gaussian_aniso.max_cutoff()); +#endif + #ifdef ESPRESSO_BMHTF_NACL consider(max_cut_current, mask, PairPotential::BMHTF, data.bmhtf.max_cutoff()); diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index d09379aee4..acea9017e4 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -134,6 +134,21 @@ struct Gaussian_Parameters { double max_cutoff() const { return cut; } }; +struct GaussianAniso_Parameters { + double eps = 0.0; + double sig_x = 1.0; + double sig_y = 1.0; + double sig_z = 1.0; + double cut = inactive_cutoff; + + GaussianAniso_Parameters() = default; + + GaussianAniso_Parameters(double eps, double sig_x, double sig_y, double sig_z, + double cutoff); + + double max_cutoff() const { return cut; } +}; + /** BMHTF NaCl potential */ struct BMHTF_Parameters { double A = 0.0; @@ -289,6 +304,7 @@ enum class PairPotential : unsigned { SmoothStep, Hertzian, Gaussian, + GaussianAniso, BMHTF, Buckingham, Morse, @@ -343,6 +359,10 @@ struct IA_parameters { Gaussian_Parameters gaussian; #endif +#ifdef ESPRESSO_GAUSSIAN_ANISO + GaussianAniso_Parameters gaussian_aniso; +#endif + #ifdef ESPRESSO_BMHTF_NACL BMHTF_Parameters bmhtf; #endif diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index fd341a1f85..fd7cf34973 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -706,6 +706,40 @@ def default_params(self): return {} +@script_interface_register +class GaussianAnisoInteraction(NonBondedInteraction): + """Anisotropic Gaussian interaction. + + Methods + ------- + set_params() + Set new parameters for the interaction. + + Parameters + ---------- + eps : :obj:`float` + Amplitude of the Gaussian interaction. + sig_x : :obj:`float` + Width of the Gaussian interaction along x. + sig_y : :obj:`float` + Width of the Gaussian interaction along y. + sig_z : :obj:`float` + Width of the Gaussian interaction along z. + cutoff : :obj:`float` + Radial cutoff distance of the interaction. + + """ + + _so_name = "Interactions::InteractionGaussianAniso" + _so_feature = "GAUSSIAN_ANISO" + + def default_params(self): + """Python dictionary of default parameters. + + """ + return {} + + @script_interface_register class TholeInteraction(NonBondedInteraction): """Thole interaction. diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index edd4fd13c0..7f867efef0 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/src/script_interface/interactions/NonBondedInteraction.hpp @@ -400,6 +400,34 @@ class InteractionGaussian }; #endif // ESPRESSO_GAUSSIAN +#ifdef ESPRESSO_GAUSSIAN_ANISO +class InteractionGaussianAniso + : public InteractionPotentialInterface<::GaussianAniso_Parameters> { +protected: + CoreInteraction IA_parameters::*get_ptr_offset() const override { + return &::IA_parameters::gaussian_aniso; + } + +public: + InteractionGaussianAniso() { + add_parameters({ + make_autoparameter(&CoreInteraction::eps, "eps"), + make_autoparameter(&CoreInteraction::sig_x, "sig_x"), + make_autoparameter(&CoreInteraction::sig_y, "sig_y"), + make_autoparameter(&CoreInteraction::sig_z, "sig_z"), + make_autoparameter(&CoreInteraction::cut, "cutoff"), + }); + } + +private: + void make_new_instance(VariantMap const ¶ms) override { + m_handle = make_shared_from_args( + params, "eps", "sig_x", "sig_y", "sig_z", "cutoff"); + } +}; +#endif // ESPRESSO_GAUSSIAN_ANISO + #ifdef ESPRESSO_BMHTF_NACL class InteractionBMHTF : public InteractionPotentialInterface<::BMHTF_Parameters> { @@ -737,6 +765,9 @@ class NonBondedInteractionHandle #ifdef ESPRESSO_GAUSSIAN std::shared_ptr m_gaussian; #endif +#ifdef ESPRESSO_GAUSSIAN_ANISO + std::shared_ptr m_gaussian_aniso; +#endif #ifdef ESPRESSO_BMHTF_NACL std::shared_ptr m_bmhtf; #endif @@ -856,6 +887,10 @@ class NonBondedInteractionHandle #ifdef ESPRESSO_GAUSSIAN fun(m_gaussian, "gaussian", "Interactions::InteractionGaussian"); #endif +#ifdef ESPRESSO_GAUSSIAN_ANISO + fun(m_gaussian_aniso, "gaussian_aniso", + "Interactions::InteractionGaussianAniso"); +#endif #ifdef ESPRESSO_BMHTF_NACL fun(m_bmhtf, "bmhtf", "Interactions::InteractionBMHTF"); #endif diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index 5d6e4a22de..c6d302e1a8 100644 --- a/src/script_interface/interactions/initialize.cpp +++ b/src/script_interface/interactions/initialize.cpp @@ -77,6 +77,10 @@ void initialize(Utils::Factory *om) { #ifdef ESPRESSO_GAUSSIAN om->register_new("Interactions::InteractionGaussian"); #endif +#ifdef ESPRESSO_GAUSSIAN_ANISO + om->register_new( + "Interactions::InteractionGaussianAniso"); +#endif #ifdef ESPRESSO_BMHTF_NACL om->register_new("Interactions::InteractionBMHTF"); #endif diff --git a/test_codes/test_gauss_aniso.py b/test_codes/test_gauss_aniso.py new file mode 100644 index 0000000000..c8bd1ff126 --- /dev/null +++ b/test_codes/test_gauss_aniso.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jun 10 17:00:50 2026 + +@author: abhinav +""" + +import sys +# sys.path.insert(0, "/home/abhinav/espresso_github/build/src/python") + +from pathlib import Path +# /home/abku051h/espresso_int if script is in codes/ +ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(ROOT / "build/src/python")) + +import numpy as np +import espressomd + +# ------------------------------------------------------------ +# Two-particle test for anisotropic Gaussian interaction +# ------------------------------------------------------------ + +system = espressomd.System(box_l=[10.0, 10.0, 10.0]) +system.time_step = 0.01 +system.cell_system.skin = 0.4 + +# Clear particles/interactions if needed +system.part.clear() +system.non_bonded_inter.reset() + +# Parameters +eps = -2.0 +sig_x = 0.5 +sig_y = 1.0 +sig_z = 2.0 +cutoff = 1.0 + +# Particle positions +pos0 = np.array([1.0, 1.0, 1.0]) +pos1 = np.array([1.2, 1.0, 1.0]) + +d = pos1 - pos0 +dx, dy, dz = d +r = np.linalg.norm(d) + +# Add particles +p0 = system.part.add(pos=pos0, type=0) +p1 = system.part.add(pos=pos1, type=0) + +# Set interaction +system.non_bonded_inter[0, 0].gaussian_aniso.set_params( + eps=eps, + sig_x=sig_x, + sig_y=sig_y, + sig_z=sig_z, + cutoff=cutoff, +) + + +print("features:", espressomd.features()) +print("params:", system.non_bonded_inter[0, 0].gaussian_aniso.get_params()) +print("p0 type:", p0.type, "p1 type:", p1.type) + +# Compute forces and energy without integrating +# system.integrator.run(0) +system.integrator.run(0, recalc_forces=True) + +# Analytic energy +expected_energy = eps * np.exp( + -0.5 * ( + dx**2 / sig_x**2 + + dy**2 / sig_y**2 + + dz**2 / sig_z**2 + ) +) + +expected_force_on_0 = -expected_energy * np.array([ + dx / sig_x**2, + dy / sig_y**2, + dz / sig_z**2, +]) + +energy_dict = system.analysis.energy() +measured_energy = energy_dict["non_bonded"] + +f0 = np.array(p0.f) +f1 = np.array(p1.f) + +print("------------------------------------------------------------") +print("Two-particle anisotropic Gaussian test") +print("------------------------------------------------------------") +print(f"d = {d}") +print(f"r = {r}") +print() +print("Energy:") +print(f" expected = {expected_energy:.16e}") +print(f" measured = {measured_energy:.16e}") +print(f" diff = {measured_energy - expected_energy:.16e}") +print() +print("Force on particle 0:") +print(f" expected = {expected_force_on_0}") +print(f" measured = {f0}") +print(f" diff = {f0 - expected_force_on_0}") +print() +print("Force on particle 1:") +print(f" expected = {-expected_force_on_0}") +print(f" measured = {f1}") +print(f" diff = {f1 + expected_force_on_0}") +print() +print("Newton's third law check:") +print(f" f0 + f1 = {f0 + f1}") +print() + +# Assertions +np.testing.assert_allclose( + measured_energy, expected_energy, rtol=1e-12, atol=1e-12) +np.testing.assert_allclose(f0, expected_force_on_0, rtol=1e-12, atol=1e-12) +np.testing.assert_allclose(f1, -expected_force_on_0, rtol=1e-12, atol=1e-12) +np.testing.assert_allclose(f0 + f1, np.zeros(3), rtol=1e-12, atol=1e-12) + +print("PASS: energy and force match analytic anisotropic Gaussian.")