Skip to content

Commit

Permalink
Bond cutoff (#3443)
Browse files Browse the repository at this point in the history 
Fixes #3442

Description of changes:
 - Define and consider cutoff for all bond types.
  • Loading branch information
@kodiakhq
kodiakhq[bot] authored Jan 30, 2020
2 parents 846690f + a757875 commit 780c7fd
Show file tree
Hide file tree
Showing 3 changed files with 127 additions and 61 deletions.
116 changes: 72 additions & 44 deletions src/core/bonded_interactions/bonded_interaction_data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,59 +19,86 @@
#include "bonded_interaction_data.hpp"
#include "communication.hpp"

#include <boost/range/numeric.hpp>
#include <utils/constants.hpp>

std::vector<Bonded_ia_parameters> bonded_ia_params;

auto cutoff(int type, Bond_parameters const &bp) {
switch (type) {
case BONDED_IA_NONE:
return -1.;
case BONDED_IA_FENE:
return bp.fene.cutoff();
case BONDED_IA_HARMONIC:
return bp.harmonic.cutoff();
case BONDED_IA_HARMONIC_DUMBBELL:
return bp.harmonic_dumbbell.cutoff();
case BONDED_IA_QUARTIC:
return bp.quartic.cutoff();
case BONDED_IA_BONDED_COULOMB:
return bp.bonded_coulomb.cutoff();
case BONDED_IA_BONDED_COULOMB_SR:
return bp.bonded_coulomb_sr.cutoff();
case BONDED_IA_DIHEDRAL:
return bp.dihedral.cutoff();
case BONDED_IA_TABULATED_DISTANCE:
case BONDED_IA_TABULATED_ANGLE:
case BONDED_IA_TABULATED_DIHEDRAL:
return bp.tab.cutoff();
case BONDED_IA_SUBT_LJ:
return bp.subt_lj.cutoff();
case BONDED_IA_RIGID_BOND:
return bp.rigid_bond.cutoff();
case BONDED_IA_VIRTUAL_BOND:
return bp.virt.cutoff();
case BONDED_IA_ANGLE_HARMONIC:
return bp.angle_harmonic.cutoff();
case BONDED_IA_ANGLE_COSINE:
return bp.angle_cosine.cutoff();
case BONDED_IA_ANGLE_COSSQUARE:
return bp.angle_cossquare.cutoff();
case BONDED_IA_OIF_LOCAL_FORCES:
return bp.oif_local_forces.cutoff();
case BONDED_IA_OIF_GLOBAL_FORCES:
return bp.oif_global_forces.cutoff();
case BONDED_IA_IBM_TRIEL:
return bp.ibm_triel.cutoff();
case BONDED_IA_IBM_VOLUME_CONSERVATION:
return bp.ibmVolConsParameters.cutoff();
case BONDED_IA_IBM_TRIBEND:
return bp.ibm_tribend.cutoff();
case BONDED_IA_UMBRELLA:
return bp.umbrella.cutoff();
case BONDED_IA_THERMALIZED_DIST:
return bp.thermalized_bond.cutoff();
default:
throw std::runtime_error("Unknown bond type.");
}
}

double recalc_maximal_cutoff_bonded() {
auto max_cut_bonded = -1.;
auto const max_cut_bonded =
boost::accumulate(bonded_ia_params, -1.,
[](auto max_cut, Bonded_ia_parameters const &bond) {
return std::max(max_cut, cutoff(bond.type, bond.p));
});

for (auto const &bonded_ia_param : bonded_ia_params) {
switch (bonded_ia_param.type) {
case BONDED_IA_FENE:
max_cut_bonded =
std::max(max_cut_bonded,
bonded_ia_param.p.fene.r0 + bonded_ia_param.p.fene.drmax);
break;
case BONDED_IA_HARMONIC:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.harmonic.r_cut);
break;
case BONDED_IA_THERMALIZED_DIST:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.thermalized_bond.r_cut);
break;
case BONDED_IA_RIGID_BOND:
max_cut_bonded =
std::max(max_cut_bonded, std::sqrt(bonded_ia_param.p.rigid_bond.d2));
break;
case BONDED_IA_TABULATED_DISTANCE:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.tab.pot->cutoff());
break;
case BONDED_IA_IBM_TRIEL:
max_cut_bonded =
std::max(max_cut_bonded, bonded_ia_param.p.ibm_triel.maxDist);
break;
default:
break;
}
}
/* Check if there are dihedrals */
auto const any_dihedrals = std::any_of(
bonded_ia_params.begin(), bonded_ia_params.end(), [](auto const &bond) {
switch (bond.type) {
case BONDED_IA_DIHEDRAL:
case BONDED_IA_TABULATED_DIHEDRAL:
return true;
default:
return false;
}
});

