Permalink
Browse files

Merge pull request #1484 from pkreissl/test_interaction_potentials

Test interaction potentials
  • Loading branch information...
fweik committed Oct 11, 2017
2 parents 3f30deb + d19af93 commit e99423cf44eb015c5ef04f5ccb7b63ec544c9b22
View
@@ -232,9 +232,9 @@ section :ref:`Capping the force during warmup`
The optional ``LJGEN_SOFTCORE`` feature activates a softcore version of
the potential, where the following transformations apply:
:math:`\epsilon \rightarrow \lambda \epsilon` and
:math:`r-r_\mathrm{off} \rightarrow \sqrt{(r-r_\mathrm{off})^2 -
(1-\lambda) \delta \sigma^2}`. allows to tune the strength of the
interaction, while varies how smoothly the potential goes to zero as
:math:`r-r_\mathrm{off} \rightarrow \sqrt{(r-r_\mathrm{off})^2 +
(1-\lambda) \delta \sigma^2}`. :math:`\lambda` allows to tune the strength of the
interaction, while :math:`\delta` varies how smoothly the potential goes to zero as
:math:`\lambda\rightarrow 0`. Such a feature allows one to perform
alchemical transformations, where a group of atoms can be slowly turned
on/off during a simulation.
@@ -218,8 +218,8 @@ typedef struct {
double LJGEN_shift;
double LJGEN_offset;
double LJGEN_capradius;
int LJGEN_a1;
int LJGEN_a2;
double LJGEN_a1;
double LJGEN_a2;
double LJGEN_b1;
double LJGEN_b2;
double LJGEN_lambda;
View
@@ -40,7 +40,7 @@
int ljgen_set_params(int part_type_a, int part_type_b,
double eps, double sig, double cut,
double shift, double offset,
int a1, int a2, double b1, double b2,
double a1, double a2, double b1, double b2,
double cap_radius
#ifdef LJGEN_SOFTCORE
, double lambda, double softrad
View
@@ -41,7 +41,7 @@
#include "utils.hpp"
int ljgen_set_params(int part_type_a, int part_type_b, double eps, double sig,
double cut, double shift, double offset, int a1, int a2,
double cut, double shift, double offset, double a1, double a2,
double b1, double b2, double cap_radius
#ifdef LJGEN_SOFTCORE
,
@@ -59,14 +59,14 @@ inline void add_ljgen_pair_force(const Particle *const p1,
double r_off, frac, fac = 0.0;
r_off = dist - ia_params->LJGEN_offset;
#ifdef LJGEN_SOFTCORE
r_off *= r_off;
#ifdef LJGEN_SOFTCORE
r_off += pow(ia_params->LJGEN_sig, 2) * (1.0 - ia_params->LJGEN_lambda) *
ia_params->LJGEN_softrad;
#endif
/* Taking a square root is not optimal, but we can't prevent the user from
using an odd m, n coefficient. */
r_off = sqrt(r_off);
#endif
/* normal case: resulting force/energy smaller than capping. */
if ((ia_params->LJGEN_capradius == 0) ||
(r_off > ia_params->LJGEN_capradius)) {
@@ -149,14 +149,14 @@ inline double ljgen_pair_energy(Particle *p1, Particle *p2,
if ((dist < ia_params->LJGEN_cut + ia_params->LJGEN_offset)) {
r_off = dist - ia_params->LJGEN_offset;
#ifdef LJGEN_SOFTCORE
r_off *= r_off;
#ifdef LJGEN_SOFTCORE
r_off += pow(ia_params->LJGEN_sig, 2) * (1.0 - ia_params->LJGEN_lambda) *
ia_params->LJGEN_softrad;
#endif
/* Taking a square root is not optimal, but we can't prevent the user from
using an odd m, n coefficient. */
r_off = sqrt(r_off);
#endif
/* normal case: resulting force/energy smaller than capping. */
if (r_off > ia_params->LJGEN_capradius) {
frac = ia_params->LJGEN_sig / r_off;
@@ -40,8 +40,8 @@ cdef extern from "interaction_data.hpp":
double LJGEN_shift
double LJGEN_offset
double LJGEN_capradius
int LJGEN_a1
int LJGEN_a2
double LJGEN_a1
double LJGEN_a2
double LJGEN_b1
double LJGEN_b2
double LJGEN_lambda
@@ -88,14 +88,14 @@ cdef extern from "ljgen.hpp":
cdef int ljgen_set_params(int part_type_a, int part_type_b,
double eps, double sig, double cut,
double shift, double offset,
int a1, int a2, double b1, double b2,
double a1, double a2, double b1, double b2,
double cap_radius,
double genlj_lambda, double softrad)
ELSE:
cdef int ljgen_set_params(int part_type_a, int part_type_b,
double eps, double sig, double cut,
double shift, double offset,
int a1, int a2, double b1, double b2,
double a1, double a2, double b1, double b2,
double cap_radius)
@@ -374,7 +374,7 @@ IF LENNARD_JONES_GENERIC == 1:
"e2": ia_params.LJGEN_a2,
"b1": ia_params.LJGEN_b1,
"b2": ia_params.LJGEN_b2,
"lambda": ia_params.LJGEN_lambda,
"lam": ia_params.LJGEN_lambda,
"delta": ia_params.LJGEN_softrad
}
@@ -399,7 +399,7 @@ IF LENNARD_JONES_GENERIC == 1:
self._params["b1"],
self._params["b2"],
0.0,
self._params["lambda"],
self._params["lam"],
self._params["delta"]):
raise Exception(
"Could not set Generic Lennard Jones parameters")
@@ -430,7 +430,7 @@ IF LENNARD_JONES_GENERIC == 1:
"b1": 0.,
"b2": 0.,
"delta": 0.,
"lambda": 0.}
"lam": 0.}
def type_name(self):
return "GenericLennardJones"
@@ -462,14 +462,14 @@ IF LENNARD_JONES_GENERIC == 1:
delta : float, optional
LJGEN_SOFTCORE parameter. Allows control over how smoothly
the potential drops to zero as lambda approaches zero.
lambda : float, optional
LJGEN_SOFTCORE parameter. Tune the strength of the
lam : float, optional
LJGEN_SOFTCORE parameter lambda. Tune the strength of the
interaction.
"""
super(GenericLennardJonesInteraction, self).set_params(**kwargs)
def valid_keys(self):
return "epsilon", "sigma", "cutoff", "shift", "offset", "e1", "e2", "b1", "b2", "delta", "lambda"
return "epsilon", "sigma", "cutoff", "shift", "offset", "e1", "e2", "b1", "b2", "delta", "lam"
def required_keys(self):
return "epsilon", "sigma", "cutoff", "shift", "offset", "e1", "e2", "b1", "b2"
View
@@ -23,7 +23,6 @@ function(PYTHON_TEST)
set(python_tests ${python_tests} ${TEST_FILE} PARENT_SCOPE)
endfunction(PYTHON_TEST)
python_test(FILE bondedInteractions.py MAX_NUM_PROC 4)
python_test(FILE cellsystem.py MAX_NUM_PROC 4)
python_test(FILE constraint_homogeneous_magnetic_field.py MAX_NUM_PROC 4)
python_test(FILE constraint_shape_based.py MAX_NUM_PROC 4)
@@ -39,7 +38,11 @@ python_test(FILE ewald_gpu.py MAX_NUM_PROC 4)
python_test(FILE icc.py MAX_NUM_PROC 4)
python_test(FILE magnetostaticInteractions.py MAX_NUM_PROC 1)
python_test(FILE mass-and-rinertia_per_particle.py MAX_NUM_PROC 1)
python_test(FILE nonBondedInteractions.py MAX_NUM_PROC 4)
python_test(FILE interactions_bond_angle.py MAX_NUM_PROC 4)
python_test(FILE interactions_bonded_interface.py MAX_NUM_PROC 4)
python_test(FILE interactions_bonded.py MAX_NUM_PROC 4)
python_test(FILE interactions_non-bonded_interface.py MAX_NUM_PROC 4)
python_test(FILE interactions_non-bonded.py MAX_NUM_PROC 4)
python_test(FILE observables.py MAX_NUM_PROC 4)
python_test(FILE p3m_gpu.py MAX_NUM_PROC 4)
python_test(FILE particle.py MAX_NUM_PROC 4)
@@ -0,0 +1,103 @@
#
# Copyright (C) 2013,2014,2015,2016 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/>.
#
from __future__ import print_function
import espressomd
import numpy
import unittest as ut
class InteractionsNonBondedTest(ut.TestCase):
system = espressomd.System()
box_l = 10.
start_pos = numpy.random.rand(3) * box_l
axis = numpy.random.rand(3)
axis /= numpy.linalg.norm(axis)
rel_pos = numpy.cross(numpy.random.rand(3), axis)
rel_pos /= numpy.linalg.norm(rel_pos)
def setUp(self):
self.system.box_l = [self.box_l] * 3
self.system.cell_system.skin = 0.4
self.system.time_step = 1.
self.system.part.add(id=0, pos=self.start_pos, type=0)
self.system.part.add(id=1, pos=self.start_pos + self.rel_pos, type=0)
self.system.part.add(id=2, pos=self.start_pos + self.rel_pos, type=0)
def tearDown(self):
self.system.part.clear()
def rotate_vector(self, v, k, phi):
"""Rotates vector v around unit vector k by angle phi.
Uses Rodrigues' rotation formula."""
vrot = v * numpy.cos(phi) + numpy.cross(k, v) * \
numpy.sin(phi) + k * numpy.dot(k, v) * (1 - numpy.cos(phi))
return vrot
# Required, since assertAlmostEqual does NOT check significant places
def assertFractionAlmostEqual(self, a, b, places=10):
if abs(b) < 1E-8:
self.assertAlmostEqual(a, b)
else:
self.assertAlmostEqual(a / b, 1., places=4)
def assertItemsFractionAlmostEqual(self, a, b):
for i, ai in enumerate(a):
self.assertFractionAlmostEqual(ai, b[i])
# Analytical Expression
def angle_harmonic_potential(self, phi, bend=1.0, phi0=numpy.pi):
return 0.5 * bend * numpy.power(phi - phi0, 2)
# Test Angle Harmonic Potential
@ut.skipIf(not espressomd.has_features(["BOND_ANGLE"]),
"Features not available, skipping test!")
def test_angle_harmonic(self):
ah_bend = 1.
ah_phi0 = 0.4327 * numpy.pi
angle_harmonic = espressomd.interactions.Angle_Harmonic(
bend=ah_bend, phi0=ah_phi0)
self.system.bonded_inter.add(angle_harmonic)
self.system.part[0].add_bond((angle_harmonic, 1, 2))
N = 111
d_phi = numpy.pi / N
for i in range(N):
self.system.part[2].pos = self.start_pos + \
self.rotate_vector(self.rel_pos, self.axis, i * d_phi)
self.system.integrator.run(recalc_forces=True, steps=0)
# Calculate energies
E_sim = self.system.analysis.energy()["bonded"]
E_ref = self.angle_harmonic_potential(
phi=i * d_phi, bend=ah_bend, phi0=ah_phi0)
# Check that energies match
self.assertFractionAlmostEqual(E_sim, E_ref)
if __name__ == '__main__':
print("Features: ", espressomd.features())
ut.main()
Oops, something went wrong.

0 comments on commit e99423c

Please sign in to comment.