From fa1e6a3c926a8bf98ca0e422e4f6b5f48db5bdd5 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 30 Jan 2020 15:34:19 +0100 Subject: [PATCH 1/2] core: Removed unused bonded ia type flag --- src/core/bonded_interactions/bonded_interaction_data.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/core/bonded_interactions/bonded_interaction_data.hpp b/src/core/bonded_interactions/bonded_interaction_data.hpp index 1f8ed52b7d0..aea34580745 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.hpp +++ b/src/core/bonded_interactions/bonded_interaction_data.hpp @@ -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 From a757875f82f5be6816515032e2d1c093fe81338b Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 30 Jan 2020 18:02:20 +0100 Subject: [PATCH 2/2] core: Define cutoff for all bond types --- .../bonded_interaction_data.cpp | 116 +++++++++++------- .../bonded_interaction_data.hpp | 67 ++++++++-- src/core/bonded_interactions/bonded_tab.cpp | 1 + 3 files changed, 127 insertions(+), 57 deletions(-) diff --git a/src/core/bonded_interactions/bonded_interaction_data.cpp b/src/core/bonded_interactions/bonded_interaction_data.cpp index b7d4fb6a4d2..efb1b8439b7 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.cpp +++ b/src/core/bonded_interactions/bonded_interaction_data.cpp @@ -19,59 +19,86 @@ #include "bonded_interaction_data.hpp" #include "communication.hpp" +#include #include std::vector 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) { @@ -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); diff --git a/src/core/bonded_interactions/bonded_interaction_data.hpp b/src/core/bonded_interactions/bonded_interaction_data.hpp index aea34580745..dc4508e8565 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.hpp +++ b/src/core/bonded_interactions/bonded_interaction_data.hpp @@ -106,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 @@ -122,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 @@ -147,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 */ @@ -157,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 **/ @@ -170,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 */ @@ -183,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). */ @@ -211,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). */ @@ -223,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). */ @@ -233,6 +253,8 @@ 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). */ @@ -240,24 +262,38 @@ 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::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 { @@ -269,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 }; @@ -295,6 +333,8 @@ struct IBM_Triel_Parameters { tElasticLaw elasticLaw; double k1; double k2; + + double cutoff() const { return maxDist; } }; /** Parameters for IBM volume conservation bond **/ @@ -305,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 **/ @@ -318,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 @@ -327,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; @@ -338,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. */ diff --git a/src/core/bonded_interactions/bonded_tab.cpp b/src/core/bonded_interactions/bonded_tab.cpp index b2df85025e6..c0ffee2a69e 100644 --- a/src/core/bonded_interactions/bonded_tab.cpp +++ b/src/core/bonded_interactions/bonded_tab.cpp @@ -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) {