/* dihedrals: the central particle is indirectly connected to the fourth
* particle via the third particle, so we have to double the cutoff */
for (auto const &bonded_ia_param : bonded_ia_params) {
switch (bonded_ia_param.type) {
case BONDED_IA_DIHEDRAL:
case BONDED_IA_TABULATED_DIHEDRAL:
max_cut_bonded *= 2;
break;
default:
break;
}
}

return max_cut_bonded;
return (any_dihedrals) ? 2 * max_cut_bonded : max_cut_bonded;
}

void make_bond_type_exist(int type) {
Expand All @@ -95,6 +122,7 @@ int virtual_set_params(int bond_type) {

bonded_ia_params[bond_type].type = BONDED_IA_VIRTUAL_BOND;
bonded_ia_params[bond_type].num = 1;
bonded_ia_params[bond_type].p.virt = VirtualBond_Parameters{};

/* broadcast interaction parameters */
mpi_bcast_ia_params(bond_type, -1);
Expand Down
71 changes: 54 additions & 17 deletions src/core/bonded_interactions/bonded_interaction_data.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,6 @@ enum BondedInteraction {
BONDED_IA_OIF_LOCAL_FORCES,
/** Type of bonded interaction: OIF global forces. */
BONDED_IA_OIF_GLOBAL_FORCES,
/** Type of bonded interaction: determining outward direction of OIF membrane
* (not associated to a parameter struct).
*/
BONDED_IA_OIF_OUT_DIRECTION,
/** Type of bonded interaction is a wall repulsion (immersed boundary). */
BONDED_IA_IBM_TRIEL,
/** Type of bonded interaction is volume conservation force (immersed
Expand Down Expand Up @@ -110,6 +106,8 @@ struct Fene_bond_parameters {
double drmax2;
/** inverse square of @p drmax (internal parameter) */
double drmax2i;

double cutoff() const { return r0 + drmax; }
};

/** Parameters for OIF global forces
Expand All @@ -126,6 +124,8 @@ struct Oif_global_forces_bond_parameters {
double V0;
/** Volume coefficient */
double kv;

double cutoff() const { return -1.; }
};

/** Parameters for OIF local forces
Expand All @@ -151,6 +151,8 @@ struct Oif_local_forces_bond_parameters {
double kal;
/** Viscous coefficient of the triangle vertices */
double kvisc;

double cutoff() const { return -1.; }
};

/** Parameters for harmonic bond Potential */
Expand All @@ -161,6 +163,8 @@ struct Harmonic_bond_parameters {
double r;
/** cutoff bond length */
double r_cut;

double cutoff() const { return r_cut; }
};

/** Parameters for Thermalized bond **/
Expand All @@ -174,9 +178,10 @@ struct Thermalized_bond_parameters {
double pref2_com;
double pref1_dist;
double pref2_dist;

double cutoff() const { return r_cut; }
};

