From 484bc3d0afe37853308717d020db570e1a90a6b1 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Fri, 7 Aug 2020 15:52:21 -0400 Subject: [PATCH 01/15] Add support to cache which atoms within/outside volumetric maps The refactoring of the NAMD proxy backend leads to failing the 007_map_total_internal regtest. This is due once again to the combination of compiler optimizations and GridForceGrid being based on float instead of double. Updating the test's reference files in a separate commit. --- namd/src/colvarproxy_namd.C | 78 ++++++++++++++++++++++++++++++------- namd/src/colvarproxy_namd.h | 9 +++-- src/colvarcomp.h | 3 ++ src/colvarcomp_volmaps.cpp | 17 +++++--- src/colvarproxy_volmaps.cpp | 3 +- src/colvarproxy_volmaps.h | 8 +++- vmd/src/colvarproxy_vmd.C | 6 ++- vmd/src/colvarproxy_vmd.h | 3 +- 8 files changed, 100 insertions(+), 27 deletions(-) diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index dbcadbceb..5adc9d9d2 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -1435,18 +1435,38 @@ void colvarproxy_namd::GridForceGridLoop(T const *g, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) + cvm::real *atom_field, + int *inside) { float V = 0.0f; Vector dV(0.0); - int i = 0; + int i = 0, status = -1; cvm::atom_iter ai = atom_begin; for ( ; ai != atom_end; ai++, i++) { + + if ((flags & volmap_flag_use_atomlist) && + !(flags & volmap_flag_rebuild_atomlist)) { + if (inside[i] == 0) { + // Skip atom according to precomputed list + continue; + } + } + + // TODO look into compute_V() to skip gradient computation + if (g->compute_VdV(Position(ai->pos.x, ai->pos.y, ai->pos.z), V, dV)) { // out-of-bounds atom + if (flags & volmap_flag_rebuild_atomlist) { + inside[i] = 0; + } V = 0.0f; dV = 0.0; } else { + + if (flags & volmap_flag_rebuild_atomlist) { + inside[i] = 1; + } + if (flags & volmap_flag_use_atom_field) { *value += V * atom_field[i]; if (flags & volmap_flag_gradients) { @@ -1469,16 +1489,47 @@ void colvarproxy_namd::getGridForceGridValue(int flags, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) + cvm::real *atom_field, + int *inside) { - if (flags & volmap_flag_use_atom_field) { - int const new_flags = volmap_flag_use_atom_field | volmap_flag_gradients; - GridForceGridLoop(g, atom_begin, atom_end, - value, atom_field); + if (atom_field) { + + if (flags & volmap_flag_use_atomlist) { + if (flags & volmap_flag_rebuild_atomlist) { + int const new_flags = volmap_flag_use_atom_field | + volmap_flag_use_atomlist | volmap_flag_rebuild_atomlist; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } else { + int const new_flags = volmap_flag_use_atom_field | + volmap_flag_use_atomlist; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } + } else { + int const new_flags = volmap_flag_use_atom_field; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } + } else { - int const new_flags = volmap_flag_gradients; - GridForceGridLoop(g, atom_begin, atom_end, - value, atom_field); + + if (flags & volmap_flag_use_atomlist) { + if (flags & volmap_flag_rebuild_atomlist) { + int const new_flags = volmap_flag_use_atomlist | + volmap_flag_rebuild_atomlist; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } else { + int const new_flags = volmap_flag_use_atomlist; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } + } else { + int const new_flags = volmap_flag_null; + GridForceGridLoop(g, atom_begin, atom_end, + value, atom_field, inside); + } } } @@ -1488,7 +1539,8 @@ int colvarproxy_namd::compute_volmap(int flags, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) + cvm::real *atom_field, + int *inside) { Molecule *mol = Node::Object()->molecule; GridforceGrid *grid = mol->get_gridfrc_grid(volmap_id); @@ -1496,11 +1548,11 @@ int colvarproxy_namd::compute_volmap(int flags, if (grid->get_grid_type() == GridforceGrid::GridforceGridTypeFull) { GridforceFullMainGrid *g = dynamic_cast(grid); getGridForceGridValue(flags, g, atom_begin, atom_end, - value, atom_field); + value, atom_field, inside); } else if (grid->get_grid_type() == GridforceGrid::GridforceGridTypeLite) { GridforceLiteGrid *g = dynamic_cast(grid); getGridForceGridValue(flags, g, atom_begin, atom_end, - value, atom_field); + value, atom_field, inside); } return COLVARS_OK; } diff --git a/namd/src/colvarproxy_namd.h b/namd/src/colvarproxy_namd.h index 9b4eba093..2bf74cbc7 100644 --- a/namd/src/colvarproxy_namd.h +++ b/namd/src/colvarproxy_namd.h @@ -242,7 +242,8 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster { cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) override; + cvm::real *atom_field, + int *inside) override; /// Abstraction of the two types of NAMD volumetric maps template @@ -251,7 +252,8 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster { cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field); + cvm::real *atom_field, + int *inside); /// Implementation of inner loop; allows for atom list computation and use template @@ -259,7 +261,8 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster { cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field); + cvm::real *atom_field, + int *inside); #endif diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 53755576c..8c6639d6e 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -1571,6 +1571,9 @@ class colvar::map_total /// Weights assigned to each atom (default: uniform weights) std::vector atom_weights; + + /// Flags recording whether the atoms are within the map or not + std::vector atoms_inside; }; diff --git a/src/colvarcomp_volmaps.cpp b/src/colvarcomp_volmaps.cpp index a15501ab6..640b070b4 100644 --- a/src/colvarcomp_volmaps.cpp +++ b/src/colvarcomp_volmaps.cpp @@ -21,6 +21,7 @@ colvar::map_total::map_total() } + int colvar::map_total::init(std::string const &conf) { int error_code = cvc::init(conf); @@ -90,17 +91,23 @@ void colvar::map_total::calc_value() int flags = is_enabled(f_cvc_gradient) ? colvarproxy::volmap_flag_gradients : colvarproxy::volmap_flag_null; - if (atoms != NULL) { + if (atoms) { // Compute the map inside Colvars x.real_value = 0.0; - cvm::real *w = NULL; + if ((flags & colvarproxy::volmap_flag_use_atomlist) && (atoms_inside.size() != atoms->size())) { + atoms_inside.resize(atoms->size()); + } + + cvm::real *w = nullptr; if (atom_weights.size() > 0) { flags |= colvarproxy::volmap_flag_use_atom_field; - w = &(atom_weights[0]); + w = atom_weights.data(); } - proxy->compute_volmap(flags, volmap_id, atoms->begin(), atoms->end(), - &(x.real_value), w); + + proxy->compute_volmap(flags, volmap_id, atoms->begin(), atoms->end(), &(x.real_value), w, + atoms_inside.data()); + } else { // Get the externally computed value x.real_value = proxy->get_volmap_value(volmap_index); diff --git a/src/colvarproxy_volmaps.cpp b/src/colvarproxy_volmaps.cpp index 3d02d8085..d9865d7d0 100644 --- a/src/colvarproxy_volmaps.cpp +++ b/src/colvarproxy_volmaps.cpp @@ -114,7 +114,8 @@ int colvarproxy_volmaps::compute_volmap(int /* flags */, cvm::atom_iter /* atom_begin */, cvm::atom_iter /* atom_end */, cvm::real * /* value */, - cvm::real * /* atom_field */) + cvm::real * /* atom_field */, + int * /* inside */) { return COLVARS_NOT_IMPLEMENTED; } diff --git a/src/colvarproxy_volmaps.h b/src/colvarproxy_volmaps.h index 6e1ffcc7c..88a9c9d54 100644 --- a/src/colvarproxy_volmaps.h +++ b/src/colvarproxy_volmaps.h @@ -81,18 +81,22 @@ class colvarproxy_volmaps { /// \param atom_end Iterator pointing past the last atom /// \param value Pointer to location of total to increment /// \param atom_field Array of atomic field values (if NULL, ones are used) + /// \param inside Array of flags (1 = inside map, 0 = outside) virtual int compute_volmap(int flags, int volmap_id, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field); + cvm::real *atom_field, + int *inside); /// Flags controlling what computation is done on the map enum { volmap_flag_null = 0, volmap_flag_gradients = 1, - volmap_flag_use_atom_field = (1<<8) + volmap_flag_use_atom_field = (1<<8), + volmap_flag_use_atomlist = (1<<9), + volmap_flag_rebuild_atomlist = (1<<10) }; /// Compute the root-mean-square of the applied forces diff --git a/vmd/src/colvarproxy_vmd.C b/vmd/src/colvarproxy_vmd.C index 9ec1d948c..d6b425bdb 100644 --- a/vmd/src/colvarproxy_vmd.C +++ b/vmd/src/colvarproxy_vmd.C @@ -808,7 +808,8 @@ void colvarproxy_vmd::compute_voldata(VolumetricData const *voldata, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) + cvm::real *atom_field, + int * /* inside */) { int i = 0; float coord[3], voxcoord[3], grad[3]; @@ -872,7 +873,8 @@ int colvarproxy_vmd::compute_volmap(int flags, cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field) + cvm::real *atom_field, + int * /* inside */) { int error_code = COLVARS_OK; VolumetricData const *voldata = vmdmol->get_volume_data(volmap_id); diff --git a/vmd/src/colvarproxy_vmd.h b/vmd/src/colvarproxy_vmd.h index e56eacd89..469e81789 100644 --- a/vmd/src/colvarproxy_vmd.h +++ b/vmd/src/colvarproxy_vmd.h @@ -126,7 +126,8 @@ class colvarproxy_vmd : public colvarproxy { cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field); + cvm::real *atom_field, + int *inside); /// Get value of alchemical lambda parameter from back-end (stub) int get_alch_lambda(cvm::real* lambda) { From d7c18fb36c9368d02c4a06b796e80a881ae4faad Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 10 Aug 2020 19:06:48 -0400 Subject: [PATCH 02/15] Add atomListFrequency keyword to mapTotal --- src/colvarcomp.cpp | 5 +++++ src/colvarcomp.h | 3 +++ src/colvarcomp_volmaps.cpp | 15 ++++++++++++++- src/colvardeps.h | 2 ++ 4 files changed, 24 insertions(+), 1 deletion(-) diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 9a056f7dd..673265903 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -274,6 +274,8 @@ int colvar::cvc::init_dependencies() { init_feature(f_cvc_one_site_total_force, "total_force_from_one_group", f_type_user); require_feature_self(f_cvc_one_site_total_force, f_cvc_com_based); + init_feature(f_cvc_dynamic_atom_list, "dynamic_atom_list", f_type_dynamic); + init_feature(f_cvc_com_based, "function_of_centers_of_mass", f_type_static); init_feature(f_cvc_pbc_minimum_image, "use_minimum-image_with_PBCs", f_type_user); @@ -333,6 +335,9 @@ int colvar::cvc::init_dependencies() { // Features that are implemented by default if their requirements are feature_states[f_cvc_one_site_total_force].available = true; + // By default the list of atoms contributing to the CVC is fixed + feature_states[f_cvc_dynamic_atom_list].available = false; + // Features That are implemented only for certain simulation engine configurations feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK); feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available; diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 8c6639d6e..d0a74da0f 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -300,6 +300,9 @@ class colvar::cvc /// \brief CVC-specific default colvar width (default: not provided) cvm::real width = 0.0; + + /// Frequency at which the list of contributing atoms will be updated + size_t atom_list_freq = 0; }; diff --git a/src/colvarcomp_volmaps.cpp b/src/colvarcomp_volmaps.cpp index 640b070b4..c34b683c2 100644 --- a/src/colvarcomp_volmaps.cpp +++ b/src/colvarcomp_volmaps.cpp @@ -17,11 +17,11 @@ colvar::map_total::map_total() { set_function_type("mapTotal"); + provide(f_cvc_dynamic_atom_list); x.type(colvarvalue::type_scalar); } - int colvar::map_total::init(std::string const &conf) { int error_code = cvc::init(conf); @@ -49,6 +49,11 @@ int colvar::map_total::init(std::string const &conf) error_code |= proxy->check_volmap_by_id(volmap_id); } + get_keyval(conf, "atomListFrequency", atom_list_freq, atom_list_freq); + if (atom_list_freq > 0) { + enable(f_cvc_dynamic_atom_list); + } + } else { // Using selection from the MD engine @@ -87,10 +92,18 @@ int colvar::map_total::init(std::string const &conf) void colvar::map_total::calc_value() { + colvarproxy *proxy = cvm::main()->proxy; int flags = is_enabled(f_cvc_gradient) ? colvarproxy::volmap_flag_gradients : colvarproxy::volmap_flag_null; + if (atom_list_freq > 0) { + flags |= colvarproxy::volmap_flag_use_atomlist; + if (cvm::step_relative() % atom_list_freq == 0) { + flags |= colvarproxy::volmap_flag_rebuild_atomlist; + } + } + if (atoms) { // Compute the map inside Colvars x.real_value = 0.0; diff --git a/src/colvardeps.h b/src/colvardeps.h index 92e7a8832..b51be88b2 100644 --- a/src/colvardeps.h +++ b/src/colvardeps.h @@ -386,6 +386,8 @@ class colvardeps { /// With PBCs, minimum-image convention will be used for distances /// (does not affect the periodicity of CVC values, e.g. angles) f_cvc_pbc_minimum_image, + /// The list of atoms contributing being processed may be updated + f_cvc_dynamic_atom_list, /// This CVC is a function of centers of mass f_cvc_com_based, /// This CVC can be computed in parallel From 307c188cc4ceadd7c9beacbc4eb77ee2ba1cd0d8 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 17:15:14 -0400 Subject: [PATCH 03/15] Allow clear_atom() to return error code --- namd/src/colvarproxy_namd.C | 4 ++-- namd/src/colvarproxy_namd.h | 2 +- src/colvarproxy.cpp | 13 +++++++++---- src/colvarproxy.h | 5 ++--- 4 files changed, 14 insertions(+), 10 deletions(-) diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index 5adc9d9d2..25424d90e 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -836,9 +836,9 @@ int colvarproxy_namd::init_atom(cvm::residue_id const &residue, } -void colvarproxy_namd::clear_atom(int index) +int colvarproxy_namd::clear_atom(int index) { - colvarproxy::clear_atom(index); + return colvarproxy::clear_atom(index); // TODO remove it from GlobalMaster arrays? } diff --git a/namd/src/colvarproxy_namd.h b/namd/src/colvarproxy_namd.h index 2bf74cbc7..5bd6b59fe 100644 --- a/namd/src/colvarproxy_namd.h +++ b/namd/src/colvarproxy_namd.h @@ -193,7 +193,7 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster { int check_atom_id(cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id) override; - void clear_atom(int index) override; + int clear_atom(int index) override; void update_atom_properties(int index); diff --git a/src/colvarproxy.cpp b/src/colvarproxy.cpp index 1ed7a5555..25d1375cc 100644 --- a/src/colvarproxy.cpp +++ b/src/colvarproxy.cpp @@ -99,15 +99,20 @@ int colvarproxy_atoms::check_atom_id(cvm::residue_id const &residue, } -void colvarproxy_atoms::clear_atom(int index) +int colvarproxy_atoms::clear_atom(int index) { - if (((size_t) index) >= atoms_ids.size()) { - cvm::error("Error: trying to disable an atom that was not previously requested.\n", - COLVARS_INPUT_ERROR); + if (static_cast(index) >= atoms_ids.size()) { + return cvm::error("Error: trying to unrequest an atom that was not " + "previously requested.\n", COLVARS_BUG_ERROR); + } + if (index < 0) { + return cvm::error("Error: invalid argument to clear_atom().\n", + COLVARS_BUG_ERROR); } if (atoms_refcount[index] > 0) { atoms_refcount[index] -= 1; } + return COLVARS_OK; } diff --git a/src/colvarproxy.h b/src/colvarproxy.h index 353f354ef..90ae30606 100644 --- a/src/colvarproxy.h +++ b/src/colvarproxy.h @@ -73,9 +73,8 @@ class colvarproxy_atoms { std::string const &atom_name, std::string const &segment_id); - /// \brief Used by the atom class destructor: rather than deleting the array slot - /// (costly) set the corresponding atoms_refcount to zero - virtual void clear_atom(int index); + /// Used by the atom class destructor: set atoms_refcount entry to zero + virtual int clear_atom(int index); /// Clear atomic data int reset(); From 8cbbdee3cdba82395e74a0cc0d547ee6e411d69a Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 17:16:21 -0400 Subject: [PATCH 04/15] Don't send forces to NAMD for atoms with zero refcount --- namd/src/colvarproxy_namd.C | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index 25424d90e..5a9591bc3 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -566,11 +566,14 @@ void colvarproxy_namd::calculate() print_output_atomic_data(); } - // communicate all forces to the MD integrator + // communicate all forces to NAMD for (size_t i = 0; i < atoms_ids.size(); i++) { - cvm::rvector const &f = atoms_new_colvar_forces[i]; - modifyForcedAtoms().add(atoms_ids[i]); - modifyAppliedForces().add(Vector(f.x, f.y, f.z)); + if (atoms_refcount[i] > 0) { + // Add a force only if the atom is currently used + cvm::rvector const &f = atoms_new_colvar_forces[i]; + modifyForcedAtoms().add(atoms_ids[i]); + modifyAppliedForces().add(Vector(f.x, f.y, f.z)); + } } if (atom_groups_new_colvar_forces.size() > 0) { From 62e1c97abec2251dc6e3d3af142098205dea8f00 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 17:21:03 -0400 Subject: [PATCH 05/15] Unrequest atom from NAMD when not needed any more --- namd/src/colvarproxy_namd.C | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index 5a9591bc3..19d952cd3 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -841,8 +841,18 @@ int colvarproxy_namd::init_atom(cvm::residue_id const &residue, int colvarproxy_namd::clear_atom(int index) { - return colvarproxy::clear_atom(index); - // TODO remove it from GlobalMaster arrays? + int error_code = colvarproxy::clear_atom(index); + if (error_code == COLVARS_OK) { + if (atoms_refcount[index] == 0) { + // Clear this atom entry from the requested atoms + int aid = atoms_ids[index]; + int const aid_index_in_gm = modifyRequestedAtoms().find(aid); + if (aid_index_in_gm >= 0) { + modifyRequestedAtoms().del(aid_index_in_gm, 1); + } + } + } + return error_code; } From 0a19c105537e5e0bbe0b13a89ef7da79263fbb52 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 18:24:19 -0400 Subject: [PATCH 06/15] Add support for reupdating the atom list in the proxy --- namd/src/colvarproxy_namd.C | 42 ++++++++++++++++++++++++++++++++++++- namd/src/colvarproxy_namd.h | 8 +++++-- src/colvaratoms.h | 16 +++++++++----- src/colvarproxy.cpp | 30 ++++++++++++++++++++++++++ src/colvarproxy.h | 23 +++++++++++++++++--- vmd/src/colvarproxy_vmd.h | 4 +++- 6 files changed, 111 insertions(+), 12 deletions(-) diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index 19d952cd3..a18fdc166 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -202,6 +202,10 @@ int colvarproxy_namd::update_atoms_map(AtomIDList::const_iterator begin, init_atoms_map(); } + if (cvm::debug()) { + cvm::log("Updating atoms_map for "+cvm::to_str(begin - end)+" atoms.\n"); + } + for (AtomIDList::const_iterator a_i = begin; a_i != end; a_i++) { if (atoms_map[*a_i] >= 0) continue; @@ -215,7 +219,7 @@ int colvarproxy_namd::update_atoms_map(AtomIDList::const_iterator begin, if (atoms_map[*a_i] < 0) { // this atom is probably managed by another GlobalMaster: - // add it here anyway to avoid having to test for array boundaries at each step + // add it here anyway so that Colvars can ensure it is requested int const index = add_atom_slot(*a_i); atoms_map[*a_i] = index; modifyRequestedAtoms().add(*a_i); @@ -615,6 +619,17 @@ void colvarproxy_namd::calculate() #endif #endif + if (atom_list_frequency() > 0) { + if (((cvm::step_relative()+1) % atom_list_frequency()) == 0) { + // Before all-atom evaluation + update_requested_atoms(); + } + if ((cvm::step_relative() % atom_list_frequency()) == 0) { + // After all-atom computation + update_requested_atoms(); + } + } + // NAMD does not destruct GlobalMaster objects, so we must remember // to write all output files at the end of a run if (step == simparams->N) { @@ -622,6 +637,30 @@ void colvarproxy_namd::calculate() } } + +int colvarproxy_namd::update_requested_atoms() +{ + int error_code = COLVARS_OK; + if (cvm::debug()) { + cvm::log("Updating list of requested atoms from NAMD.\n"); + cvm::log("Before: "+cvm::to_str(modifyRequestedAtoms().size())+ + " elements.\n"); + } + modifyRequestedAtoms().clear(); + for (size_t i = 0; i < atoms_ids.size(); i++) { + if (atoms_refcount[i] > 0) { + modifyRequestedAtoms().add(atoms_ids[i]); + } + } + if (cvm::debug()) { + cvm::log("After: "+cvm::to_str(modifyRequestedAtoms().size())+ + " elements.\n"); + } + + return COLVARS_OK; +} + + void colvarproxy_namd::update_accelMD_info() { // This aMD factor is from previous step! amd_weight_factor = std::exp(controller->accelMDdV / (target_temperature() * boltzmann())); @@ -1471,6 +1510,7 @@ void colvarproxy_namd::GridForceGridLoop(T const *g, // out-of-bounds atom if (flags & volmap_flag_rebuild_atomlist) { inside[i] = 0; + decrease_refcount(ai->array_index()); } V = 0.0f; dV = 0.0; diff --git a/namd/src/colvarproxy_namd.h b/namd/src/colvarproxy_namd.h index 5bd6b59fe..8cfbf5220 100644 --- a/namd/src/colvarproxy_namd.h +++ b/namd/src/colvarproxy_namd.h @@ -87,11 +87,15 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster { /// Allocate an atoms map with the same size as the NAMD topology void init_atoms_map(); - // synchronize the local arrays with requested or forced atoms + /// Rebuild the list of requested atoms based on Colvars-internal refcounts + int update_requested_atoms(); + + /// Synchronize the local arrays with requested or forced atoms int update_atoms_map(AtomIDList::const_iterator begin, AtomIDList::const_iterator end); - void calculate(); + /// Overrides GlobalMaster::calculate() + void calculate() override; void log(std::string const &message) override; void error(std::string const &message) override; diff --git a/src/colvaratoms.h b/src/colvaratoms.h index 528e849df..5c26233d4 100644 --- a/src/colvaratoms.h +++ b/src/colvaratoms.h @@ -30,11 +30,6 @@ struct rotation_derivative; class colvarmodule::atom { -protected: - - /// Index in the colvarproxy arrays (\b NOT in the global topology!) - int index; - public: /// Identifier for the MD program (0-based) @@ -103,6 +98,12 @@ class colvarmodule::atom { vel = grad = total_force = cvm::rvector(0.0); } + /// Index of this atom in the internal arrays + inline int array_index() const + { + return index; + } + /// Get the latest value of the mass inline void update_mass() { @@ -148,6 +149,11 @@ class colvarmodule::atom { { (cvm::proxy)->apply_atom_force(index, new_force); } + +protected: + + /// Index in the colvarproxy arrays (\b NOT in the global topology!) + int index; }; diff --git a/src/colvarproxy.cpp b/src/colvarproxy.cpp index 25d1375cc..47a872cac 100644 --- a/src/colvarproxy.cpp +++ b/src/colvarproxy.cpp @@ -26,6 +26,7 @@ colvarproxy_atoms::colvarproxy_atoms() atoms_max_applied_force_id_ = -1; modified_atom_list_ = false; updated_masses_ = updated_charges_ = false; + atom_list_freq_ = 0; } @@ -126,6 +127,35 @@ size_t colvarproxy_atoms::get_num_active_atoms() const } +int colvarproxy_atoms::set_atom_list_frequency(int atom_list_freq) +{ + int const proxy_freq = atom_list_frequency(); + if (proxy_freq > 0) { + if (atom_list_freq < proxy_freq) { + // Use the shortest number of steps among all variables + if ((atom_list_freq % proxy_freq) != 0) { + return cvm::error("Error: all values of atomListFrequency should be " + "multiples of the same number (proxy_freq = "+ + cvm::to_str(proxy_freq)+").\n", + COLVARS_INPUT_ERROR); + } else { + atom_list_frequency() = atom_list_freq; + } + } else { + if ((proxy_freq % atom_list_freq) != 0) { + return cvm::error("Error: all values of atomListFrequency should be " + "multiples of the same number (proxy_freq = "+ + cvm::to_str(proxy_freq)+").\n", + COLVARS_INPUT_ERROR); + } + } + } else { + atom_list_frequency() = atom_list_freq; + } + return COLVARS_OK; +} + + void colvarproxy_atoms::compute_rms_atoms_applied_force() { atoms_rms_applied_force_ = diff --git a/src/colvarproxy.h b/src/colvarproxy.h index 90ae30606..b17230cd5 100644 --- a/src/colvarproxy.h +++ b/src/colvarproxy.h @@ -100,6 +100,15 @@ class colvarproxy_atoms { atoms_refcount[index] += 1; } + /// Decrease the reference count of the given atom + /// \param index Internal index in the Colvars arrays + inline void decrease_refcount(int index) + { + if (atoms_refcount[index] > 0) { + atoms_refcount[index] -= 1; + } + } + /// Get the charge of the given atom /// \param index Internal index in the Colvars arrays inline cvm::real get_atom_charge(int index) const @@ -248,12 +257,17 @@ class colvarproxy_atoms { return updated_charges_; } + /// Request/unrequest atoms at this frequency + inline int & atom_list_frequency() + { + return atom_list_freq_; + } + protected: - /// \brief Array of 0-based integers used to uniquely associate atoms - /// within the host program + /// Array of integers used to identify atoms in the engine std::vector atoms_ids; - /// \brief Keep track of how many times each atom is used by a separate colvar object + /// Keep track of how many times each atom is used by a separate colvar object std::vector atoms_refcount; /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) std::vector atoms_masses; @@ -281,6 +295,9 @@ class colvarproxy_atoms { /// Whether the masses and charges have been updated from the host code bool updated_masses_, updated_charges_; + /// Request/unrequest atoms at this frequency + int atom_list_freq_; + /// Used by all init_atom() functions: create a slot for an atom not /// requested yet; returns the index in the arrays int add_atom_slot(int atom_id); diff --git a/vmd/src/colvarproxy_vmd.h b/vmd/src/colvarproxy_vmd.h index 469e81789..3bde72255 100644 --- a/vmd/src/colvarproxy_vmd.h +++ b/vmd/src/colvarproxy_vmd.h @@ -119,8 +119,10 @@ class colvarproxy_vmd : public colvarproxy { cvm::atom_iter atom_begin, cvm::atom_iter atom_end, cvm::real *value, - cvm::real *atom_field); + cvm::real *atom_field, + int *inside); + /// Implements loops (does not support yet atom lists) template void compute_voldata(VolumetricData const *voldata, cvm::atom_iter atom_begin, From 45da9ede296356909289db1c3e0b3304a7cba987 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 17:24:57 -0400 Subject: [PATCH 07/15] Define atom_list_freq for all CVCs, let them coordinate with the proxy --- src/colvar.cpp | 11 +++++++ src/colvar.h | 3 ++ src/colvarcomp.cpp | 64 ++++++++++++++++++++++++++++++++++++-- src/colvarcomp.h | 11 ++++++- src/colvarcomp_volmaps.cpp | 6 +--- src/colvarmodule.cpp | 5 +++ src/colvarproxy.h | 3 ++ 7 files changed, 95 insertions(+), 8 deletions(-) diff --git a/src/colvar.cpp b/src/colvar.cpp index 136aa2afd..bf83162bc 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -1739,6 +1739,17 @@ int colvar::collect_cvc_Jacobians() } +int colvar::update_requested_atoms() +{ + int error_code = COLVARS_OK; + size_t i, cvc_count; + for (size_t i = 0; i < cvcs.size(); i++) { + error_code |= (cvcs[i])->update_all_requested_atoms(); + } + return error_code; +} + + int colvar::calc_colvar_properties() { if (is_enabled(f_cv_fdiff_velocity)) { diff --git a/src/colvar.h b/src/colvar.h index 1db313f41..228d7967e 100644 --- a/src/colvar.h +++ b/src/colvar.h @@ -333,6 +333,9 @@ class colvar : public colvarparse, public colvardeps { /// \brief Calculate the quantities associated to the colvar (but not to the CVCs) int calc_colvar_properties(); + /// If the variable doesn't need all atoms at every step, update their list + int update_requested_atoms(); + /// Get the current applied force inline colvarvalue const applied_force() const { diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 673265903..f8fdd3f18 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -124,8 +124,19 @@ int colvar::cvc::init(std::string const &conf) enable(f_cvc_pbc_minimum_image); } - // Attempt scalable calculations when in parallel? (By default yes, if available) - get_keyval(conf, "scalable", b_try_scalable, b_try_scalable); + if (is_available(f_cvc_scalable)) { + // Attempt scalable calculations when in parallel? (By default yes) + get_keyval(conf, "scalable", b_try_scalable, b_try_scalable); + } + + if (is_available(f_cvc_dynamic_atom_list)) { + get_keyval(conf, "atomListFrequency", atom_list_freq, atom_list_freq); + if (atom_list_freq > 0) { + enable(f_cvc_dynamic_atom_list); + // Set the parameter and enable its dependencies + error_code |= set_atom_list_frequency(atom_list_freq); + } + } if (cvm::debug()) cvm::log("Done initializing cvc base object.\n"); @@ -134,6 +145,55 @@ int colvar::cvc::init(std::string const &conf) } +int colvar::cvc::set_atom_list_frequency(int new_frequency) +{ + atom_list_freq = new_frequency; + return cvm::main()->proxy->set_atom_list_frequency(atom_list_freq); +} + + +int colvar::cvc::update_requested_atoms(cvm::atom_group *dyn_atoms) +{ + int error_code = COLVARS_OK; + colvarproxy *proxy = cvm::main()->proxy; + + if (atom_list_freq > 0) { + + // Reenable all atoms for the next step + if (((cvm::step_absolute()+1) % atom_list_freq) == 0) { + for (cvm::atom_iter ai = dyn_atoms->begin(); + ai != dyn_atoms->end(); ai++) { + proxy->increase_refcount(ai->array_index()); + } + } + + if (!is_enabled(f_cvc_dynamic_atom_list)) { + // If the CVC is not enabling/disabling atoms on its own, then disable + // them all for the next step + if (((cvm::step_absolute()) % atom_list_freq) == 0) { + for (cvm::atom_iter ai = dyn_atoms->begin(); + ai != dyn_atoms->end(); ai++) { + proxy->decrease_refcount(ai->array_index()); + } + } + } + } + + return error_code; +} + + +int colvar::cvc::update_all_requested_atoms() +{ + int error_code = COLVARS_OK; + for (std::vector::iterator agi = atom_groups.begin(); + agi != atom_groups.end(); agi++) { + error_code |= update_requested_atoms(*agi); + } + return COLVARS_OK; +} + + int colvar::cvc::init_total_force_params(std::string const &conf) { if (cvm::get_error()) return COLVARS_ERROR; diff --git a/src/colvarcomp.h b/src/colvarcomp.h index d0a74da0f..d45dca06d 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -243,6 +243,15 @@ class colvar::cvc } } + /// Set frequency at which the list of contributing atoms will be updated + virtual int set_atom_list_frequency(int new_frequency); + + /// If needed, update the requested list of atoms from this group + virtual int update_requested_atoms(cvm::atom_group *dyn_atoms); + + /// If needed, update the requested list of atoms from all groups + virtual int update_all_requested_atoms(); + protected: /// Set the value of \link function_type \endlink and its dependencies @@ -302,7 +311,7 @@ class colvar::cvc cvm::real width = 0.0; /// Frequency at which the list of contributing atoms will be updated - size_t atom_list_freq = 0; + int atom_list_freq = 0; }; diff --git a/src/colvarcomp_volmaps.cpp b/src/colvarcomp_volmaps.cpp index c34b683c2..3babf9258 100644 --- a/src/colvarcomp_volmaps.cpp +++ b/src/colvarcomp_volmaps.cpp @@ -49,11 +49,6 @@ int colvar::map_total::init(std::string const &conf) error_code |= proxy->check_volmap_by_id(volmap_id); } - get_keyval(conf, "atomListFrequency", atom_list_freq, atom_list_freq); - if (atom_list_freq > 0) { - enable(f_cvc_dynamic_atom_list); - } - } else { // Using selection from the MD engine @@ -118,6 +113,7 @@ void colvar::map_total::calc_value() w = atom_weights.data(); } + // Note: this function also updates the list of requested atoms proxy->compute_volmap(flags, volmap_id, atoms->begin(), atoms->end(), &(x.real_value), w, atoms_inside.data()); diff --git a/src/colvarmodule.cpp b/src/colvarmodule.cpp index 34485d788..79f53826d 100644 --- a/src/colvarmodule.cpp +++ b/src/colvarmodule.cpp @@ -1008,6 +1008,11 @@ int colvarmodule::calc_colvars() cvm::decrease_depth(); } + // Update requested atoms for variables that support it (no-op otherwise) + for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) { + error_code |= (*cvi)->update_requested_atoms(); + } + error_code |= cvm::get_error(); return error_code; } diff --git a/src/colvarproxy.h b/src/colvarproxy.h index b17230cd5..c4a09e094 100644 --- a/src/colvarproxy.h +++ b/src/colvarproxy.h @@ -76,6 +76,9 @@ class colvarproxy_atoms { /// Used by the atom class destructor: set atoms_refcount entry to zero virtual int clear_atom(int index); + /// Check that all calls use multiples of the same, use that number + virtual int set_atom_list_frequency(int atom_list_freq); + /// Clear atomic data int reset(); From 96e18ab8d5de864b13c4069c34f1136de80c2d49 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 15 Mar 2021 00:10:31 -0400 Subject: [PATCH 08/15] Use atom_list_freq to support timeStepFactor for CVCs --- src/colvar.cpp | 30 ++++++++++++++++++++++-------- src/colvar.h | 5 ++++- src/colvarcomp.cpp | 6 ++++++ 3 files changed, 32 insertions(+), 9 deletions(-) diff --git a/src/colvar.cpp b/src/colvar.cpp index bf83162bc..1fb2b0077 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -298,14 +298,7 @@ int colvar::init(std::string const &conf) reset_bias_force(); - get_keyval(conf, "timeStepFactor", time_step_factor, 1); - if (time_step_factor < 0) { - cvm::error("Error: timeStepFactor must be positive.\n"); - return COLVARS_ERROR; - } - if (time_step_factor != 1) { - enable(f_cv_multiple_ts); - } + error_code |= init_mts_parameters(conf); error_code |= init_grid_parameters(conf); @@ -510,6 +503,27 @@ int colvar::init_custom_function(std::string const &conf) #endif // #ifdef LEPTON +int colvar::init_mts_parameters(std::string const &conf) +{ + int error_code = COLVARS_OK; + get_keyval(conf, "timeStepFactor", time_step_factor, 1); + if (time_step_factor <= 0) { + return cvm::error("Error: timeStepFactor must be positive.\n", + COLVARS_INPUT_ERROR); + } else { + if (time_step_factor != 1) { + enable(f_cv_multiple_ts); + for (size_t i = 0; i < cvcs.size(); i++) { + // Tell CVCs that atom lists will be refreshed only every + // time_step_factor steps + error_code |= (cvcs[i])->set_atom_list_frequency(time_step_factor); + } + } + } + return error_code; +} + + int colvar::init_grid_parameters(std::string const &conf) { int error_code = COLVARS_OK; diff --git a/src/colvar.h b/src/colvar.h index 228d7967e..0ad2e8643 100644 --- a/src/colvar.h +++ b/src/colvar.h @@ -260,7 +260,10 @@ class colvar : public colvarparse, public colvardeps { /// Parse parameters for custom function with Lepton int init_custom_function(std::string const &conf); - /// Init defaults for grid options + /// Parse parameters for multiple-time stepping + int init_mts_parameters(std::string const &conf); + + /// Parse parameters for grid definitions int init_grid_parameters(std::string const &conf); /// Consistency check for the grid paramaters diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index f8fdd3f18..00bef34fa 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -147,6 +147,12 @@ int colvar::cvc::init(std::string const &conf) int colvar::cvc::set_atom_list_frequency(int new_frequency) { + if (atom_list_freq > 0) { + if (atom_list_freq != new_frequency) { + return cvm::error("Error: inconsistent definitions for " + "atomListFrequency.\n", COLVARS_INPUT_ERROR); + } + } atom_list_freq = new_frequency; return cvm::main()->proxy->set_atom_list_frequency(atom_list_freq); } From b3df2115946a0bb2354ecd0bbd019fdf1472ef16 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 15 Mar 2021 00:43:04 -0400 Subject: [PATCH 09/15] Fix value of timeStepFactor in regtest to reflect doc --- namd/tests/library/016_mts_bias/test.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/namd/tests/library/016_mts_bias/test.in b/namd/tests/library/016_mts_bias/test.in index 27e4a017d..b59b05c48 100644 --- a/namd/tests/library/016_mts_bias/test.in +++ b/namd/tests/library/016_mts_bias/test.in @@ -7,7 +7,7 @@ colvar { name one outputAppliedForce on - timeStepFactor 0 + timeStepFactor 4 width 0.5 From aee991e7c53fa3a6f479c922ca92d025a44e5424 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 15 Mar 2021 00:56:01 -0400 Subject: [PATCH 10/15] Add test for RMSD with timeStepFactor (explicit atoms + dynamic atom list combo) --- .../AutoDiff/test.colvars.out | 135 +++++++++++++++++ .../AutoDiff/test.colvars.state.stripped | 17 +++ .../AutoDiff/test.colvars.traj | 22 +++ .../AutoDiff/test.restart.colvars.out | 140 ++++++++++++++++++ .../test.restart.colvars.state.stripped | 17 +++ .../AutoDiff/test.restart.colvars.traj | 22 +++ .../namd-version.txt | 3 + .../000_rmsd-mts_harmonic-fixed/test.in | 28 ++++ tests/build_tests.sh | 2 +- tests/rmsd-mts.in | 17 +++ 10 files changed, 402 insertions(+), 1 deletion(-) create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.out create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.state.stripped create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.traj create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.out create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.state.stripped create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.traj create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/namd-version.txt create mode 100644 namd/tests/library/000_rmsd-mts_harmonic-fixed/test.in create mode 100644 tests/rmsd-mts.in diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.out b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.out new file mode 100644 index 000000000..07089ff51 --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.out @@ -0,0 +1,135 @@ +colvars: ---------------------------------------------------------------------- +colvars: Please cite Fiorin et al, Mol Phys 2013: +colvars: https://dx.doi.org/10.1080/00268976.2013.813594 +colvars: in any publication based on this calculation. +colvars: SMP parallelism is enabled; if needed, use "smp off" to override this. +colvars: This version was built with the C++11 standard or higher. +colvars: ---------------------------------------------------------------------- +colvars: Reading new configuration from file "test.in": +colvars: # units = "" [default] +colvars: # indexFile = "index.ndx" +colvars: The following index groups are currently defined: +colvars: Protein (104 atoms) +colvars: Protein_noH (51 atoms) +colvars: Protein_Backbone (40 atoms) +colvars: Protein_C-alpha (10 atoms) +colvars: RMSD_atoms (10 atoms) +colvars: Protein_C-alpha_1_2 (2 atoms) +colvars: Protein_C-alpha_9_10 (2 atoms) +colvars: Protein_C-alpha_1 (1 atoms) +colvars: group1 (4 atoms) +colvars: Protein_C-alpha_2 (1 atoms) +colvars: group2 (4 atoms) +colvars: Protein_C-alpha_3 (1 atoms) +colvars: group3 (4 atoms) +colvars: Protein_C-alpha_4 (1 atoms) +colvars: group4 (4 atoms) +colvars: Protein_C-alpha_5 (1 atoms) +colvars: group5 (4 atoms) +colvars: Protein_C-alpha_6 (1 atoms) +colvars: group6 (4 atoms) +colvars: Protein_C-alpha_7 (1 atoms) +colvars: group7 (4 atoms) +colvars: Protein_C-alpha_8 (1 atoms) +colvars: group8 (4 atoms) +colvars: Protein_C-alpha_9 (1 atoms) +colvars: group9 (4 atoms) +colvars: Protein_C-alpha_10 (1 atoms) +colvars: group10 (4 atoms) +colvars: heavy_atoms (51 atoms) +colvars: # smp = on [default] +colvars: # colvarsTrajFrequency = 1 +colvars: # colvarsRestartFrequency = 10 +colvars: # scriptedColvarForces = off [default] +colvars: # scriptingAfterBiases = off [default] +colvars: ---------------------------------------------------------------------- +colvars: Initializing a new collective variable. +colvars: # name = "one" +colvars: Initializing a new "rmsd" component. +colvars: # name = "" [default] +colvars: # componentCoeff = 1 [default] +colvars: # componentExp = 1 [default] +colvars: # period = 0 [default] +colvars: # wrapAround = 0 [default] +colvars: # forceNoPBC = off [default] +colvars: # scalable = on [default] +colvars: Initializing atom group "atoms". +colvars: # name = "" [default] +colvars: # centerReference = off [default] +colvars: # rotateReference = off [default] +colvars: # atomsOfGroup = "" [default] +colvars: # indexGroup = "RMSD_atoms" +colvars: # psfSegID = [default] +colvars: # atomsFile = "" [default] +colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] +colvars: # enableForces = on [default] +colvars: # enableFitGradients = on [default] +colvars: # printAtomIDs = off [default] +colvars: Atom group "atoms" defined with 10 atoms requested: total mass = 120.11, total charge = 0.53. +colvars: # refPositions = [default] +colvars: # refPositionsFile = "rmsd_atoms_refpos.xyz" +colvars: # refPositionsCol = "" [default] +colvars: Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units. +colvars: Enabling "centerReference" and "rotateReference", to minimize RMSD before calculating it as a variable: if this is not the desired behavior, disable them explicitly within the "atoms" block. +colvars: This is a standard minimum RMSD, derivatives of the optimal rotation will not be computed as they cancel out in the gradients. +colvars: All components initialized. +colvars: # timeStepFactor = 4 +colvars: # width = 0.5 +colvars: # lowerBoundary = 0 [default] +colvars: # upperBoundary = 0 [default] +colvars: # hardLowerBoundary = on [default] +colvars: # hardUpperBoundary = off [default] +colvars: # expandBoundaries = off [default] +colvars: # extendedLagrangian = off [default] +colvars: # outputValue = on [default] +colvars: # outputVelocity = off [default] +colvars: # outputTotalForce = off [default] +colvars: # outputAppliedForce = on +colvars: # subtractAppliedForce = off [default] +colvars: # runAve = off [default] +colvars: # corrFunc = off [default] +colvars: ---------------------------------------------------------------------- +colvars: Collective variables initialized, 1 in total. +colvars: ---------------------------------------------------------------------- +colvars: Initializing a new "harmonic" instance. +colvars: # name = "harmonic1" [default] +colvars: # colvars = { one } +colvars: # stepZeroData = off [default] +colvars: # outputEnergy = off [default] +colvars: # outputFreq = 10 [default] +colvars: # timeStepFactor = 4 +colvars: # writeTISamples = off [default] +colvars: # writeTIPMF = off [default] +colvars: # centers = { 0.1 } +colvars: # targetCenters = { 0.1 } [default] +colvars: # outputCenters = off [default] +colvars: # forceConstant = 0.001 +colvars: # targetForceConstant = -1 [default] +colvars: The force constant for colvar "one" will be rescaled to 0.004 according to the specified width (0.5). +colvars: ---------------------------------------------------------------------- +colvars: Collective variables biases initialized, 1 in total. +colvars: ---------------------------------------------------------------------- +colvars: Collective variables module (re)initialized. +colvars: ---------------------------------------------------------------------- +colvars: Updating NAMD interface: +colvars: updating atomic data (10 atoms). +colvars: updating group data (0 scalable groups, 0 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: Re-initialized atom group for variable "one":0/0. 10 atoms: total mass = 120.11, total charge = 0.53. +colvars: The restart output state file will be "test.tmp.colvars.state". +colvars: The final output state file will be "test.colvars.state". +colvars: Opening trajectory file "test.colvars.traj". +colvars: Redefining the Tcl "cv" command to the new script interface. +colvars: Updating NAMD interface: +colvars: updating atomic data (10 atoms). +colvars: updating group data (0 scalable groups, 0 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: Re-initialized atom group for variable "one":0/0. 10 atoms: total mass = 120.11, total charge = 0.53. +colvars: The restart output state file will be "test.tmp.colvars.state". +colvars: The final output state file will be "test.colvars.state". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". +colvars: Saving collective variables state to "test.tmp.colvars.state". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". +colvars: Saving collective variables state to "test.tmp.colvars.state". +colvars: Saving collective variables state to "test.colvars.state". diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.state.stripped b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.state.stripped new file mode 100644 index 000000000..afb0eb97f --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.state.stripped @@ -0,0 +1,17 @@ +configuration { + step 20 + dt 1.000000e+00 +} + +colvar { + name one + x 6.34780007519980e-01 +} + +restraint { + configuration { + step 20 + name harmonic1 + } +} + diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.traj b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.traj new file mode 100644 index 000000000..99efeb74f --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.colvars.traj @@ -0,0 +1,22 @@ +# step one fa_one + 0 6.13624515542194e-01 -8.21799224867510e-03 + 1 6.13624515542194e-01 0.00000000000000e+00 + 2 6.13624515542194e-01 0.00000000000000e+00 + 3 6.13624515542194e-01 0.00000000000000e+00 + 4 6.08443936340912e-01 -8.13510298145459e-03 + 5 6.08443936340912e-01 0.00000000000000e+00 + 6 6.08443936340912e-01 0.00000000000000e+00 + 7 6.08443936340912e-01 0.00000000000000e+00 + 8 6.05743810280055e-01 -8.09190096448088e-03 + 9 6.05743810280055e-01 0.00000000000000e+00 + 10 6.05743810280055e-01 0.00000000000000e+00 + 11 6.05743810280055e-01 0.00000000000000e+00 + 12 6.10959614033100e-01 -8.17535382452960e-03 + 13 6.10959614033100e-01 0.00000000000000e+00 + 14 6.10959614033100e-01 0.00000000000000e+00 + 15 6.10959614033100e-01 0.00000000000000e+00 + 16 6.23032917931524e-01 -8.36852668690438e-03 + 17 6.23032917931524e-01 0.00000000000000e+00 + 18 6.23032917931524e-01 0.00000000000000e+00 + 19 6.23032917931524e-01 0.00000000000000e+00 + 20 6.34780007519980e-01 -8.55648012031968e-03 diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.out b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.out new file mode 100644 index 000000000..5ba5e62c9 --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.out @@ -0,0 +1,140 @@ +colvars: ---------------------------------------------------------------------- +colvars: Please cite Fiorin et al, Mol Phys 2013: +colvars: https://dx.doi.org/10.1080/00268976.2013.813594 +colvars: in any publication based on this calculation. +colvars: SMP parallelism is enabled; if needed, use "smp off" to override this. +colvars: This version was built with the C++11 standard or higher. +colvars: ---------------------------------------------------------------------- +colvars: Reading new configuration from file "test.in": +colvars: # units = "" [default] +colvars: # indexFile = "index.ndx" +colvars: The following index groups are currently defined: +colvars: Protein (104 atoms) +colvars: Protein_noH (51 atoms) +colvars: Protein_Backbone (40 atoms) +colvars: Protein_C-alpha (10 atoms) +colvars: RMSD_atoms (10 atoms) +colvars: Protein_C-alpha_1_2 (2 atoms) +colvars: Protein_C-alpha_9_10 (2 atoms) +colvars: Protein_C-alpha_1 (1 atoms) +colvars: group1 (4 atoms) +colvars: Protein_C-alpha_2 (1 atoms) +colvars: group2 (4 atoms) +colvars: Protein_C-alpha_3 (1 atoms) +colvars: group3 (4 atoms) +colvars: Protein_C-alpha_4 (1 atoms) +colvars: group4 (4 atoms) +colvars: Protein_C-alpha_5 (1 atoms) +colvars: group5 (4 atoms) +colvars: Protein_C-alpha_6 (1 atoms) +colvars: group6 (4 atoms) +colvars: Protein_C-alpha_7 (1 atoms) +colvars: group7 (4 atoms) +colvars: Protein_C-alpha_8 (1 atoms) +colvars: group8 (4 atoms) +colvars: Protein_C-alpha_9 (1 atoms) +colvars: group9 (4 atoms) +colvars: Protein_C-alpha_10 (1 atoms) +colvars: group10 (4 atoms) +colvars: heavy_atoms (51 atoms) +colvars: # smp = on [default] +colvars: # colvarsTrajFrequency = 1 +colvars: # colvarsRestartFrequency = 10 +colvars: # scriptedColvarForces = off [default] +colvars: # scriptingAfterBiases = off [default] +colvars: ---------------------------------------------------------------------- +colvars: Initializing a new collective variable. +colvars: # name = "one" +colvars: Initializing a new "rmsd" component. +colvars: # name = "" [default] +colvars: # componentCoeff = 1 [default] +colvars: # componentExp = 1 [default] +colvars: # period = 0 [default] +colvars: # wrapAround = 0 [default] +colvars: # forceNoPBC = off [default] +colvars: # scalable = on [default] +colvars: Initializing atom group "atoms". +colvars: # name = "" [default] +colvars: # centerReference = off [default] +colvars: # rotateReference = off [default] +colvars: # atomsOfGroup = "" [default] +colvars: # indexGroup = "RMSD_atoms" +colvars: # psfSegID = [default] +colvars: # atomsFile = "" [default] +colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] +colvars: # enableForces = on [default] +colvars: # enableFitGradients = on [default] +colvars: # printAtomIDs = off [default] +colvars: Atom group "atoms" defined with 10 atoms requested: total mass = 120.11, total charge = 0.53. +colvars: # refPositions = [default] +colvars: # refPositionsFile = "rmsd_atoms_refpos.xyz" +colvars: # refPositionsCol = "" [default] +colvars: Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units. +colvars: Enabling "centerReference" and "rotateReference", to minimize RMSD before calculating it as a variable: if this is not the desired behavior, disable them explicitly within the "atoms" block. +colvars: This is a standard minimum RMSD, derivatives of the optimal rotation will not be computed as they cancel out in the gradients. +colvars: All components initialized. +colvars: # timeStepFactor = 4 +colvars: # width = 0.5 +colvars: # lowerBoundary = 0 [default] +colvars: # upperBoundary = 0 [default] +colvars: # hardLowerBoundary = on [default] +colvars: # hardUpperBoundary = off [default] +colvars: # expandBoundaries = off [default] +colvars: # extendedLagrangian = off [default] +colvars: # outputValue = on [default] +colvars: # outputVelocity = off [default] +colvars: # outputTotalForce = off [default] +colvars: # outputAppliedForce = on +colvars: # subtractAppliedForce = off [default] +colvars: # runAve = off [default] +colvars: # corrFunc = off [default] +colvars: ---------------------------------------------------------------------- +colvars: Collective variables initialized, 1 in total. +colvars: ---------------------------------------------------------------------- +colvars: Initializing a new "harmonic" instance. +colvars: # name = "harmonic1" [default] +colvars: # colvars = { one } +colvars: # stepZeroData = off [default] +colvars: # outputEnergy = off [default] +colvars: # outputFreq = 10 [default] +colvars: # timeStepFactor = 4 +colvars: # writeTISamples = off [default] +colvars: # writeTIPMF = off [default] +colvars: # centers = { 0.1 } +colvars: # targetCenters = { 0.1 } [default] +colvars: # outputCenters = off [default] +colvars: # forceConstant = 0.001 +colvars: # targetForceConstant = -1 [default] +colvars: The force constant for colvar "one" will be rescaled to 0.004 according to the specified width (0.5). +colvars: ---------------------------------------------------------------------- +colvars: Collective variables biases initialized, 1 in total. +colvars: ---------------------------------------------------------------------- +colvars: Collective variables module (re)initialized. +colvars: ---------------------------------------------------------------------- +colvars: Updating NAMD interface: +colvars: updating atomic data (10 atoms). +colvars: updating group data (0 scalable groups, 0 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: Re-initialized atom group for variable "one":0/0. 10 atoms: total mass = 120.11, total charge = 0.53. +colvars: ---------------------------------------------------------------------- +colvars: Loading state from file "test.colvars.state". +colvars: Restarting collective variable "one" from value: 0.63478 +colvars: Restarting harmonic bias "harmonic1" from step number 20. +colvars: ---------------------------------------------------------------------- +colvars: The restart output state file will be "test.restart.tmp.colvars.state". +colvars: The final output state file will be "test.restart.colvars.state". +colvars: Opening trajectory file "test.restart.colvars.traj". +colvars: Redefining the Tcl "cv" command to the new script interface. +colvars: Updating NAMD interface: +colvars: updating atomic data (10 atoms). +colvars: updating group data (0 scalable groups, 0 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: Re-initialized atom group for variable "one":0/0. 10 atoms: total mass = 120.11, total charge = 0.53. +colvars: The restart output state file will be "test.restart.tmp.colvars.state". +colvars: The final output state file will be "test.restart.colvars.state". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj". +colvars: Saving collective variables state to "test.restart.tmp.colvars.state". +colvars: Synchronizing (emptying the buffer of) trajectory file "test.restart.colvars.traj". +colvars: Saving collective variables state to "test.restart.tmp.colvars.state". +colvars: Saving collective variables state to "test.restart.colvars.state". diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.state.stripped b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.state.stripped new file mode 100644 index 000000000..4f859d49c --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.state.stripped @@ -0,0 +1,17 @@ +configuration { + step 40 + dt 1.000000e+00 +} + +colvar { + name one + x 6.57869113123593e-01 +} + +restraint { + configuration { + step 40 + name harmonic1 + } +} + diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.traj b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.traj new file mode 100644 index 000000000..1ce74e8bc --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/AutoDiff/test.restart.colvars.traj @@ -0,0 +1,22 @@ +# step one fa_one + 20 6.34780007519980e-01 -8.55648012031968e-03 + 21 6.34780007519980e-01 0.00000000000000e+00 + 22 6.34780007519980e-01 0.00000000000000e+00 + 23 6.34780007519980e-01 0.00000000000000e+00 + 24 6.41589733432378e-01 -8.66543573491804e-03 + 25 6.41589733432378e-01 0.00000000000000e+00 + 26 6.41589733432378e-01 0.00000000000000e+00 + 27 6.41589733432378e-01 0.00000000000000e+00 + 28 6.45047786395882e-01 -8.72076458233411e-03 + 29 6.45047786395882e-01 0.00000000000000e+00 + 30 6.45047786395882e-01 0.00000000000000e+00 + 31 6.45047786395882e-01 0.00000000000000e+00 + 32 6.48532740533488e-01 -8.77652384853580e-03 + 33 6.48532740533488e-01 0.00000000000000e+00 + 34 6.48532740533488e-01 0.00000000000000e+00 + 35 6.48532740533488e-01 0.00000000000000e+00 + 36 6.53382083727100e-01 -8.85411333963360e-03 + 37 6.53382083727100e-01 0.00000000000000e+00 + 38 6.53382083727100e-01 0.00000000000000e+00 + 39 6.53382083727100e-01 0.00000000000000e+00 + 40 6.57869113123593e-01 -8.92590580997749e-03 diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/namd-version.txt b/namd/tests/library/000_rmsd-mts_harmonic-fixed/namd-version.txt new file mode 100644 index 000000000..617210efd --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/namd-version.txt @@ -0,0 +1,3 @@ +Info: NAMD 2.14 for Linux-x86_64-multicore +colvars: Initializing the collective variables module, version "2020-07-07". +colvars: Using NAMD interface, version "2020-05-04". diff --git a/namd/tests/library/000_rmsd-mts_harmonic-fixed/test.in b/namd/tests/library/000_rmsd-mts_harmonic-fixed/test.in new file mode 100644 index 000000000..e51888e79 --- /dev/null +++ b/namd/tests/library/000_rmsd-mts_harmonic-fixed/test.in @@ -0,0 +1,28 @@ +colvarsTrajFrequency 1 +colvarsRestartFrequency 10 +indexFile index.ndx + +colvar { + + name one + + outputAppliedForce on + + width 0.5 + + timeStepFactor 4 + + rmsd { + atoms { + indexGroup RMSD_atoms + } + refPositionsFile rmsd_atoms_refpos.xyz + } +} + +harmonic { + colvars one + centers 0.1 + forceConstant 0.001 + timeStepFactor 4 +} diff --git a/tests/build_tests.sh b/tests/build_tests.sh index 5020e77d0..40e0efa37 100755 --- a/tests/build_tests.sh +++ b/tests/build_tests.sh @@ -218,7 +218,7 @@ for colvar in \ "hbond" \ "inertia" \ "inertiaz" \ - "rmsd" \ + "rmsd" "rmsd-mts" \ "eigenvector" \ "eigenvector-difference" \ "eigenvector-normalized" \ diff --git a/tests/rmsd-mts.in b/tests/rmsd-mts.in new file mode 100644 index 000000000..ddd9d77d9 --- /dev/null +++ b/tests/rmsd-mts.in @@ -0,0 +1,17 @@ +colvar { + + name one + + outputAppliedForce on + + width 0.5 + + timeStepFactor 4 + + rmsd { + atoms { + indexGroup RMSD_atoms + } + refPositionsFile rmsd_atoms_refpos.xyz + } +} From ec3f7087d0548ef79b7a35feb69086f05364f6d0 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 26 Jan 2023 16:25:34 -0500 Subject: [PATCH 11/15] Silence unused-variable warnings --- src/colvar.cpp | 9 ++++----- src/colvarcomp.cpp | 2 +- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/colvar.cpp b/src/colvar.cpp index 1fb2b0077..27b1c2582 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -73,10 +73,10 @@ int colvar::init(std::string const &conf) if ((cvm::colvar_by_name(this->name) != NULL) && (cvm::colvar_by_name(this->name) != this)) { - cvm::error("Error: this colvar cannot have the same name, \""+this->name+ - "\", as another colvar.\n", - COLVARS_INPUT_ERROR); - return COLVARS_INPUT_ERROR; + error_code |= cvm::error("Error: this colvar cannot have the same name, \""+ + this->name+ + "\", as another colvar.\n", + COLVARS_INPUT_ERROR); } // Initialize dependency members @@ -1756,7 +1756,6 @@ int colvar::collect_cvc_Jacobians() int colvar::update_requested_atoms() { int error_code = COLVARS_OK; - size_t i, cvc_count; for (size_t i = 0; i < cvcs.size(); i++) { error_code |= (cvcs[i])->update_all_requested_atoms(); } diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 00bef34fa..4435d373d 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -196,7 +196,7 @@ int colvar::cvc::update_all_requested_atoms() agi != atom_groups.end(); agi++) { error_code |= update_requested_atoms(*agi); } - return COLVARS_OK; + return error_code; } From d221a0277d521fc8cd6662be6a27d6f63d82c5be Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Apr 2023 21:54:37 -0400 Subject: [PATCH 12/15] Use step number when reporting NAMD properties in log --- namd/tests/library/Common/measure_net_force_torque.tcl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/namd/tests/library/Common/measure_net_force_torque.tcl b/namd/tests/library/Common/measure_net_force_torque.tcl index 2ebf337d9..7b35eb452 100644 --- a/namd/tests/library/Common/measure_net_force_torque.tcl +++ b/namd/tests/library/Common/measure_net_force_torque.tcl @@ -36,11 +36,11 @@ tclForcesScript { } } # print "[expr $t-1] Force: ( $tot ) Torque: ( $torque )" - print "NUM. ATOMS REQUESTED: [cv getnumatoms]" - print "NUM. ATOMS USED: [cv getnumactiveatoms]" - print "NUM. ATOM GROUPS USED: [cv getnumactiveatomgroups]" - print [format "TOTAL FORCE: \{ %.16f %.16f %.16f \}" {*}$tot] - print [format "TOTAL TORQUE: \{ %.16f %.16f %.16f \}" {*}$torque] + print "STEP $t: NUM. ATOMS: [cv getnumatoms] (TOTAL)" + print "STEP $t: NUM. ATOMS: [cv getnumactiveatoms] (ACTIVE)" + print "STEP $t: NUM. ATOM GROUPS: [cv getnumactiveatomgroups]" + print [format "STEP $t: TOTAL FORCE: \{ %.16f %.16f %.16f \}" {*}$tot] + print [format "STEP $t: TOTAL TORQUE: \{ %.16f %.16f %.16f \}" {*}$torque] array set prev_c [array get c] } incr t From 93e29eaefb34a08e8698d6ebb0e79fdbbea263eb Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Thu, 27 Aug 2020 18:24:19 -0400 Subject: [PATCH 13/15] Add support for reupdating the atom list in the proxy --- src/colvarproxy.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/colvarproxy.h b/src/colvarproxy.h index c4a09e094..3f3cce803 100644 --- a/src/colvarproxy.h +++ b/src/colvarproxy.h @@ -61,9 +61,8 @@ class colvarproxy_atoms { /// Check whether it is possible to select atoms by residue number name virtual int check_atom_name_selections_available(); - /// Select this atom for collective variables calculation, using name and - /// residue number. Not all programs support this: leave this function as - /// is in those cases. + /// Prepare this atom using name and residue number (not all engines support + /// it, and this function does not necessarily need reimplementation. virtual int init_atom(cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id); From 2f43d648770cee8be7c8248f0be5053870051459 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Fri, 4 Apr 2025 17:46:51 -0400 Subject: [PATCH 14/15] Small fix --- devel-tools/compile-vmd.sh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/devel-tools/compile-vmd.sh b/devel-tools/compile-vmd.sh index 8cfba9b0c..6073e316c 100755 --- a/devel-tools/compile-vmd.sh +++ b/devel-tools/compile-vmd.sh @@ -246,6 +246,8 @@ compile_vmd_target() { rm -f plugins ln -s ${VMDPLUGINDIR} plugins + retcode=1 + if VMDINSTALLBINDIR=${VMDINSTALLBINDIR} \ VMDINSTALLLIBRARYDIR=${VMDINSTALLLIBRARYDIR} \ ./configure LINUXAMD64 ${VMD_OPTS[@]} ; then From 5731262485e6213dbf071518d83d43c23bdec4d6 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Fri, 4 Apr 2025 17:47:08 -0400 Subject: [PATCH 15/15] Fix merge error --- vmd/src/colvarproxy_vmd.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vmd/src/colvarproxy_vmd.C b/vmd/src/colvarproxy_vmd.C index d6b425bdb..942628f33 100644 --- a/vmd/src/colvarproxy_vmd.C +++ b/vmd/src/colvarproxy_vmd.C @@ -886,11 +886,11 @@ int colvarproxy_vmd::compute_volmap(int flags, int const new_flags = volmap_flag_gradients | volmap_flag_use_atom_field; compute_voldata(voldata, atom_begin, atom_end, - value, atom_field); + value, atom_field, NULL); } else { int const new_flags = volmap_flag_gradients; compute_voldata(voldata, atom_begin, atom_end, - value, NULL); + value, NULL, NULL); } } else { @@ -898,11 +898,11 @@ int colvarproxy_vmd::compute_volmap(int flags, if (flags & volmap_flag_use_atom_field) { int const new_flags = volmap_flag_use_atom_field; compute_voldata(voldata, atom_begin, atom_end, - value, atom_field); + value, atom_field, NULL); } else { int const new_flags = volmap_flag_null; compute_voldata(voldata, atom_begin, atom_end, - value, NULL); + value, NULL, NULL); } } } else {