Skip to content

Commit

Permalink
Merge #2718
Browse files Browse the repository at this point in the history 
2718: 2684 coulomb sr test r=fweik a=christophlohrmann

Fixes #2684 and also the way charges and the coulomb prefactor are used in Coulomb energy/force pair calculations 

Description of changes:
 - 


PR Checklist
------------
 - [x] Tests?
   - [x] Interface
   - [ ] Core 
 - [x] Docs?


Co-authored-by: Christoph Lohrmann <clohrmann@icp.uni-stuttgart.de>
Co-authored-by: Florian Weik <florianweik@gmail.com>
Co-authored-by: Jean-Noel Grad <jgrad@icp.uni-stuttgart.de>
  • Loading branch information
@bors @christophlohrmann
@fweik @jngrad
4 people committed Apr 12, 2019
2 parents 2007aaa + 742cc6b commit 51e344c
Show file tree
Hide file tree
Showing 20 changed files with 175 additions and 161 deletions.
2 changes: 1 addition & 1 deletion src/core/CMakeLists.txt
View file
Expand Up @@ -58,7 +58,7 @@ set(EspressoCore_SRC
bonded_interactions/angle_cossquare.cpp
bonded_interactions/angle_harmonic.cpp
bonded_interactions/bonded_coulomb.cpp
bonded_interactions/bonded_coulomb_p3m_sr.cpp
bonded_interactions/bonded_coulomb_sr.cpp
bonded_interactions/bonded_interaction_data.cpp
bonded_interactions/bonded_tab.cpp
bonded_interactions/dihedral.cpp
Expand Down
12 changes: 6 additions & 6 deletions ...ed_interactions/bonded_coulomb_p3m_sr.cpp → ...bonded_interactions/bonded_coulomb_sr.cpp
View file
Expand Up @@ -20,21 +20,21 @@
*/
/** \file
*
* Implementation of \ref bonded_coulomb_p3m_sr.hpp
* Implementation of \ref bonded_coulomb_sr.hpp
*/
#include "bonded_coulomb_p3m_sr.hpp"
#include "bonded_coulomb_sr.hpp"
#include "communication.hpp"

#ifdef P3M
#ifdef ELECTROSTATICS

int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2) {
int bonded_coulomb_sr_set_params(int bond_type, double q1q2) {
if (bond_type < 0)
return ES_ERROR;

make_bond_type_exist(bond_type);

bonded_ia_params[bond_type].p.bonded_coulomb_p3m_sr.q1q2 = q1q2;
bonded_ia_params[bond_type].type = BONDED_IA_BONDED_COULOMB_P3M_SR;
bonded_ia_params[bond_type].p.bonded_coulomb_sr.q1q2 = q1q2;
bonded_ia_params[bond_type].type = BONDED_IA_BONDED_COULOMB_SR;
bonded_ia_params[bond_type].num = 1;

/* broadcast interaction parameters */
Expand Down
56 changes: 29 additions & 27 deletions ...ed_interactions/bonded_coulomb_p3m_sr.hpp → ...bonded_interactions/bonded_coulomb_sr.hpp
View file
Expand Up @@ -18,24 +18,24 @@
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 _BONDED_COULOMB_P3M_SR_HPP
#define _BONDED_COULOMB_P3M_SR_HPP
#ifndef _BONDED_COULOMB_SR_HPP
#define _BONDED_COULOMB_SR_HPP
/** \file
* Routines to calculate the BONDED_COULOMB_P3M_SR Energy or/and
* BONDED_COULOMB_P3M_SR force for a particle pair. This is only the shortrange
* part of P3M and first used to subtract certain intramolecular interactions in
* combination with Thole damping \ref forces.cpp
* Routines to calculate the BONDED_COULOMB_SR Energy or/and
* BONDED_COULOMB_SR force for a particle pair. This is only the shortrange
* part of any coulomb interaction and first used to subtract certain
* intramolecular interactions in combination with Thole damping \ref forces.cpp
*/

/************************************************************/

#include "config.hpp"

#ifdef P3M
#ifdef ELECTROSTATICS

#include "bonded_interaction_data.hpp"
#include "debug.hpp"
#include "electrostatics_magnetostatics/p3m.hpp"
#include "electrostatics_magnetostatics/coulomb_inline.hpp"
#include "particle_data.hpp"
#include "utils.hpp"

Expand All @@ -44,9 +44,9 @@
* @retval ES_OK on success
* @retval ES_ERROR on error
*/
int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2);
int bonded_coulomb_sr_set_params(int bond_type, double q1q2);

/** Computes the BONDED_COULOMB_P3M_SR pair force.
/** Computes the BONDED_COULOMB_SR pair force.
* @param[in] p1 First particle.
* @param[in] p2 Second particle.
* @param[in] iaparams Interaction parameters.
Expand All @@ -55,49 +55,51 @@ int bonded_coulomb_p3m_sr_set_params(int bond_type, double q1q2);
* @retval 0
*/
inline int
calc_bonded_coulomb_p3m_sr_pair_force(Particle const *p1, Particle const *p2,
Bonded_ia_parameters const *iaparams,
double const dx[3], double force[3]) {
calc_bonded_coulomb_sr_pair_force(Particle *p1, Particle *p2,
Bonded_ia_parameters const *iaparams,
double dx[3], double force[3]) {
double dist2 = sqrlen(dx);
double dist = sqrt(dist2);
if (dist < p3m.params.r_cut) {
// Set to zero because p3m adds forces
force[0] = force[1] = force[2] = 0.;
if (dist < coulomb_cutoff) {
// TODO ugly workaround
Vector3d forcevec{};

p3m_add_pair_force(iaparams->p.bonded_coulomb_p3m_sr.q1q2, dx, dist2, dist,
force);
Coulomb::calc_pair_force(p1, p2, iaparams->p.bonded_coulomb_sr.q1q2, dx,
dist, dist2, forcevec);
force[0] = forcevec[0];
force[1] = forcevec[1];
force[2] = forcevec[2];

ONEPART_TRACE(if (p1->p.identity == check_id) fprintf(
stderr,
"%d: OPT: BONDED_COULOMB_P3M_SR f = (%.3e,%.3e,%.3e) with part id=%d "
"%d: OPT: BONDED_COULOMB_SR f = (%.3e,%.3e,%.3e) with part id=%d "
"at dist %f\n",
this_node, p1->f.f[0], p1->f.f[1], p1->f.f[2], p2->p.identity, dist2));
ONEPART_TRACE(if (p2->p.identity == check_id) fprintf(
stderr,
"%d: OPT: BONDED_COULOMB_P3M_SR f = (%.3e,%.3e,%.3e) with part id=%d "
"%d: OPT: BONDED_COULOMB_SR f = (%.3e,%.3e,%.3e) with part id=%d "
"at dist %f\n",
this_node, p2->f.f[0], p2->f.f[1], p2->f.f[2], p1->p.identity, dist2));
}

return 0;
}

/** Computes the BONDED_COULOMB_P3M_SR pair energy.
/** Computes the BONDED_COULOMB_SR pair energy.
* @param[in] p1 First particle.
* @param[in] p2 Second particle.
* @param[in] iaparams Interaction parameters.
* @param[in] dx %Distance between the particles.
* @param[out] _energy Energy.
* @retval 0
*/
inline int
bonded_coulomb_p3m_sr_pair_energy(Particle const *p1, Particle const *p2,
Bonded_ia_parameters const *iaparams,
double const dx[3], double *_energy) {
inline int bonded_coulomb_sr_pair_energy(Particle *p1, Particle *p2,
Bonded_ia_parameters const *iaparams,
double dx[3], double *_energy) {
double dist2 = sqrlen(dx);
double dist = sqrt(dist2);
*_energy = p3m_pair_energy(iaparams->p.bonded_coulomb_p3m_sr.q1q2, dist);

*_energy = Coulomb::add_pair_energy(
p1, p2, iaparams->p.bonded_coulomb_sr.q1q2, dx, dist, dist2);
return 0;
}

Expand Down
14 changes: 5 additions & 9 deletions src/core/bonded_interactions/bonded_interaction_data.hpp
View file
Expand Up @@ -23,6 +23,8 @@ enum BondedInteraction {
BONDED_IA_QUARTIC,
/** Type of bonded interaction is a BONDED_COULOMB */
BONDED_IA_BONDED_COULOMB,
/** Type of bonded interaction is a BONDED_COULOMB_SR */
BONDED_IA_BONDED_COULOMB_SR,
/** Type of bonded interaction is a dihedral potential. */
BONDED_IA_DIHEDRAL,
/** Type of tabulated bonded interaction potential,
Expand Down Expand Up @@ -58,8 +60,6 @@ enum BondedInteraction {
BONDED_IA_UMBRELLA,
/** Type of bonded interaction is thermalized distance bond. */
BONDED_IA_THERMALIZED_DIST,
/** Type of bonded interaction is a BONDED_COULOMB_P3M_SR */
BONDED_IA_BONDED_COULOMB_P3M_SR,
};

/** Specify tabulated bonded interactions */
Expand Down Expand Up @@ -156,13 +156,11 @@ struct Bonded_coulomb_bond_parameters {
double prefactor;
};

#ifdef P3M
/** Parameters for Coulomb bond p3m short-range Potential */
struct Bonded_coulomb_p3m_sr_bond_parameters {
/** Parameters for Coulomb bond short-range Potential */
struct Bonded_coulomb_sr_bond_parameters {
/** charge factor */
double q1q2;
};
#endif

/** Parameters for three-body angular potential.
* @note
Expand Down Expand Up @@ -309,6 +307,7 @@ union Bond_parameters {
#endif
Quartic_bond_parameters quartic;
Bonded_coulomb_bond_parameters bonded_coulomb;
Bonded_coulomb_sr_bond_parameters bonded_coulomb_sr;
Angle_bond_parameters angle;
Angle_harmonic_bond_parameters angle_harmonic;
Angle_cosine_bond_parameters angle_cosine;
Expand All @@ -321,9 +320,6 @@ union Bond_parameters {
Umbrella_bond_parameters umbrella;
#endif
Thermalized_bond_parameters thermalized_bond;
#ifdef P3M
Bonded_coulomb_p3m_sr_bond_parameters bonded_coulomb_p3m_sr;
#endif
Subt_lj_bond_parameters subt_lj;
Rigid_bond_parameters rigid_bond;
IBM_Triel_Parameters ibm_triel;
Expand Down
36 changes: 18 additions & 18 deletions src/core/electrostatics_magnetostatics/coulomb_inline.hpp
View file
Expand Up @@ -15,9 +15,9 @@

namespace Coulomb {
// forces_inline
inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
double dist2, Vector3d &force) {
auto const q1q2 = p1->p.q * p2->p.q;
inline void calc_pair_force(Particle *p1, Particle *p2, double const q1q2,
double *d, double dist, double const dist2,
Vector3d &force) {

if (q1q2 != 0) {
Vector3d f{};
Expand Down Expand Up @@ -54,10 +54,10 @@ inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
add_mmm2d_coulomb_pair_force(q1q2, d, dist2, dist, f.data());
break;
case COULOMB_DH:
add_dh_coulomb_pair_force(p1, p2, d, dist, f.data());
add_dh_coulomb_pair_force(q1q2, d, dist, f.data());
break;
case COULOMB_RF:
add_rf_coulomb_pair_force(p1, p2, d, dist, f.data());
add_rf_coulomb_pair_force(q1q2, d, dist, f.data());
break;
#ifdef SCAFACOS
case COULOMB_SCAFACOS:
Expand All @@ -73,8 +73,8 @@ inline void calc_pair_force(Particle *p1, Particle *p2, double *d, double dist,
}

// pressure_inline.hpp
inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
double dist, double dist2,
inline void add_pair_pressure(Particle *p1, Particle *p2, double q1q2,
double *d, double dist, double dist2,
Observable_stat &virials,
Observable_stat &p_tensor) {
switch (coulomb.method) {
Expand All @@ -88,7 +88,7 @@ inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
case COULOMB_DH:
case COULOMB_RF: {
Vector3d force{};
calc_pair_force(p1, p2, d, dist, dist2, force);
calc_pair_force(p1, p2, q1q2, d, dist, dist2, force);

/* Calculate the virial pressure */
for (int k = 0; k < 3; k++) {
Expand All @@ -107,40 +107,40 @@ inline void add_pair_pressure(Particle *p1, Particle *p2, double *d,
}

// energy_inline
inline double add_pair_energy(Particle *p1, Particle *p2, double *d,
double dist, double dist2) {
inline double add_pair_energy(Particle *p1, Particle *p2, double const q1q2,
double *d, double dist, double dist2) {
/* real space Coulomb */
auto E = [&]() {
switch (coulomb.method) {
#ifdef P3M
case COULOMB_P3M_GPU:
case COULOMB_P3M:
return p3m_pair_energy(p1->p.q * p2->p.q, dist);
// TODO some energy functions include the prefactor, some don't
return p3m_pair_energy(q1q2, dist);
case COULOMB_ELC_P3M:
if (elc_params.dielectric_contrast_on) {
return 0.5 * ELC_P3M_dielectric_layers_energy_contribution(p1, p2) +
p3m_pair_energy(p1->p.q * p2->p.q, dist);
p3m_pair_energy(q1q2, dist);
} else {
return p3m_pair_energy(p1->p.q * p2->p.q, dist);
return p3m_pair_energy(q1q2, dist);
}
#endif
#ifdef SCAFACOS
case COULOMB_SCAFACOS:
return Scafacos::pair_energy(p1, p2, dist);
#endif
case COULOMB_DH:
return dh_coulomb_pair_energy(p1, p2, dist);
return dh_coulomb_pair_energy(q1q2, dist);
case COULOMB_RF:
return rf_coulomb_pair_energy(p1, p2, dist);
return rf_coulomb_pair_energy(q1q2, dist);
case COULOMB_MMM1D:
return mmm1d_coulomb_pair_energy(p1, p2, d, dist2, dist);
return mmm1d_coulomb_pair_energy(q1q2, d, dist2, dist);
case COULOMB_MMM2D:
return mmm2d_coulomb_pair_energy(p1->p.q * p2->p.q, d, dist2, dist);
return mmm2d_coulomb_pair_energy(q1q2, d, dist2, dist);
default:
return 0.;
}
}();

return coulomb.prefactor * E;
}
} // namespace Coulomb
Expand Down
19 changes: 8 additions & 11 deletions src/core/electrostatics_magnetostatics/debye_hueckel.hpp
View file
Expand Up @@ -25,7 +25,6 @@
* for a particle pair.
*/
#include "config.hpp"
#include "electrostatics_magnetostatics/coulomb.hpp"

#ifdef ELECTROSTATICS

Expand Down Expand Up @@ -56,32 +55,30 @@ int dh_set_params(double kappa, double r_cut);
@param dist Distance between p1 and p2.
@param force returns the force on particle 1.
*/
inline void add_dh_coulomb_pair_force(Particle *p1, Particle *p2,
double const d[3], double dist,
double force[3]) {
inline void add_dh_coulomb_pair_force(double const q1q2, double const d[3],
double const dist, double force[3]) {
if (dist < dh_params.r_cut) {
double fac;
if (dh_params.kappa > 0.0) {
/* debye hueckel case: */
double kappa_dist = dh_params.kappa * dist;
fac = coulomb.prefactor * p1->p.q * p2->p.q *
(exp(-kappa_dist) / (dist * dist * dist)) * (1.0 + kappa_dist);
fac =
q1q2 * (exp(-kappa_dist) / (dist * dist * dist)) * (1.0 + kappa_dist);
} else {
/* pure Coulomb case: */
fac = coulomb.prefactor * p1->p.q * p2->p.q / (dist * dist * dist);
fac = q1q2 / (dist * dist * dist);
}
for (int j = 0; j < 3; j++)
force[j] += fac * d[j];
}
}

inline double dh_coulomb_pair_energy(Particle *p1, Particle *p2, double dist) {
inline double dh_coulomb_pair_energy(double const q1q2, double const dist) {
if (dist < dh_params.r_cut) {
if (dh_params.kappa > 0.0)
return coulomb.prefactor * p1->p.q * p2->p.q *
exp(-dh_params.kappa * dist) / dist;
return q1q2 * exp(-dh_params.kappa * dist) / dist;

return coulomb.prefactor * p1->p.q * p2->p.q / dist;
return q1q2 / dist;
}
return 0.0;
}
Expand Down
2 changes: 1 addition & 1 deletion src/core/electrostatics_magnetostatics/icc.cpp
View file
Expand Up @@ -70,7 +70,7 @@ inline void add_non_bonded_pair_force_iccp3m(Particle *p1, Particle *p2,
double d[3], double dist,
double dist2) {
Vector3d force{};
Coulomb::calc_pair_force(p1, p2, d, dist, dist2, force);
Coulomb::calc_pair_force(p1, p2, p1->p.q * p2->p.q, d, dist, dist2, force);

p1->f.f += force;
p2->f.f -= force;
Expand Down

0 comments on commit 51e344c

Please sign in to comment.