#ifdef ROTATION
/** Parameters for harmonic dumbbell bond Potential */
struct Harmonic_dumbbell_bond_parameters {
/** spring constant */
Expand All @@ -187,26 +192,33 @@ struct Harmonic_dumbbell_bond_parameters {
double r;
/** cutoff bond length */
double r_cut;

double cutoff() const { return r_cut; }
};
#endif

/** Parameters for quartic bond Potential */
struct Quartic_bond_parameters {
double k0, k1;
double r;
double r_cut;

double cutoff() const { return r_cut; }
};

/** Parameters for %Coulomb bond Potential */
struct Bonded_coulomb_bond_parameters {
/** %Coulomb prefactor */
double prefactor;

double cutoff() const { return -1.; }
};

/** Parameters for %Coulomb bond short-range Potential */
struct Bonded_coulomb_sr_bond_parameters {
/** charge factor */
double q1q2;

double cutoff() const { return -1.; }
};

/** Parameters for three-body angular potential (harmonic). */
Expand All @@ -215,6 +227,8 @@ struct Angle_harmonic_bond_parameters {
double bend;
/** equilibrium angle (default is 180 degrees) */
double phi0;

double cutoff() const { return -1.; }
};

/** Parameters for three-body angular potential (cosine). */
Expand All @@ -227,6 +241,8 @@ struct Angle_cosine_bond_parameters {
double cos_phi0;
/** sine of @p phi0 (internal parameter) */
double sin_phi0;

double cutoff() const { return -1.; }
};

/** Parameters for three-body angular potential (cossquare). */
Expand All @@ -237,31 +253,47 @@ struct Angle_cossquare_bond_parameters {
double phi0;
/** cosine of @p phi0 (internal parameter) */
double cos_phi0;

double cutoff() const { return -1.; }
};

/** Parameters for four-body angular potential (dihedral-angle potentials). */
struct Dihedral_bond_parameters {
double mult;
double bend;
double phase;

double cutoff() const { return -1.; }
};

/** Parameters for n-body tabulated potential (n=2,3,4). */
struct Tabulated_bond_parameters {
TabulatedBondedInteraction type;
TabulatedPotential *pot;

double cutoff() const {
switch (type) {
case TAB_BOND_LENGTH:
return assert(pot), pot->cutoff();
default:
return -1.;
};
}
};

#ifdef UMBRELLA
/** Parameters for umbrella potential */
struct Umbrella_bond_parameters {
double k;
int dir;
double r;

double cutoff() const { return std::numeric_limits<double>::infinity(); }
};
#endif

/** Dummy parameters for subtracted-LJ Potential */
struct Subt_lj_bond_parameters {};
struct Subt_lj_bond_parameters {
double cutoff() const { return -1.; }
};

/** Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM */
struct Rigid_bond_parameters {
Expand All @@ -273,6 +305,8 @@ struct Rigid_bond_parameters {
/**Velocity Tolerance/Accuracy for termination of RATTLE/SHAKE iterations
* during velocity corrections */
double v_tol;

double cutoff() const { return std::sqrt(d2); }
};

enum class tElasticLaw { NeoHookean, Skalak };
Expand All @@ -299,6 +333,8 @@ struct IBM_Triel_Parameters {
tElasticLaw elasticLaw;
double k1;
double k2;

double cutoff() const { return maxDist; }
};

/** Parameters for IBM volume conservation bond **/
Expand All @@ -309,10 +345,8 @@ struct IBM_VolCons_Parameters {
double volRef;
/** Spring constant for volume force */
double kappaV;
// Whether to write out center-of-mass at each time step
// Actually this is more of an analysis function and does not strictly belong
// to volume conservation
// bool writeCOM;

double cutoff() const { return -1.; }
};

/** Parameters for IBM tribend **/
Expand All @@ -322,6 +356,12 @@ struct IBM_Tribend_Parameters {

/** Reference angle */
double theta0;

double cutoff() const { return -1.; }
};

struct VirtualBond_Parameters {
double cutoff() const { return -1.; }
};

/** Union in which to store the parameters of an individual bonded interaction
Expand All @@ -331,9 +371,7 @@ union Bond_parameters {
Oif_global_forces_bond_parameters oif_global_forces;
Oif_local_forces_bond_parameters oif_local_forces;
Harmonic_bond_parameters harmonic;
#ifdef ROTATION
Harmonic_dumbbell_bond_parameters harmonic_dumbbell;
#endif
Quartic_bond_parameters quartic;
Bonded_coulomb_bond_parameters bonded_coulomb;
Bonded_coulomb_sr_bond_parameters bonded_coulomb_sr;
Expand All @@ -342,15 +380,14 @@ union Bond_parameters {
Angle_cossquare_bond_parameters angle_cossquare;
Dihedral_bond_parameters dihedral;
Tabulated_bond_parameters tab;
#ifdef UMBRELLA
Umbrella_bond_parameters umbrella;
#endif
Thermalized_bond_parameters thermalized_bond;
Subt_lj_bond_parameters subt_lj;
Rigid_bond_parameters rigid_bond;
IBM_Triel_Parameters ibm_triel;
IBM_VolCons_Parameters ibmVolConsParameters;
IBM_Tribend_Parameters ibm_tribend;
VirtualBond_Parameters virt;
};

/** Defines parameters for a bonded interaction. */
Expand Down
1 change: 1 addition & 0 deletions src/core/bonded_interactions/bonded_tab.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ int tabulated_bonded_set_params(int bond_type,

/* set types */
auto tab_pot = bonded_ia_params[bond_type].p.tab.pot = new TabulatedPotential;
bonded_ia_params[bond_type].p.tab.type = tab_type;

/* set number of interaction partners */
switch (tab_type) {
Expand Down

0 comments on commit 780c7fd

Please sign in to comment.