Permalink
Browse files

core: Explicit WCA potential.

  • Loading branch information...
fweik committed Oct 2, 2018
1 parent c7e6d9b commit 521515792bd1249e2aa5dde693294e78a802c48e
@@ -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-Anderson interaction:
Weeks-Chandler-Anderson interaction
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. note::
Feature ``WCA`` required.
The interface for the Weeks-Chandler-Anderson 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-Anderson 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
@@ -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 */
@@ -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);
@@ -17,6 +17,7 @@ add_library(nonbonded_interactions SHARED
soft_sphere.cpp
steppot.cpp
thole.cpp
wca.cpp
)
add_dependencies(nonbonded_interactions EspressoConfig)
@@ -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));
@@ -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
@@ -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 */
@@ -0,0 +1,60 @@
/*
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
View
@@ -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
@@ -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,
@@ -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)
Oops, something went wrong.

0 comments on commit 5215157

Please sign in to comment.