Built new custom gaussian potential which is anisotropic

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

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

src/config/features.def

@@ -72,6 +72,7 @@ LJGEN_SOFTCORE                  implies LENNARD_JONES_GENERIC
 SMOOTH_STEP
 HERTZIAN
 GAUSSIAN
+GAUSSIAN_ANISO
 BMHTF_NACL
 MORSE
 BUCKINGHAM

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)) {

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,

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)) {

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

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

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 <http://www.gnu.org/licenses/>.
+ */
+/** \file
+ *
+ *  Implementation of \ref gaussian_aniso.hpp
+ */
+
+#include "gaussian_aniso.hpp"
+
+#ifdef ESPRESSO_GAUSSIAN_ANISO
+#include "nonbonded_interaction_data.hpp"
+
+#include <stdexcept>
+
+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

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 <http://www.gnu.org/licenses/>.
+ */
+#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 <cmath>
+#include <utils/Vector.hpp>
+#include <utils/math/sqr.hpp>
+
+#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 */