Merge pull request #2292 from fweik/wca

Weeks-Chandler-Andersen Potential

doc/sphinx/inter_non-bonded.rst

@@ -216,6 +216,40 @@ interaction, while :math:`\delta` varies how smoothly the potential goes to zero
 alchemical transformations, where a group of atoms can be slowly turned
 on/off during a simulation.
 
+.. _Weeks-Chandler-Andersen interaction:
+
+Weeks-Chandler-Andersen interaction
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. note::
+    Feature ``WCA`` required.
+
+
+The interface for the Weeks-Chandler-Andersen interactions is implemented in
+:class:`espressomd.interactions.WCAInteraction`. They
+are configured via the syntax::
+
+    system.non_bonded_inter[type1, type2].wca.set_params(**kwargs)
+
+This command defines a Weeks-Chandler-Andersen interaction between particles of the
+types ``type1`` and ``type2``. The potential is defined by
+
+.. math::
+
+   \label{eq:wca}
+     V_\mathrm{WCA}(r) =
+       \begin{cases}
+         4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12}
+         - \left(\frac{\sigma}{r}\right)^6 + \frac{1}{4} \right]
+         & \mathrm{if~} r < \sigma 2^{\frac{1}{6}}\\
+         0
+         & \mathrm{otherwise}
+       \end{cases}.
+
+The net force on a particle can be capped by using
+force capping ``system.non_bonded_inter.set_force_cap(max)``, see
+section :ref:`Capping the force during warmup`
+
 .. _Lennard-Jones cosine interaction:
 
 Lennard-Jones cosine interaction

doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb

