diff --git a/src/core/bonded_interactions/CMakeLists.txt b/src/core/bonded_interactions/CMakeLists.txt index 5352ed4d079..6412bf026ce 100644 --- a/src/core/bonded_interactions/CMakeLists.txt +++ b/src/core/bonded_interactions/CMakeLists.txt @@ -1,7 +1,6 @@ add_library(bonded_interactions SHARED angle_cosine.cpp angle_cossquare.cpp - angle_dist.cpp angle_harmonic.cpp bonded_coulomb.cpp bonded_coulomb_p3m_sr.cpp diff --git a/src/core/bonded_interactions/angle_dist.cpp b/src/core/bonded_interactions/angle_dist.cpp deleted file mode 100644 index ca37fbad9c2..00000000000 --- a/src/core/bonded_interactions/angle_dist.cpp +++ /dev/null @@ -1,208 +0,0 @@ -/* - Copyright (C) 2010-2018 The ESPResSo project - Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - Max-Planck-Institute for Polymer Research, Theory Group - - 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 . -*/ -/** \file - * - * Implementation of \ref angle_dist.hpp - */ -#include "angle_dist.hpp" - -#ifdef BOND_ANGLEDIST - -#include - -#include "communication.hpp" -#include "constraints.hpp" -#include "constraints/ShapeBasedConstraint.hpp" -#include "grid.hpp" -#include "shapes/Wall.hpp" -#include - -int angledist_set_params(int bond_type, double bend, double phimin, - double distmin, double phimax, double distmax) { - if (bond_type < 0) - return ES_ERROR; - - make_bond_type_exist(bond_type); - - bonded_ia_params[bond_type].p.angledist.bend = bend; - bonded_ia_params[bond_type].p.angledist.phimin = phimin; - bonded_ia_params[bond_type].p.angledist.distmin = distmin; - bonded_ia_params[bond_type].p.angledist.phimax = phimax; - bonded_ia_params[bond_type].p.angledist.distmax = distmax; - bonded_ia_params[bond_type].type = BONDED_IA_ANGLEDIST; - bonded_ia_params[bond_type].num = 2; - - /* broadcast interaction parameters */ - mpi_bcast_ia_params(bond_type, -1); - - return ES_OK; -} - -/** Function to calculate wall distance and phi0(dist). - Called by \ref calc_angledist_force */ -static double calc_angledist_param(Particle *p_mid, Particle *p_left, - Particle *p_right, - Bonded_ia_parameters *iaparams) { - double vec1[3], vec2[3], d2i = 0.0, dist2 = 0.0; - - /* vector from p_left to p_mid */ - get_mi_vector(vec1, p_mid->r.p, p_left->r.p); - const double dist1 = sqrlen(vec1); - const double d1i = 1.0 / sqrt(dist1); - for (int j = 0; j < 3; j++) - vec1[j] *= d1i; - /* vector from p_mid to p_right */ - get_mi_vector(vec2, p_right->r.p, p_mid->r.p); - dist2 = sqrlen(vec2); - d2i = 1.0 / sqrt(dist2); - for (int j = 0; j < 3; j++) - vec2[j] *= d2i; - - const double phimn = iaparams->p.angledist.phimin; - const double distmn = iaparams->p.angledist.distmin; - const double phimx = iaparams->p.angledist.phimax; - const double distmx = iaparams->p.angledist.distmax; - - /* folds coordinates of p_mid into original box */ - Vector3d folded_pos = folded_position(*p_mid); - - double pwdistmin = std::numeric_limits::infinity(); - double pwdistmin_d = 0.0; - - for (auto const &c : Constraints::constraints) { - auto cs = - std::dynamic_pointer_cast(c); - if (cs) { - if (dynamic_cast(&(cs->shape()))) { - const Shapes::Wall *wall = - dynamic_cast(&(cs->shape())); - double dist = -wall->d(); - for (int j = 0; j < 3; j++) { - dist += folded_pos[j] * (wall->n())[j]; - } - - if (dist < pwdistmin) { - pwdistmin = dist; - pwdistmin_d = wall->d(); - } - } - } - } - - /*get phi0(z)*/ - if (pwdistmin <= distmn) { - return phimn; - } else if (pwdistmin >= distmx && - pwdistmin <= box_l[2] - pwdistmin_d - distmx) { - return phimx; - } else { - const double drange = (pwdistmin - distmn) * PI / (distmx - distmn); - return ((cos(drange - PI) + 1.0) * (phimx - phimn)) * 0.5 + phimn; - } -} - -int calc_angledist_force(Particle *p_mid, Particle *p_left, Particle *p_right, - Bonded_ia_parameters *iaparams, double force1[3], - double force2[3]) { - double cosine = 0.0, vec1[3], vec2[3], fac = 0.0; - - /* vector from p_left to p_mid */ - get_mi_vector(vec1, p_mid->r.p, p_left->r.p); - const double dist12 = sqrlen(vec1); - const double d1i = 1.0 / sqrt(dist12); - for (int j = 0; j < 3; j++) - vec1[j] *= d1i; - - /* vector from p_mid to p_right */ - get_mi_vector(vec2, p_right->r.p, p_mid->r.p); - const double dist22 = sqrlen(vec2); - const double d2i = 1.0 / sqrt(dist22); - for (int j = 0; j < 3; j++) - vec2[j] *= d2i; - /* scalar product of vec1 and vec2 */ - cosine = scalar(vec1, vec2); - fac = iaparams->p.angledist.bend; - -#ifdef BOND_ANGLEDIST_HARMONIC - { - const double phi0 = calc_angledist_param(p_mid, p_left, p_right, iaparams); - - if (cosine > TINY_COS_VALUE) - cosine = TINY_COS_VALUE; - if (cosine < -TINY_COS_VALUE) - cosine = -TINY_COS_VALUE; - const double phi = acos(-cosine); - - double sinphi = sin(phi); - if (sinphi < TINY_SIN_VALUE) - sinphi = TINY_SIN_VALUE; - fac *= (phi - phi0) / sinphi; - } -#endif - - for (int j = 0; j < 3; j++) { - double f1 = fac * (cosine * vec1[j] - vec2[j]) * d1i; - double f2 = fac * (cosine * vec2[j] - vec1[j]) * d2i; - - force1[j] = (f1 - f2); - force2[j] = -f1; - } - return 0; -} - -#ifdef BOND_ANGLEDIST_HARMONIC -int angledist_energy(Particle *p_mid, Particle *p_left, Particle *p_right, - Bonded_ia_parameters *iaparams, double *_energy) { - int j; - double cosine = 0.0, d1i, d2i, dist1, dist2; - double vec1[3], vec2[3]; - - /* vector from p_mid to p_left */ - get_mi_vector(vec1, p_mid->r.p, p_left->r.p); - dist1 = sqrlen(vec1); - d1i = 1.0 / sqrt(dist1); - for (j = 0; j < 3; j++) - vec1[j] *= d1i; - /* vector from p_right to p_mid */ - get_mi_vector(vec2, p_right->r.p, p_mid->r.p); - dist2 = sqrlen(vec2); - d2i = 1.0 / sqrt(dist2); - for (j = 0; j < 3; j++) - vec2[j] *= d2i; - /* scalar product of vec1 and vec2 */ - cosine = scalar(vec1, vec2); - if (cosine > TINY_COS_VALUE) - cosine = TINY_COS_VALUE; - if (cosine < -TINY_COS_VALUE) - cosine = -TINY_COS_VALUE; - - { - double phi; - double phi0 = calc_angledist_param(p_mid, p_left, p_right, iaparams); - phi = acos(-cosine); - *_energy = 0.5 * iaparams->p.angledist.bend * Utils::sqr(phi - phi0); - } - - return 0; -} -#endif - -#endif diff --git a/src/core/bonded_interactions/angle_dist.hpp b/src/core/bonded_interactions/angle_dist.hpp deleted file mode 100644 index de85f011a8d..00000000000 --- a/src/core/bonded_interactions/angle_dist.hpp +++ /dev/null @@ -1,62 +0,0 @@ -/* - Copyright (C) 2010-2018 The ESPResSo project - Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - Max-Planck-Institute for Polymer Research, Theory Group - - 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 . -*/ -#ifndef _ANGLEDIST_H -#define _ANGLEDIST_H -/** \file - * Routines to calculate the angle and distance dependent (from a constraint) - * energy or/and and force - * for a particle triple. - * \ref forces.cpp - */ - -#include "config.hpp" - -#ifdef BOND_ANGLEDIST - -#include "bonded_interaction_data.hpp" -#include "particle_data.hpp" -#include "utils.hpp" - -/** set parameters for the angledist potential. The type of the - angledist potential is chosen via config.hpp and cannot be changed - at runtime. -**/ -int angledist_set_params(int bond_type, double bend, double phimin, - double distmin, double phimax, double distmax); - -/// -int calc_angledist_force(Particle *p_mid, Particle *p_left, Particle *p_right, - Bonded_ia_parameters *iaparams, double force1[3], - double force2[3]); - -/** Computes the three body angle interaction energy. - @param p_mid Pointer to second/middle particle. - @param p_left Pointer to first particle. - @param p_right Pointer to third particle. - @param iaparams bond type number of the angle interaction. - @param _energy return energy pointer. - @return 0. -*/ -int angledist_energy(Particle *p_mid, Particle *p_left, Particle *p_right, - Bonded_ia_parameters *iaparams, double *_energy); - -#endif /* BOND_ANGLEDIST */ -#endif /* ANGLEDIST_H */ diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 1813e6fe4ba..32d92523e42 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -28,7 +28,6 @@ #include "bonded_interactions/angle_cosine.hpp" #include "bonded_interactions/angle_cossquare.hpp" -#include "bonded_interactions/angle_dist.hpp" #include "bonded_interactions/angle_harmonic.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "bonded_interactions/bonded_tab.hpp" @@ -394,11 +393,6 @@ inline void add_bonded_energy(Particle *p1) { bond_broken = angle_cossquare_energy(p1, p2, p3, iaparams, &ret); break; #endif -#ifdef BOND_ANGLEDIST - case BONDED_IA_ANGLEDIST: - bond_broken = angledist_energy(p1, p2, p3, iaparams, &ret); - break; -#endif #ifdef TABULATED case BONDED_IA_TABULATED: if (iaparams->num == 2) diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 1e4994d851c..ed887cfb668 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -25,7 +25,6 @@ #include "bonded_interactions/angle_cosine.hpp" #include "bonded_interactions/angle_cossquare.hpp" -#include "bonded_interactions/angle_dist.hpp" #include "bonded_interactions/angle_harmonic.hpp" #include "bonded_interactions/bonded_tab.hpp" #include "bonded_interactions/dihedral.hpp" @@ -617,11 +616,6 @@ inline void add_bonded_force(Particle *p1) { calc_tab_angle_force(p1, p2, p3, iaparams, force, force2); break; #endif -#ifdef BOND_ANGLEDIST - case BONDED_IA_ANGLEDIST: - bond_broken = calc_angledist_force(p1, p2, p3, iaparams, force, force2); - break; -#endif #ifdef IMMERSED_BOUNDARY case BONDED_IA_IBM_TRIEL: bond_broken = IBM_Triel_CalcForce(p1, p2, p3, iaparams); diff --git a/src/features.def b/src/features.def index d92471959b4..ec43dcca0a4 100644 --- a/src/features.def +++ b/src/features.def @@ -116,9 +116,6 @@ MEMBRANE_COLLISION BOND_ANGLE requires not BOND_ANGLE_OLD -BOND_ANGLEDIST -BOND_ANGLEDIST_HARMONIC implies BOND_ANGLEDIST - /* Immersed-Boundary Bayreuth version */ IMMERSED_BOUNDARY IMMERSED_BOUNDARY requires VIRTUAL_SITES