@@ -18,7 +18,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "One of the typical applications of **ESPResSo** is the simulation of polymer chains with a bead-spring-model. For this we need a repulsive interaction between all beads, for which one usually takes a shifted and truncated Lennard-Jones (so called Weeks-Chandler-Anderson) interaction, and additionally a bonded interaction between adjacent beads to hold the polymer together. You have already learned that the command"
+    "One of the typical applications of **ESPResSo** is the simulation of polymer chains with a bead-spring-model. For this we need a repulsive interaction between all beads, for which one usually takes a shifted and truncated Lennard-Jones (so called Weeks-Chandler-Andersen) interaction, and additionally a bonded interaction between adjacent beads to hold the polymer together. You have already learned that the command"
    ]
   },
   {

src/core/cuda_interface.hpp

@@ -94,8 +94,6 @@ struct CUDA_particle_data {
   float mass;
 #endif
 
-  unsigned int fixed;
-
 #ifdef VIRTUAL_SITES
   bool is_virtual;
 #endif

src/core/energy_inline.hpp

@@ -58,6 +58,7 @@
 #include "nonbonded_interactions/soft_sphere.hpp"
 #include "nonbonded_interactions/steppot.hpp"
 #include "nonbonded_interactions/thole.hpp"
+#include "nonbonded_interactions/wca.hpp"
 #ifdef ELECTROSTATICS
 #include "bonded_interactions/bonded_coulomb.hpp"
 #endif
@@ -100,6 +101,10 @@ inline double calc_non_bonded_pair_energy(const Particle *p1,
   /* Lennard-Jones */
   ret += lj_pair_energy(p1, p2, ia_params, d, dist);
 #endif
+#ifdef WCA
+  /* WCA */
+  ret += wca_pair_energy(p1, p2, ia_params, d, dist);
+#endif
 
 #ifdef LENNARD_JONES_GENERIC
   /* Generic Lennard-Jones */

src/core/forces_inline.hpp

@@ -61,6 +61,7 @@
 #include "nonbonded_interactions/soft_sphere.hpp"
 #include "nonbonded_interactions/steppot.hpp"
 #include "nonbonded_interactions/thole.hpp"
+#include "nonbonded_interactions/wca.hpp"
 #include "npt.hpp"
 #include "object-in-fluid/affinity.hpp"
 #include "object-in-fluid/membrane_collision.hpp"
@@ -185,6 +186,10 @@ inline void calc_non_bonded_pair_force_parts(
 #ifdef LENNARD_JONES
   add_lj_pair_force(p1, p2, ia_params, d, dist, force);
 #endif
+/* WCA */
+#ifdef WCA
+  add_wca_pair_force(p1, p2, ia_params, d, dist, force);
+#endif
 /* Lennard-Jones generic */
 #ifdef LENNARD_JONES_GENERIC
   add_ljgen_pair_force(p1, p2, ia_params, d, dist, force);

src/core/nonbonded_interactions/CMakeLists.txt

@@ -17,6 +17,7 @@ add_library(nonbonded_interactions SHARED
     soft_sphere.cpp
     steppot.cpp
     thole.cpp
+    wca.cpp
 )
 
 add_dependencies(nonbonded_interactions EspressoConfig)

src/core/nonbonded_interactions/nonbonded_interaction_data.cpp

@@ -255,6 +255,10 @@ static void recalc_maximal_cutoff_nonbonded() {
         max_cut_current = (data->LJ_cut + data->LJ_offset);
 #endif
 
+#ifdef WCA
+      max_cut_current = std::max(max_cut_current, data->WCA_cut);
+#endif
+
 #ifdef DPD
       max_cut_current = std::max(max_cut_current,
                                  std::max(data->dpd_r_cut, data->dpd_tr_cut));

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -120,6 +120,12 @@ struct IA_parameters {
 
 #endif
 
+#ifdef WCA
+  double WCA_eps = 0.0;
+  double WCA_sig = 0.0;
+  double WCA_cut = INACTIVE_CUTOFF;
+#endif
+
   /** flag that tells whether there is any short-ranged interaction,
       i.e. one that contributes to the "nonbonded" section of the
       energy/pressure. Note that even if there is no short-ranged

src/core/nonbonded_interactions/wca.cpp

@@ -0,0 +1,40 @@
+/*
+  Copyright (C) 2018 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.hpp"
+
+#ifdef WCA
+#include "wca.hpp"
+
+#include "communication.hpp"
+
+int wca_set_params(int part_type_a, int part_type_b, double eps, double sig) {
+  IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b);
+
+  data->WCA_eps = eps;
+  data->WCA_sig = sig;
+  data->WCA_cut = sig * std::pow(2., 1. / 6.);
+
+  /* broadcast interaction parameters */
+  mpi_bcast_ia_params(part_type_a, part_type_b);
+
+  return ES_OK;
+}
+
+#endif /* ifdef WCA */

src/core/nonbonded_interactions/wca.hpp

@@ -0,0 +1,59 @@
+/*
+  Copyright (C) 2018 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 <http://www.gnu.org/licenses/>.
+*/
+#ifndef WCA_HPP
+#define WCA_HPP
+
+#include "config.hpp"
+
+#ifdef WCA
+
+#include "nonbonded_interaction_data.hpp"
+#include "particle_data.hpp"
+
+int wca_set_params(int part_type_a, int part_type_b, double eps, double sig);
+
+/** Calculate WCA force between particle p1 and p2 */
+inline void add_wca_pair_force(const Particle *const p1,
+                               const Particle *const p2,
+                               IA_parameters *ia_params, double const d[3],
+                               double dist, double force[3]) {
+  if (dist < ia_params->WCA_cut) {
+    auto const frac2 = Utils::sqr(ia_params->WCA_sig / dist);
+    auto const frac6 = frac2 * frac2 * frac2;
+    auto const fac =
+        48.0 * ia_params->WCA_eps * frac6 * (frac6 - 0.5) / (dist * dist);
+    for (int j = 0; j < 3; j++)
+      force[j] += fac * d[j];
+  }
+}
+
+/** calculate WCA energy between particle p1 and p2. */
+inline double wca_pair_energy(const Particle *p1, const Particle *p2,
+                              const IA_parameters *ia_params, const double d[3],
+                              double dist) {
+  if (dist < ia_params->WCA_cut) {
+    auto const frac2 = Utils::sqr(ia_params->WCA_sig / dist);
+    auto const frac6 = frac2 * frac2 * frac2;
+    return 4.0 * ia_params->WCA_eps * (Utils::sqr(frac6) - frac6 + .25);
+  }
+  return 0.0;
+}
+
+#endif /* ifdef WCA */
+#endif

src/features.def

@@ -85,6 +85,7 @@ SHANCHEN                        requires not ELECTROKINETICS and EXPERIMENTAL_FE
 /* Interaction features */
 TABULATED
 LENNARD_JONES
+WCA
 LJ_WARN_WHEN_CLOSE
 LENNARD_JONES_GENERIC           implies LENNARD_JONES
 LJCOS

src/python/espressomd/electrostatics.pyx

@@ -191,7 +191,7 @@ IF P3M == 1:
                 raise ValueError(
                     "P3M cao has to be an integer between -1 and 7")
 
-            if not (self._params["accuracy"] >= 0):
+            if self._params["tune"] and not (self._params["accuracy"] >= 0):
                 raise ValueError("P3M accuracy has to be positive")
 
             if self._params["epsilon"] == "metallic":

src/python/espressomd/interactions.pxd

@@ -47,6 +47,10 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
         double LJ_offset
         double LJ_min
 
+        double WCA_eps
+        double WCA_sig
+        double WCA_cut
+
         double LJCOS_eps
         double LJCOS_sig
         double LJCOS_cut
@@ -160,6 +164,10 @@ cdef extern from "nonbonded_interactions/lj.hpp":
                                       double eps, double sig, double cut,
                                       double shift, double offset,
                                       double min)
+
+cdef extern from "nonbonded_interactions/wca.hpp":
+    cdef int wca_set_params(int part_type_a, int part_type_b,
+                            double eps, double sig)
 IF LJCOS:
     cdef extern from "nonbonded_interactions/ljcos.hpp":
         cdef int ljcos_set_params(int part_type_a, int part_type_b,

src/python/espressomd/interactions.pyx

@@ -324,8 +324,87 @@ IF LENNARD_JONES == 1:
             """
             return "epsilon", "sigma", "cutoff", "shift"
 
-# Generic Lennard-Jones
+IF WCA == 1:
+    cdef class WCAInteraction(NonBondedInteraction):
 
+        def validate_params(self):
+            """Check that parameters are valid.
+
+            Raises
+            ------
+            ValueError
+                If not true.
+            """
+            if self._params["epsilon"] < 0:
+                raise ValueError("WCA eps has to be >=0")
+            if self._params["sigma"] < 0:
+                raise ValueError("WCA sigma has to be >=0")
+            return True
+
+        def _get_params_from_es_core(self):
+            cdef IA_parameters * ia_params
+            ia_params = get_ia_param_safe(
+                self._part_types[0],
+                self._part_types[1])
+            return {
+                "epsilon": ia_params.WCA_eps,
+                "sigma": ia_params.WCA_sig,
+                "cutoff": ia_params.WCA_cut}
+
+        def is_active(self):
+            """Check if interaction is active.
+
+            """
+            return (self._params["epsilon"] > 0)
+
+        def set_params(self, **kwargs):
+            """ Set parameters for the WCA interaction.
+
+            Parameters
+            ----------
+
+            epsilon : :obj:`float`
+                      The magnitude of the interaction.
+            sigma : :obj:`float`
+                    Determines the interaction length scale.
+
+            """
+            super(WCAInteraction, self).set_params(**kwargs)
+
+        def _set_params_in_es_core(self):
+            if wca_set_params(
+                    self._part_types[0], self._part_types[1],
+                    self._params["epsilon"],
+                    self._params["sigma"]):
+                raise Exception("Could not set WCA parameters")
+
+        def default_params(self):
+            """Python dictionary of default parameters.
+
+            """
+            return {
+                "epsilon": 0.,
+                "sigma": 0.}
+
+        def type_name(self):
+            """Name of interaction type.
+
+            """
+            return "WCA"
+
+        def valid_keys(self):
+            """All parameters that can be set.
+
+            """
+            return "epsilon", "sigma"
+
+        def required_keys(self):
+            """Parameters that have to be set.
+
+            """
+            return "epsilon", "sigma"
+
+# Generic Lennard-Jones
 IF LENNARD_JONES_GENERIC == 1:
 
     cdef class GenericLennardJonesInteraction(NonBondedInteraction):
@@ -1672,6 +1751,8 @@ class NonBondedInteractionHandle(object):
         # Here, add one line for each nonbonded ia
         IF LENNARD_JONES:
             self.lennard_jones = LennardJonesInteraction(_type1, _type2)
+        IF WCA:
+            self.wca = WCAInteraction(_type1, _type2)
         IF MEMBRANE_COLLISION:
             self.membrane_collision = MembraneCollisionInteraction(
                 _type1, _type2)