From 212a955285830e2957efcd2b0421eced15a0cf9c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 15 Apr 2016 16:07:01 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14829 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/Makefile.fermi | 2 +- lib/colvars/Makefile.g++ | 2 +- lib/colvars/Makefile.mingw32-cross | 2 +- lib/colvars/Makefile.mingw64-cross | 2 +- lib/colvars/colvar.cpp | 917 ++++++++++++--------------- lib/colvars/colvar.h | 162 ++--- lib/colvars/colvaratoms.cpp | 141 ++-- lib/colvars/colvaratoms.h | 16 +- lib/colvars/colvarbias.cpp | 23 +- lib/colvars/colvarbias.h | 13 +- lib/colvars/colvarbias_abf.cpp | 41 +- lib/colvars/colvarbias_abf.h | 2 +- lib/colvars/colvarbias_alb.cpp | 6 +- lib/colvars/colvarbias_alb.h | 2 +- lib/colvars/colvarbias_histogram.cpp | 14 +- lib/colvars/colvarbias_histogram.h | 2 +- lib/colvars/colvarbias_meta.cpp | 12 +- lib/colvars/colvarbias_meta.h | 4 +- lib/colvars/colvarbias_restraint.cpp | 11 +- lib/colvars/colvarbias_restraint.h | 2 +- lib/colvars/colvarcomp.cpp | 174 +++-- lib/colvars/colvarcomp.h | 83 ++- lib/colvars/colvarcomp_angles.cpp | 249 ++++---- lib/colvars/colvarcomp_coordnums.cpp | 63 +- lib/colvars/colvarcomp_distances.cpp | 587 +++++++++-------- lib/colvars/colvarcomp_rotations.cpp | 85 ++- lib/colvars/colvargrid.h | 9 - lib/colvars/colvarmodule.cpp | 412 ++++++++---- lib/colvars/colvarmodule.h | 81 ++- lib/colvars/colvarproxy.h | 75 ++- lib/colvars/colvarscript.cpp | 12 +- 31 files changed, 1713 insertions(+), 1493 deletions(-) diff --git a/lib/colvars/Makefile.fermi b/lib/colvars/Makefile.fermi index ba00a2e6d8..fe8d3ac16c 100644 --- a/lib/colvars/Makefile.fermi +++ b/lib/colvars/Makefile.fermi @@ -19,7 +19,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \ colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ + colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ colvarscript.cpp colvartypes.cpp colvarvalue.cpp LIB = libcolvars.a diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++ index 8f66b0c743..d0fdd1c325 100644 --- a/lib/colvars/Makefile.g++ +++ b/lib/colvars/Makefile.g++ @@ -18,7 +18,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \ colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ + colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ colvarscript.cpp colvartypes.cpp colvarvalue.cpp LIB = libcolvars.a diff --git a/lib/colvars/Makefile.mingw32-cross b/lib/colvars/Makefile.mingw32-cross index 381f4455d5..289b006a45 100644 --- a/lib/colvars/Makefile.mingw32-cross +++ b/lib/colvars/Makefile.mingw32-cross @@ -21,7 +21,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \ colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ + colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ colvarscript.cpp colvartypes.cpp colvarvalue.cpp DIR = Obj_mingw32/ diff --git a/lib/colvars/Makefile.mingw64-cross b/lib/colvars/Makefile.mingw64-cross index bcb23c6475..a6081273f8 100644 --- a/lib/colvars/Makefile.mingw64-cross +++ b/lib/colvars/Makefile.mingw64-cross @@ -21,7 +21,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \ colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ + colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \ colvarscript.cpp colvartypes.cpp colvarvalue.cpp DIR = Obj_mingw64/ diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index fa38c4d33f..ebc3bf44f1 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -32,10 +32,13 @@ colvar::colvar(std::string const &conf) return; } - // all tasks disabled by default - for (i = 0; i < task_ntot; i++) { - tasks[i] = false; - } + // Initialize dependency members + // Could be a function defined in a different source file, for space? + + this->description = "colvar " + this->name; + + // Initialize static array once and for all + init_cv_requires(); kinetic_energy = 0.0; potential_energy = 0.0; @@ -93,6 +96,7 @@ colvar::colvar(std::string const &conf) cvcs.back()->name = s.str(); \ /* pad cvc number for correct ordering when sorting by name */\ } \ + cvcs.back()->setup(); \ if (cvm::debug()) \ cvm::log("Done initializing a \""+ \ std::string(def_config_key)+ \ @@ -166,13 +170,22 @@ colvar::colvar(std::string const &conf) return; } + n_active_cvcs = cvcs.size(); + cvm::log("All components initialized.\n"); + // Store list of children cvcs for dependency checking purposes + for (i = 0; i < cvcs.size(); i++) { + add_child(cvcs[i]); + } + // Setup colvar as scripted function of components if (get_keyval(conf, "scriptedFunction", scripted_function, "", colvarparse::parse_silent)) { - enable(task_scripted); + // Make feature available only on user request + provide(f_cv_scripted); + enable(f_cv_scripted); cvm::log("This colvar uses scripted function \"" + scripted_function + "\"."); std::string type_str; @@ -219,16 +232,9 @@ colvar::colvar(std::string const &conf) for (i = 0; i < cvcs.size(); i++) { sorted_cvc_values.push_back(&(cvcs[i]->value())); } - - b_homogeneous = false; - // Scripted functions are deemed non-periodic - b_periodic = false; - period = 0.0; - b_inverse_gradients = false; - b_Jacobian_force = false; } - if (!tasks[task_scripted]) { + if (!is_enabled(f_cv_scripted)) { colvarvalue const &cvc_value = (cvcs[0])->value(); if (cvm::debug()) cvm::log ("This collective variable is a "+ @@ -242,51 +248,56 @@ colvar::colvar(std::string const &conf) // If using scripted biases, any colvar may receive bias forces // and will need its gradient if (cvm::scripted_forces()) { - enable(task_gradients); + enable(f_cv_gradient); } // check for linear combinations + { + bool lin = !is_enabled(f_cv_scripted); + for (i = 0; i < cvcs.size(); i++) { - b_linear = !tasks[task_scripted]; - for (i = 0; i < cvcs.size(); i++) { - if ((cvcs[i])->b_debug_gradients) - enable(task_gradients); + // FIXME this is a reverse dependency, ie. cv feature depends on cvc flag + // need to clarify this case + // if ((cvcs[i])->b_debug_gradients) + // enable(task_gradients); - if ((cvcs[i])->sup_np != 1) { - if (cvm::debug() && b_linear) - cvm::log("Warning: You are using a non-linear polynomial " - "combination to define this collective variable, " - "some biasing methods may be unavailable.\n"); - b_linear = false; + if ((cvcs[i])->sup_np != 1) { + if (cvm::debug() && lin) + cvm::log("Warning: You are using a non-linear polynomial " + "combination to define this collective variable, " + "some biasing methods may be unavailable.\n"); + lin = false; - if ((cvcs[i])->sup_np < 0) { - cvm::log("Warning: you chose a negative exponent in the combination; " - "if you apply forces, the simulation may become unstable " - "when the component \""+ - (cvcs[i])->function_type+"\" approaches zero.\n"); + if ((cvcs[i])->sup_np < 0) { + cvm::log("Warning: you chose a negative exponent in the combination; " + "if you apply forces, the simulation may become unstable " + "when the component \""+ + (cvcs[i])->function_type+"\" approaches zero.\n"); + } } } + feature_states[f_cv_linear]->enabled = lin; } // Colvar is homogeneous iff: - // - it is not scripted - // - it is linear + // - it is linear (hence not scripted) // - all cvcs have coefficient 1 or -1 // i.e. sum or difference of cvcs - - b_homogeneous = !tasks[task_scripted] && b_linear; - for (i = 0; i < cvcs.size(); i++) { - if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { - b_homogeneous = false; + { + bool homogeneous = is_enabled(f_cv_linear); + for (i = 0; i < cvcs.size(); i++) { + if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { + homogeneous = false; + } } + feature_states[f_cv_homogeneous]->enabled = homogeneous; } - // Colvar is deemed periodic iff: // - it is homogeneous // - all cvcs are periodic // - all cvcs have the same period - b_periodic = cvcs[0]->b_periodic && b_homogeneous; + b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous); period = cvcs[0]->period; for (i = 1; i < cvcs.size(); i++) { if (!cvcs[i]->b_periodic || cvcs[i]->period != period) { @@ -294,12 +305,10 @@ colvar::colvar(std::string const &conf) period = 0.0; } } + feature_states[f_cv_periodic]->enabled = b_periodic; - // these will be set to false if any of the cvcs has them false - b_inverse_gradients = !tasks[task_scripted]; - b_Jacobian_force = !tasks[task_scripted]; + // check that cvcs are compatible - // check the available features of each cvc for (i = 0; i < cvcs.size(); i++) { if ((cvcs[i])->b_periodic && !b_periodic) { cvm::log("Warning: although this component is periodic, the colvar will " @@ -308,14 +317,8 @@ colvar::colvar(std::string const &conf) "you know what you are doing!"); } - if (! (cvcs[i])->b_inverse_gradients) - b_inverse_gradients = false; - - if (! (cvcs[i])->b_Jacobian_derivative) - b_Jacobian_force = false; - // components may have different types only for scripted functions - if (!tasks[task_scripted] && (colvarvalue::check_types(cvcs[i]->value(), + if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(), cvcs[0]->value())) ) { cvm::error("ERROR: you are definining this collective variable " "by using components of different types. " @@ -325,6 +328,11 @@ colvar::colvar(std::string const &conf) } } + active_cvc_square_norm = 0.; + for (i = 0; i < cvcs.size(); i++) { + active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff; + } + // at this point, the colvar's type is defined f.type(value()); fb.type(value()); @@ -334,44 +342,51 @@ colvar::colvar(std::string const &conf) cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); } + // NOTE: not porting wall stuff to new deps, as this will change to a separate bias + // the grid functions will wait a little as well + lower_boundary.type(value()); lower_wall.type(value()); upper_boundary.type(value()); upper_wall.type(value()); - if (value().type() == colvarvalue::type_scalar) { + feature_states[f_cv_scalar]->enabled = (value().type() == colvarvalue::type_scalar); + + if (is_enabled(f_cv_scalar)) { if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) { - enable(task_lower_boundary); + provide(f_cv_lower_boundary); + enable(f_cv_lower_boundary); } get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0); if (lower_wall_k > 0.0) { get_keyval(conf, "lowerWall", lower_wall, lower_boundary); - enable(task_lower_wall); + enable(f_cv_lower_wall); } if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) { - enable(task_upper_boundary); + provide(f_cv_upper_boundary); + enable(f_cv_upper_boundary); } get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0); if (upper_wall_k > 0.0) { get_keyval(conf, "upperWall", upper_wall, upper_boundary); - enable(task_upper_wall); + enable(f_cv_upper_wall); } } - if (tasks[task_lower_boundary]) { + if (is_enabled(f_cv_lower_boundary)) { get_keyval(conf, "hardLowerBoundary", hard_lower_boundary, false); } - if (tasks[task_upper_boundary]) { + if (is_enabled(f_cv_upper_boundary)) { get_keyval(conf, "hardUpperBoundary", hard_upper_boundary, false); } // consistency checks for boundaries and walls - if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) { + if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) { if (lower_boundary >= upper_boundary) { cvm::error("Error: the upper boundary, "+ cvm::to_str(upper_boundary)+ @@ -381,7 +396,7 @@ colvar::colvar(std::string const &conf) } } - if (tasks[task_lower_wall] && tasks[task_upper_wall]) { + if (is_enabled(f_cv_lower_wall) && is_enabled(f_cv_upper_wall)) { if (lower_wall >= upper_wall) { cvm::error("Error: the upper wall, "+ cvm::to_str(upper_wall)+ @@ -389,13 +404,6 @@ colvar::colvar(std::string const &conf) cvm::to_str(lower_wall)+".\n", INPUT_ERROR); } - - if (dist2(lower_wall, upper_wall) < 1.0E-12) { - cvm::log("Lower wall and upper wall are equal " - "in the periodic domain of the colvar: disabling walls.\n"); - disable(task_lower_wall); - disable(task_upper_wall); - } } get_keyval(conf, "expandBoundaries", expand_boundaries, false); @@ -412,16 +420,18 @@ colvar::colvar(std::string const &conf) } { - bool b_extended_lagrangian; - get_keyval(conf, "extendedLagrangian", b_extended_lagrangian, false); + bool b_extended_Lagrangian; + get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false); - if (b_extended_lagrangian) { + if (b_extended_Lagrangian) { cvm::real temp, tolerance, period; cvm::log("Enabling the extended Lagrangian term for colvar \""+ this->name+"\".\n"); - enable(task_extended_lagrangian); + // Make feature available only on user request + provide(f_cv_extended_Lagrangian); + enable(f_cv_extended_Lagrangian); xr.type(value()); vr.type(value()); @@ -456,7 +466,7 @@ colvar::colvar(std::string const &conf) bool b_output_energy; get_keyval(conf, "outputEnergy", b_output_energy, false); if (b_output_energy) { - enable(task_output_energy); + enable(f_cv_output_energy); } } @@ -465,7 +475,7 @@ colvar::colvar(std::string const &conf) cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR); } if (ext_gamma != 0.0) { - enable(task_langevin); + enable(f_cv_Langevin); ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1 ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt()); } @@ -476,7 +486,7 @@ colvar::colvar(std::string const &conf) bool b_output_value; get_keyval(conf, "outputValue", b_output_value, true); if (b_output_value) { - enable(task_output_value); + enable(f_cv_output_value); } } @@ -484,7 +494,7 @@ colvar::colvar(std::string const &conf) bool b_output_velocity; get_keyval(conf, "outputVelocity", b_output_velocity, false); if (b_output_velocity) { - enable(task_output_velocity); + enable(f_cv_output_velocity); } } @@ -492,7 +502,7 @@ colvar::colvar(std::string const &conf) bool b_output_system_force; get_keyval(conf, "outputSystemForce", b_output_system_force, false); if (b_output_system_force) { - enable(task_output_system_force); + enable(f_cv_output_system_force); } } @@ -500,10 +510,22 @@ colvar::colvar(std::string const &conf) bool b_output_applied_force; get_keyval(conf, "outputAppliedForce", b_output_applied_force, false); if (b_output_applied_force) { - enable(task_output_applied_force); + enable(f_cv_output_applied_force); } } + // Start in active state by default + enable(f_cv_active); + // Make sure dependency side-effects are correct + refresh_deps(); + + x_old.type(value()); + v_fdiff.type(value()); + v_reported.type(value()); + fj.type(value()); + ft.type(value()); + ft_reported.type(value()); + if (cvm::b_analysis) parse_analysis(conf); @@ -512,6 +534,18 @@ colvar::colvar(std::string const &conf) } +int colvar::refresh_deps() +{ + // If enabled features are changed upstream, the features below should be refreshed + if (is_enabled(f_cv_system_force_calc)) { + cvm::request_system_force(); + } + if (is_enabled(f_cv_collect_gradient) && atom_ids.size() == 0) { + build_atom_list(); + } + return COLVARS_OK; +} + void colvar::build_atom_list(void) { @@ -520,8 +554,8 @@ void colvar::build_atom_list(void) for (size_t i = 0; i < cvcs.size(); i++) { for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { - for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + for (size_t k = 0; k < ag.size(); k++) { temp_id_list.push_back(ag[k].id); } } @@ -561,7 +595,7 @@ int colvar::parse_analysis(std::string const &conf) bool b_runave = false; if (get_keyval(conf, "runAve", b_runave) && b_runave) { - enable(task_runave); + enable(f_cv_runave); get_keyval(conf, "runAveLength", runave_length, 1000); get_keyval(conf, "runAveStride", runave_stride, 1); @@ -576,6 +610,7 @@ int colvar::parse_analysis(std::string const &conf) this->name+".runave.traj")); size_t const this_cv_width = x.output_width(cvm::cv_width); + cvm::backup_file(runave_outfile.c_str()); runave_os.open(runave_outfile.c_str()); runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " " @@ -589,7 +624,7 @@ int colvar::parse_analysis(std::string const &conf) bool b_acf = false; if (get_keyval(conf, "corrFunc", b_acf) && b_acf) { - enable(task_corrfunc); + enable(f_cv_corrfunc); std::string acf_colvar_name; get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name); @@ -606,9 +641,9 @@ int colvar::parse_analysis(std::string const &conf) acf_type = acf_coor; } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) { acf_type = acf_vel; - enable(task_fdiff_velocity); + enable(f_cv_fdiff_velocity); if (acf_colvar_name.size()) - (cvm::colvar_by_name(acf_colvar_name))->enable(task_fdiff_velocity); + (cvm::colvar_by_name(acf_colvar_name))->enable(f_cv_fdiff_velocity); } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) { acf_type = acf_p2coor; } else { @@ -634,199 +669,6 @@ int colvar::parse_analysis(std::string const &conf) } -int colvar::enable(colvar::task const &t) -{ - switch (t) { - - case task_output_system_force: - enable(task_system_force); - tasks[t] = true; - break; - - case task_report_Jacobian_force: - enable(task_Jacobian_force); - enable(task_system_force); - if (cvm::debug()) - cvm::log("Adding the Jacobian force to the system force, " - "rather than applying its opposite silently.\n"); - tasks[t] = true; - break; - - case task_Jacobian_force: - // checks below do not apply to extended-system colvars - if ( !tasks[task_extended_lagrangian] ) { - enable(task_gradients); - - if (!b_Jacobian_force) { - cvm::error("Error: colvar \""+this->name+ - "\" does not have Jacobian forces implemented.\n", - INPUT_ERROR); - } - if (!b_linear) { - cvm::error("Error: colvar \""+this->name+ - "\" must be defined as a linear combination " - "to calculate the Jacobian force.\n", - INPUT_ERROR); - } - if (cvm::debug()) - cvm::log("Enabling calculation of the Jacobian force " - "on this colvar.\n"); - } - fj.type(value()); - tasks[t] = true; - break; - - case task_system_force: - if (!tasks[task_extended_lagrangian]) { - if (!b_inverse_gradients) { - cvm::error("Error: one or more of the components of " - "colvar \""+this->name+ - "\" does not support system force calculation.\n", - INPUT_ERROR); - } - cvm::request_system_force(); - } - ft.type(value()); - ft_reported.type(value()); - tasks[t] = true; - break; - - case task_output_applied_force: - case task_lower_wall: - case task_upper_wall: - // all of the above require gradients - enable(task_gradients); - tasks[t] = true; - break; - - case task_fdiff_velocity: - x_old.type(value()); - v_fdiff.type(value()); - v_reported.type(value()); - tasks[t] = true; - break; - - case task_output_velocity: - enable(task_fdiff_velocity); - tasks[t] = true; - break; - - case task_grid: - if (value().type() != colvarvalue::type_scalar) { - cvm::error("Cannot calculate a grid for collective variable, \""+ - this->name+"\", because its value is not a scalar number.\n", - INPUT_ERROR); - } - tasks[t] = true; - break; - - case task_extended_lagrangian: - enable(task_gradients); - v_reported.type(value()); - tasks[t] = true; - break; - - case task_lower_boundary: - case task_upper_boundary: - if (value().type() != colvarvalue::type_scalar) { - cvm::error("Error: this colvar is not a scalar value " - "and cannot produce a grid.\n", - INPUT_ERROR); - } - tasks[t] = true; - break; - - case task_output_value: - case task_runave: - case task_corrfunc: - case task_langevin: - case task_output_energy: - case task_scripted: - case task_gradients: - tasks[t] = true; - break; - - case task_collect_gradients: - if (value().type() != colvarvalue::type_scalar) { - cvm::error("Collecting atomic gradients for non-scalar collective variable \""+ - this->name+"\" is not yet implemented.\n", - INPUT_ERROR); - } - - enable(task_gradients); - if (atom_ids.size() == 0) { - build_atom_list(); - } - tasks[t] = true; - break; - - default: - break; - } - - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); -} - - -void colvar::disable(colvar::task const &t) -{ - // check dependencies - switch (t) { - case task_gradients: - disable(task_upper_wall); - disable(task_lower_wall); - disable(task_output_applied_force); - disable(task_system_force); - disable(task_Jacobian_force); - tasks[t] = false; - break; - - case task_system_force: - disable(task_output_system_force); - tasks[t] = false; - break; - - case task_Jacobian_force: - disable(task_report_Jacobian_force); - tasks[t] = false; - break; - - case task_fdiff_velocity: - disable(task_output_velocity); - tasks[t] = false; - break; - - case task_lower_boundary: - case task_upper_boundary: - disable(task_grid); - tasks[t] = false; - break; - - case task_extended_lagrangian: - case task_report_Jacobian_force: - case task_output_value: - case task_output_velocity: - case task_output_applied_force: - case task_output_system_force: - case task_runave: - case task_corrfunc: - case task_grid: - case task_lower_wall: - case task_upper_wall: - case task_langevin: - case task_output_energy: - case task_collect_gradients: - case task_scripted: - tasks[t] = false; - break; - - default: - break; - } - -} - - void colvar::setup() { // loop over all components to reset masses of all groups for (size_t i = 0; i < cvcs.size(); i++) { @@ -842,8 +684,18 @@ void colvar::setup() { colvar::~colvar() { - for (size_t i = 0; i < cvcs.size(); i++) { - delete cvcs[i]; +// Clear references to this colvar's cvcs as children +// for dependency purposes + remove_all_children(); + + for (std::vector::reverse_iterator ci = cvcs.rbegin(); + ci != cvcs.rend(); + ++ci) { + // clear all children of this cvc (i.e. its atom groups) + // because the cvc base class destructor can't do it early enough + // and we don't want to have each cvc derived class do it separately + (*ci)->remove_all_children(); + delete *ci; } // remove reference to this colvar from the CVM @@ -862,91 +714,114 @@ colvar::~colvar() // ******************** CALC FUNCTIONS ******************** -void colvar::calc() +// Default schedule (everything is serialized) +int colvar::calc() { - if (cvm::debug()) - cvm::log("Calculating colvar \""+this->name+"\".\n"); - update_cvc_flags(); - calc_cvcs(); - calc_colvar_properties(); + // Note: if anything is added here, it should be added also in the SMP block of calc_colvars() + int error_code = COLVARS_OK; + if (is_enabled(f_cv_active)) { + error_code |= update_cvc_flags(); + error_code |= calc_cvcs(); + error_code |= collect_cvc_data(); + } + return COLVARS_OK; } int colvar::calc_cvcs(int first_cvc, size_t num_cvcs) { - size_t i, ig; - size_t cvc_count; - size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); - - // prepare atom groups for calculation + int error_code = COLVARS_OK; if (cvm::debug()) - cvm::log("Collecting data from atom groups.\n"); + cvm::log("Calculating colvar \""+this->name+"\", components "+ + cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n"); - if (first_cvc >= cvcs.size()) { + error_code |= check_cvc_range(first_cvc, num_cvcs); + if (error_code != COLVARS_OK) { + return error_code; + } + + error_code |= calc_cvc_values(first_cvc, num_cvcs); + error_code |= calc_cvc_gradients(first_cvc, num_cvcs); + error_code |= calc_cvc_sys_forces(first_cvc, num_cvcs); + error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs); + + if (cvm::debug()) + cvm::log("Done calculating colvar \""+this->name+"\".\n"); + + return error_code; +} + + +int colvar::collect_cvc_data() +{ + if (cvm::debug()) + cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n"); + + int error_code = COLVARS_OK; + + error_code |= collect_cvc_values(); + error_code |= collect_cvc_gradients(); + error_code |= collect_cvc_sys_forces(); + error_code |= collect_cvc_Jacobians(); + + error_code |= calc_colvar_properties(); + + if (cvm::debug()) + cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n"); + + return error_code; +} + + +int colvar::check_cvc_range(int first_cvc, size_t num_cvcs) +{ + if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) { cvm::error("Error: trying to address a component outside the " "range defined for colvar \""+name+"\".\n", BUG_ERROR); return BUG_ERROR; } + return COLVARS_OK; +} - for (i = first_cvc, cvc_count = 0; - (i < cvcs.size()) && (cvc_count < cvc_max_count); - i++) { - if (!cvcs[i]->b_enabled) continue; - cvc_count++; - for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { - cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); - atoms.reset_atoms_data(); - atoms.read_positions(); - atoms.calc_required_properties(); - // each atom group will take care of its own ref_pos_group, if defined - } - } -//// Don't try to get atom velocities, as no back-end currently implements it -// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) { -// for (i = 0; i < cvcs.size(); i++) { -// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { -// cvcs[i]->atom_groups[ig]->read_velocities(); -// } -// } -// } - - if (tasks[task_system_force]) { - for (i = first_cvc, cvc_count = 0; - (i < cvcs.size()) && (cvc_count < cvc_max_count); - i++) { - if (!cvcs[i]->b_enabled) continue; - cvc_count++; - for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { - cvcs[i]->atom_groups[ig]->read_system_forces(); - } - } - } +int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs) +{ + size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); + size_t i, cvc_count; // calculate the value of the colvar if (cvm::debug()) cvm::log("Calculating colvar components.\n"); - x.reset(); // First, calculate component values + cvm::increase_depth(); for (i = first_cvc, cvc_count = 0; (i < cvcs.size()) && (cvc_count < cvc_max_count); i++) { - if (!cvcs[i]->b_enabled) continue; + if (!cvcs[i]->is_enabled()) continue; cvc_count++; - cvm::increase_depth(); + (cvcs[i])->read_data(); (cvcs[i])->calc_value(); - cvm::decrease_depth(); if (cvm::debug()) cvm::log("Colvar component no. "+cvm::to_str(i+1)+ " within colvar \""+this->name+"\" has value "+ cvm::to_str((cvcs[i])->value(), cvm::cv_width, cvm::cv_prec)+".\n"); } + cvm::decrease_depth(); - // Then combine them appropriately, using either a scripted function or a polynomial - if (tasks[task_scripted]) { + return COLVARS_OK; +} + + +int colvar::collect_cvc_values() +{ + x.reset(); + size_t i; + + // combine them appropriately, using either a scripted function or a polynomial + if (is_enabled(f_cv_scripted)) { // cvcs combined by user script int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x); if (res == COLVARS_NOT_IMPLEMENTED) { @@ -959,23 +834,16 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs) } } else if (x.type() == colvarvalue::type_scalar) { // polynomial combination allowed - for (i = first_cvc, cvc_count = 0; - (i < cvcs.size()) && (cvc_count < cvc_max_count); - i++) { - if (!cvcs[i]->b_enabled) continue; - cvc_count++; + for (i = 0; i < cvcs.size(); i++) { + if (!cvcs[i]->is_enabled()) continue; x += (cvcs[i])->sup_coeff * ( ((cvcs[i])->sup_np != 1) ? std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : (cvcs[i])->value().real_value ); } } else { - // only linear combination allowed - for (i = first_cvc, cvc_count = 0; - (i < cvcs.size()) && (cvc_count < cvc_max_count); - i++) { - if (!cvcs[i]->b_enabled) continue; - cvc_count++; + for (i = 0; i < cvcs.size(); i++) { + if (!cvcs[i]->is_enabled()) continue; x += (cvcs[i])->sup_coeff * (cvcs[i])->value(); } } @@ -984,131 +852,209 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs) cvm::log("Colvar \""+this->name+"\" has value "+ cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n"); - if (tasks[task_gradients]) { + return COLVARS_OK; +} + + +int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs) +{ + size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); + size_t i, cvc_count; + + if (is_enabled(f_cv_gradient)) { if (cvm::debug()) cvm::log("Calculating gradients of colvar \""+this->name+"\".\n"); // calculate the gradients of each component + cvm::increase_depth(); for (i = first_cvc, cvc_count = 0; (i < cvcs.size()) && (cvc_count < cvc_max_count); i++) { - if (!cvcs[i]->b_enabled) continue; + if (!cvcs[i]->is_enabled()) continue; cvc_count++; - cvm::increase_depth(); (cvcs[i])->calc_gradients(); // if requested, propagate (via chain rule) the gradients above // to the atoms used to define the roto-translation - for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { if (cvcs[i]->atom_groups[ig]->b_fit_gradients) cvcs[i]->atom_groups[ig]->calc_fit_gradients(); } - cvm::decrease_depth(); } + cvm::decrease_depth(); if (cvm::debug()) cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n"); + } - if (tasks[task_collect_gradients]) { + return COLVARS_OK; +} - if (tasks[task_scripted]) { - cvm::error("Collecting atomic gradients is not implemented for " - "scripted colvars.", COLVARS_NOT_IMPLEMENTED); - return COLVARS_NOT_IMPLEMENTED; - } - // Collect the atomic gradients inside colvar object - for (unsigned int a = 0; a < atomic_gradients.size(); a++) { - atomic_gradients[a].reset(); - } - for (i = first_cvc, cvc_count = 0; - (i < cvcs.size()) && (cvc_count < cvc_max_count); - i++) { - if (!cvcs[i]->b_enabled) continue; - cvc_count++; - // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx - cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) * - std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); +int colvar::collect_cvc_gradients() +{ + size_t i; - for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { + if (is_enabled(f_cv_collect_gradient)) { - cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + if (is_enabled(f_cv_scripted)) { + cvm::error("Collecting atomic gradients is not implemented for " + "scripted colvars.", COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } - // If necessary, apply inverse rotation to get atomic - // gradient in the laboratory frame - if (cvcs[i]->atom_groups[j]->b_rotate) { - cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse(); + // Collect the atomic gradients inside colvar object + for (unsigned int a = 0; a < atomic_gradients.size(); a++) { + atomic_gradients[a].reset(); + } + for (i = 0; i < cvcs.size(); i++) { + if (!cvcs[i]->is_enabled()) continue; + // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx + cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) * + std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); - for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - ag[k].id) - atom_ids.begin(); - atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad); - } + for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { - } else { + cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); - for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - ag[k].id) - atom_ids.begin(); - atomic_gradients[a] += coeff * ag[k].grad; - } + // If necessary, apply inverse rotation to get atomic + // gradient in the laboratory frame + if (ag.b_rotate) { + cvm::rotation const rot_inv = ag.rot.inverse(); + + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad); + } + + } else { + + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; } } } } } + return COLVARS_OK; +} - if (tasks[task_system_force] && !tasks[task_extended_lagrangian]) { + +int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs) +{ + size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); + size_t i, cvc_count; + + if (is_enabled(f_cv_system_force) && !is_enabled(f_cv_extended_Lagrangian)) { // If extended Lagrangian is enabled, system force calculation is trivial // and done together with integration of the extended coordinate. - if (tasks[task_scripted]) { + if (is_enabled(f_cv_scripted)) { // TODO see if this could reasonably be done in a generic way // from generic inverse gradients cvm::error("System force is not implemented for " "scripted colvars.", COLVARS_NOT_IMPLEMENTED); return COLVARS_NOT_IMPLEMENTED; } + if (cvm::debug()) cvm::log("Calculating system force of colvar \""+this->name+"\".\n"); - ft.reset(); - // if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { - // Disabled check to allow for explicit system force calculation + // Disabled check to allow for explicit system force calculation // even with extended Lagrangian if (cvm::step_relative() > 0) { + cvm::increase_depth(); // get from the cvcs the system forces from the PREVIOUS step for (i = first_cvc, cvc_count = 0; (i < cvcs.size()) && (cvc_count < cvc_max_count); i++) { - if (!cvcs[i]->b_enabled) continue; + if (!cvcs[i]->is_enabled()) continue; cvc_count++; (cvcs[i])->calc_force_invgrads(); - // linear combination is assumed - cvm::increase_depth(); - ft += (cvcs[i])->system_force() / ((cvcs[i])->sup_coeff * cvm::real(cvcs.size())); - cvm::decrease_depth(); } - } - - if (tasks[task_report_Jacobian_force]) { - // add the Jacobian force to the system force, and don't apply any silent - // correction internally: biases such as colvarbias_abf will handle it - ft += fj; + cvm::decrease_depth(); } if (cvm::debug()) cvm::log("Done calculating system force of colvar \""+this->name+"\".\n"); } + + return COLVARS_OK; +} + + +int colvar::collect_cvc_sys_forces() +{ + if (is_enabled(f_cv_system_force) && !is_enabled(f_cv_extended_Lagrangian)) { + // If extended Lagrangian is enabled, system force calculation is trivial + // and done together with integration of the extended coordinate. + ft.reset(); + + if (cvm::step_relative() > 0) { + // get from the cvcs the system forces from the PREVIOUS step + for (size_t i = 0; i < cvcs.size(); i++) { + if (!cvcs[i]->is_enabled()) continue; + // linear combination is assumed + ft += (cvcs[i])->system_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm; + } + } + + if (!is_enabled(f_cv_hide_Jacobian)) { + // add the Jacobian force to the system force, and don't apply any silent + // correction internally: biases such as colvarbias_abf will handle it + ft += fj; + } + } + + return COLVARS_OK; +} + + +int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs) +{ + size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); + + if (is_enabled(f_cv_Jacobian)) { + cvm::increase_depth(); + size_t i, cvc_count; + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { + if (!cvcs[i]->is_enabled()) continue; + cvc_count++; + (cvcs[i])->calc_Jacobian_derivative(); + } + cvm::decrease_depth(); + } + + return COLVARS_OK; +} + + +int colvar::collect_cvc_Jacobians() +{ + if (is_enabled(f_cv_Jacobian)) { + fj.reset(); + for (size_t i = 0; i < cvcs.size(); i++) { + if (!cvcs[i]->is_enabled()) continue; + // linear combination is assumed + fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm; + } + fj *= cvm::boltzmann() * cvm::temperature(); + } + return COLVARS_OK; } int colvar::calc_colvar_properties() { - if (tasks[task_fdiff_velocity]) { + if (is_enabled(f_cv_fdiff_velocity)) { // calculate the velocity by finite differences if (cvm::step_relative() == 0) x_old = x; @@ -1118,7 +1064,7 @@ int colvar::calc_colvar_properties() } } - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { // initialize the restraint center in the first step to the value // just calculated from the cvcs @@ -1142,14 +1088,11 @@ int colvar::calc_colvar_properties() ft_reported = ft; } - if (cvm::debug()) - cvm::log("Done calculating colvar \""+this->name+"\".\n"); - return COLVARS_OK; } -cvm::real colvar::update() +cvm::real colvar::update_forces_energy() { if (cvm::debug()) cvm::log("Updating colvar \""+this->name+"\".\n"); @@ -1162,35 +1105,14 @@ cvm::real colvar::update() // been summed over each bias using this colvar f += fb; - - if (tasks[task_Jacobian_force]) { - size_t i; - for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; - cvm::increase_depth(); - (cvcs[i])->calc_Jacobian_derivative(); - cvm::decrease_depth(); - } - - size_t ncvc = 0; - fj.reset(); - for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; - // linear combination is assumed - fj += 1.0 / cvm::real((cvcs[i])->sup_coeff) * - (cvcs[i])->Jacobian_derivative(); - ncvc++; - } - fj *= (1.0/cvm::real(ncvc)) * cvm::boltzmann() * cvm::temperature(); - + if (is_enabled(f_cv_Jacobian)) { // the instantaneous Jacobian force was not included in the reported system force; // instead, it is subtracted from the applied force (silent Jacobian correction) - if (! tasks[task_report_Jacobian_force]) + if (is_enabled(f_cv_hide_Jacobian)) f -= fj; } - - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { cvm::real dt = cvm::dt(); cvm::real f_ext; @@ -1212,7 +1134,7 @@ cvm::real colvar::update() kinetic_energy = 0.5 * ext_mass * vr * vr; potential_energy = 0.5 * ext_force_k * this->dist2(xr, x); // leap to v_(i+1/2) - if (tasks[task_langevin]) { + if (is_enabled(f_cv_Langevin)) { vr -= dt * ext_gamma * vr.real_value; vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass; } @@ -1224,7 +1146,7 @@ cvm::real colvar::update() // Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use - if (tasks[task_lower_wall] || tasks[task_upper_wall]) { + if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) { // Wall force colvarvalue fw(x); @@ -1235,7 +1157,7 @@ cvm::real colvar::update() // For a periodic colvar, both walls may be applicable at the same time // in which case we pick the closer one - if ( (!tasks[task_upper_wall]) || + if ( (!is_enabled(f_cv_upper_wall)) || (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) { cvm::real const grad = this->dist2_lgrad(x, lower_wall); @@ -1261,7 +1183,7 @@ cvm::real colvar::update() } - if (tasks[task_fdiff_velocity]) { + if (is_enabled(f_cv_fdiff_velocity)) { // set it for the next step x_old = x; } @@ -1278,11 +1200,11 @@ void colvar::communicate_forces() if (cvm::debug()) cvm::log("Communicating forces from colvar \""+this->name+"\".\n"); - if (tasks[task_scripted]) { + if (is_enabled(f_cv_scripted)) { std::vector > func_grads; func_grads.reserve(cvcs.size()); for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; + if (!cvcs[i]->is_enabled()) continue; func_grads.push_back(cvm::matrix2d (x.size(), cvcs[i]->value().size())); } @@ -1299,33 +1221,27 @@ void colvar::communicate_forces() int grad_index = 0; // index in the scripted gradients, to account for some components being disabled for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; - cvm::increase_depth(); + if (!cvcs[i]->is_enabled()) continue; // cvc force is colvar force times colvar/cvc Jacobian // (vector-matrix product) (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++], cvcs[i]->value().type())); - cvm::decrease_depth(); } } else if (x.type() == colvarvalue::type_scalar) { for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; - cvm::increase_depth(); + if (!cvcs[i]->is_enabled()) continue; (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) * (std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1)) ); - cvm::decrease_depth(); } } else { for (i = 0; i < cvcs.size(); i++) { - if (!cvcs[i]->b_enabled) continue; - cvm::increase_depth(); + if (!cvcs[i]->is_enabled()) continue; (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff); - cvm::decrease_depth(); } } @@ -1349,45 +1265,38 @@ int colvar::set_cvc_flags(std::vector const &flags) int colvar::update_cvc_flags() { - size_t i; // Update the enabled/disabled status of cvcs if necessary if (cvc_flags.size()) { - bool any = false; - for (i = 0; i < cvcs.size(); i++) { - cvcs[i]->b_enabled = cvc_flags[i]; - any = any || cvc_flags[i]; + n_active_cvcs = 0; + active_cvc_square_norm = 0.; + + for (size_t i = 0; i < cvcs.size(); i++) { + cvcs[i]->feature_states[f_cvc_active]->enabled = cvc_flags[i]; + if (cvcs[i]->is_enabled()) { + n_active_cvcs++; + active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff; + } } - if (!any) { + if (!n_active_cvcs) { cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n"); return COLVARS_ERROR; } cvc_flags.resize(0); } + return COLVARS_OK; } -int colvar::num_active_cvcs() const -{ - int result = 0; - for (size_t i = 0; i < cvcs.size(); i++) { - if (cvcs[i]->b_enabled) result++; - } - return result; -} - - // ******************** METRIC FUNCTIONS ******************** // Use the metrics defined by \link cvc \endlink objects bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const { - if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { - cvm::log("Error: requesting to histogram the " - "collective variable \""+this->name+"\", but a " - "pair of lower and upper boundaries must be " - "defined.\n"); + if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) { + cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" " + "requires lower and upper boundaries to be defined.\n"); cvm::set_error_bits(INPUT_ERROR); } @@ -1403,11 +1312,9 @@ bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) c bool colvar::periodic_boundaries() const { - if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { - cvm::error("Error: requesting to histogram the " - "collective variable \""+this->name+"\", but a " - "pair of lower and upper boundaries must be " - "defined.\n"); + if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) { + cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" " + "requires lower and upper boundaries to be defined.\n"); } return periodic_boundaries(lower_boundary, upper_boundary); @@ -1417,7 +1324,7 @@ bool colvar::periodic_boundaries() const cvm::real colvar::dist2(colvarvalue const &x1, colvarvalue const &x2) const { - if (b_homogeneous) { + if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2(x1, x2); } else { return x1.dist2(x2); @@ -1427,7 +1334,7 @@ cvm::real colvar::dist2(colvarvalue const &x1, colvarvalue colvar::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - if (b_homogeneous) { + if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2_lgrad(x1, x2); } else { return x1.dist2_grad(x2); @@ -1437,7 +1344,7 @@ colvarvalue colvar::dist2_lgrad(colvarvalue const &x1, colvarvalue colvar::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { - if (b_homogeneous) { + if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2_rgrad(x1, x2); } else { return x2.dist2_grad(x1); @@ -1446,7 +1353,7 @@ colvarvalue colvar::dist2_rgrad(colvarvalue const &x1, void colvar::wrap(colvarvalue &x) const { - if (b_homogeneous) { + if (is_enabled(f_cv_homogeneous)) { (cvcs[0])->wrap(x); } return; @@ -1492,7 +1399,7 @@ std::istream & colvar::read_restart(std::istream &is) cvm::to_str(x)+"\n"); } - if (tasks[colvar::task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { if ( !(get_keyval(conf, "extended_x", xr, colvarvalue(x.type()), colvarparse::parse_silent)) && @@ -1502,15 +1409,12 @@ std::istream & colvar::read_restart(std::istream &is) "\"extended_x\" or \"extended_v\" for the colvar \""+ name+"\", but you requested \"extendedLagrangian\".\n"); } - } - - if (tasks[task_extended_lagrangian]) { x_reported = xr; } else { x_reported = x; } - if (tasks[task_output_velocity]) { + if (is_enabled(f_cv_output_velocity)) { if ( !(get_keyval(conf, "v", v_fdiff, colvarvalue(x.type()), colvarparse::parse_silent)) ) { @@ -1519,7 +1423,7 @@ std::istream & colvar::read_restart(std::istream &is) name+"\", but you requested \"outputVelocity\".\n"); } - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { v_reported = vr; } else { v_reported = v_fdiff; @@ -1534,7 +1438,7 @@ std::istream & colvar::read_traj(std::istream &is) { size_t const start_pos = is.tellg(); - if (tasks[task_output_value]) { + if (is_enabled(f_cv_output_value)) { if (!(is >> x)) { cvm::log("Error: in reading the value of colvar \""+ @@ -1545,7 +1449,7 @@ std::istream & colvar::read_traj(std::istream &is) return is; } - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { is >> xr; x_reported = xr; } else { @@ -1553,11 +1457,11 @@ std::istream & colvar::read_traj(std::istream &is) } } - if (tasks[task_output_velocity]) { + if (is_enabled(f_cv_output_velocity)) { is >> v_fdiff; - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { is >> vr; v_reported = vr; } else { @@ -1565,12 +1469,12 @@ std::istream & colvar::read_traj(std::istream &is) } } - if (tasks[task_output_system_force]) { + if (is_enabled(f_cv_output_system_force)) { is >> ft; ft_reported = ft; } - if (tasks[task_output_applied_force]) { + if (is_enabled(f_cv_output_applied_force)) { is >> f; } @@ -1589,14 +1493,14 @@ std::ostream & colvar::write_restart(std::ostream &os) { << std::setw(cvm::cv_width) << x << "\n"; - if (tasks[task_output_velocity]) { + if (is_enabled(f_cv_output_velocity)) { os << " v " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << v_reported << "\n"; } - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { os << " extended_x " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) @@ -1619,43 +1523,43 @@ std::ostream & colvar::write_traj_label(std::ostream & os) os << " "; - if (tasks[task_output_value]) { + if (is_enabled(f_cv_output_value)) { os << " " << cvm::wrap_string(this->name, this_cv_width); - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { // extended DOF os << " r_" << cvm::wrap_string(this->name, this_cv_width-2); } } - if (tasks[task_output_velocity]) { + if (is_enabled(f_cv_output_velocity)) { os << " v_" << cvm::wrap_string(this->name, this_cv_width-2); - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { // extended DOF os << " vr_" << cvm::wrap_string(this->name, this_cv_width-3); } } - if (tasks[task_output_energy]) { + if (is_enabled(f_cv_output_energy)) { os << " Ep_" << cvm::wrap_string(this->name, this_cv_width-3) << " Ek_" << cvm::wrap_string(this->name, this_cv_width-3); } - if (tasks[task_output_system_force]) { + if (is_enabled(f_cv_output_system_force)) { os << " fs_" << cvm::wrap_string(this->name, this_cv_width-3); } - if (tasks[task_output_applied_force]) { + if (is_enabled(f_cv_output_applied_force)) { os << " fa_" << cvm::wrap_string(this->name, this_cv_width-3); } @@ -1668,9 +1572,9 @@ std::ostream & colvar::write_traj(std::ostream &os) { os << " "; - if (tasks[task_output_value]) { + if (is_enabled(f_cv_output_value)) { - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << x; @@ -1681,9 +1585,9 @@ std::ostream & colvar::write_traj(std::ostream &os) << x_reported; } - if (tasks[task_output_velocity]) { + if (is_enabled(f_cv_output_velocity)) { - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_extended_Lagrangian)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << v_fdiff; @@ -1694,7 +1598,7 @@ std::ostream & colvar::write_traj(std::ostream &os) << v_reported; } - if (tasks[task_output_energy]) { + if (is_enabled(f_cv_output_energy)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << potential_energy @@ -1703,14 +1607,14 @@ std::ostream & colvar::write_traj(std::ostream &os) } - if (tasks[task_output_system_force]) { + if (is_enabled(f_cv_output_system_force)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << ft_reported; } - if (tasks[task_output_applied_force]) { - if (tasks[task_extended_lagrangian]) { + if (is_enabled(f_cv_output_applied_force)) { + if (is_enabled(f_cv_extended_Lagrangian)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << fr; @@ -1731,15 +1635,16 @@ int colvar::write_output_files() if (acf.size()) { cvm::log("Writing acf to file \""+acf_outfile+"\".\n"); + cvm::backup_file(acf_outfile.c_str()); cvm::ofstream acf_os(acf_outfile.c_str()); - if (! acf_os.good()) { + if (! acf_os.is_open()) { cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR); } write_acf(acf_os); acf_os.close(); } - if (runave_os.good()) { + if (runave_os.is_open()) { runave_os.close(); } } @@ -1752,11 +1657,11 @@ int colvar::write_output_files() void colvar::analyze() { - if (tasks[task_runave]) { + if (is_enabled(f_cv_runave)) { calc_runave(); } - if (tasks[task_corrfunc]) { + if (is_enabled(f_cv_corrfunc)) { calc_acf(); } } @@ -1841,7 +1746,7 @@ int colvar::calc_acf() case acf_vel: - if (tasks[task_fdiff_velocity]) { + if (is_enabled(f_cv_fdiff_velocity)) { // calc() should do this already, but this only happens in a // simulation; better do it again in case a trajectory is // being read @@ -1874,7 +1779,7 @@ int colvar::calc_acf() } } - if (tasks[task_fdiff_velocity]) { + if (is_enabled(f_cv_fdiff_velocity)) { // set it for the next step x_old = x; } @@ -2032,4 +1937,6 @@ void colvar::calc_runave() } +// Static members +std::vector colvar::cv_features; diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 7ab4467bfb..ea36d3e443 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -10,6 +10,7 @@ #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" +#include "colvardeps.h" /// \brief A collective variable (main class); to be defined, it needs @@ -37,7 +38,7 @@ /// \link colvarvalue \endlink type, you should also add its /// initialization line in the \link colvar \endlink constructor. -class colvar : public colvarparse { +class colvar : public colvarparse, public cvm::deps { public: @@ -63,8 +64,6 @@ public: /// subtracted. colvarvalue const & system_force() const; - /// \brief - /// \brief Typical fluctuation amplitude for this collective /// variable (e.g. local width of a free energy basin) /// @@ -75,87 +74,15 @@ public: /// variable. cvm::real width; - /// \brief True if this \link colvar \endlink is a linear - /// combination of \link cvc \endlink elements - bool b_linear; + /// \brief Implementation of the feature list for colvar + static std::vector cv_features; - /// \brief True if this \link colvar \endlink is a linear - /// combination of cvcs with coefficients 1 or -1 - bool b_homogeneous; + /// \brief Implementation of the feature list accessor for colvar + std::vector &features() { + return cv_features; + } - /// \brief True if all \link cvc \endlink objects are capable - /// of calculating inverse gradients - bool b_inverse_gradients; - - /// \brief True if all \link cvc \endlink objects are capable - /// of calculating Jacobian forces - bool b_Jacobian_force; - - /// \brief Options controlling the behaviour of the colvar during - /// the simulation, which are set from outside the cvcs - enum task { - /// \brief Gradients are calculated and temporarily stored, so - /// that external forces can be applied - task_gradients, - /// \brief Collect atomic gradient data from all cvcs into vector - /// atomic_gradients - task_collect_gradients, - /// \brief Calculate the velocity with finite differences - task_fdiff_velocity, - /// \brief The system force is calculated, projecting the atomic - /// forces on the inverse gradients - task_system_force, - /// \brief The variable has a harmonic restraint around a moving - /// center with fictitious mass; bias forces will be applied to - /// the center - task_extended_lagrangian, - /// \brief The extended system coordinate undergoes Langevin - /// dynamics - task_langevin, - /// \brief Output the potential and kinetic energies - /// (for extended Lagrangian colvars only) - task_output_energy, - /// \brief Compute analytically the "force" arising from the - /// geometric entropy component (for example, that from the angular - /// states orthogonal to a distance vector) - task_Jacobian_force, - /// \brief Report the Jacobian force as part of the system force - /// if disabled, apply a correction internally to cancel it - task_report_Jacobian_force, - /// \brief Output the value to the trajectory file (on by default) - task_output_value, - /// \brief Output the velocity to the trajectory file - task_output_velocity, - /// \brief Output the applied force to the trajectory file - task_output_applied_force, - /// \brief Output the system force to the trajectory file - task_output_system_force, - /// \brief A lower boundary is defined - task_lower_boundary, - /// \brief An upper boundary is defined - task_upper_boundary, - /// \brief Provide a discretization of the values of the colvar to - /// be used by the biases or in analysis (needs lower and upper - /// boundary) - task_grid, - /// \brief Apply a restraining potential (|x-xb|^2) to the colvar - /// when it goes below the lower wall - task_lower_wall, - /// \brief Apply a restraining potential (|x-xb|^2) to the colvar - /// when it goes above the upper wall - task_upper_wall, - /// \brief Compute running average - task_runave, - /// \brief Compute time correlation function - task_corrfunc, - /// \brief Value and gradients computed by user script - task_scripted, - /// \brief Number of possible tasks - task_ntot - }; - - /// Tasks performed by this colvar - bool tasks[task_ntot]; + int refresh_deps(); /// List of biases that depend on this colvar std::vector biases; @@ -213,7 +140,7 @@ protected: /// Cached reported velocity colvarvalue v_reported; - // Options for task_extended_lagrangian + // Options for extended_lagrangian /// Restraint center colvarvalue xr; /// Velocity of the restraint center @@ -230,7 +157,7 @@ protected: /// \brief Harmonic restraint force colvarvalue fr; - /// \brief Jacobian force, when task_Jacobian_force is enabled + /// \brief Jacobian force, when Jacobian_force is enabled colvarvalue fj; /// Cached reported system force @@ -243,7 +170,7 @@ public: /// the biases are updated colvarvalue fb; - /// \brief Total \em applied force; fr (if task_extended_lagrangian + /// \brief Total \em applied force; fr (if extended_lagrangian /// is defined), fb (if biases are applied) and the walls' forces /// (if defined) contribute to it colvarvalue f; @@ -293,12 +220,6 @@ public: /// Constructor colvar(std::string const &conf); - /// Enable the specified task - int enable(colvar::task const &t); - - /// Disable the specified task - void disable(colvar::task const &t); - /// Get ready for a run and possibly re-initialize internal data void setup(); @@ -306,16 +227,46 @@ public: ~colvar(); - /// \brief Calculate the colvar value and all the other requested - /// quantities - void calc(); + /// \brief Calculate the colvar's value and related quantities + int calc(); - /// \brief Calculate a given subset of colvar components (CVCs) (default: all CVCs) + /// \brief Calculate a subset of the colvar components (CVCs) currently active + /// (default: all active CVCs) + /// Note: both arguments refer to the sect of *active* CVCs, not all CVCs int calc_cvcs(int first = 0, size_t num_cvcs = 0); + /// Ensure that the selected range of CVCs is consistent + int check_cvc_range(int first_cvc, size_t num_cvcs); + + /// \brief Calculate the values of the given subset of CVCs + int calc_cvc_values(int first, size_t num_cvcs); + /// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients + int calc_cvc_gradients(int first, size_t num_cvcs); + /// \brief Same as \link colvar::calc_cvc_values \endlink but for system forces + int calc_cvc_sys_forces(int first, size_t num_cvcs); + /// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces + int calc_cvc_Jacobians(int first, size_t num_cvcs); + + /// \brief Collect quantities from CVCs and update aggregated data for the colvar + int collect_cvc_data(); + + /// \brief Collect the values of the CVCs + int collect_cvc_values(); + /// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients + int collect_cvc_gradients(); + /// \brief Same as \link colvar::collect_cvc_values \endlink but for system forces + int collect_cvc_sys_forces(); + /// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces + int collect_cvc_Jacobians(); /// \brief Calculate the quantities associated to the colvar (but not to the CVCs) int calc_colvar_properties(); + /// Get the current biasing force + inline colvarvalue bias_force() const + { + return fb; + } + /// Set the total biasing force to zero void reset_bias_force(); @@ -326,7 +277,7 @@ public: /// equations of motion of internal degrees of freedom; see also /// colvar::communicate_forces() /// return colvar energy if extended Lagrandian active - cvm::real update(); + cvm::real update_forces_energy(); /// \brief Communicate forces (previously calculated in /// colvar::update()) to the external degrees of freedom @@ -335,11 +286,19 @@ public: /// \brief Enables and disables individual CVCs based on flags int set_cvc_flags(std::vector const &flags); - /// \brief Updates the flags in the CVC objects + /// \brief Updates the flags in the CVC objects, and their number int update_cvc_flags(); - /// \brief Return the number of CVC objects with an active flag - int num_active_cvcs() const; +protected: + /// \brief Number of CVC objects with an active flag + size_t n_active_cvcs; + + /// Sum of square coefficients for active cvcs + cvm::real active_cvc_square_norm; + +public: + /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags) + inline size_t num_active_cvcs() const { return n_active_cvcs; } /// \brief Use the internal metrics (as from \link cvc /// \endlink objects) to calculate square distances and gradients @@ -528,7 +487,7 @@ protected: std::vector cvc_flags; /// \brief Initialize the sorted list of atom IDs for atoms involved - /// in all cvcs (called when enabling task_collect_gradients) + /// in all cvcs (called when enabling f_cv_collect_gradients) void build_atom_list(void); private: @@ -580,6 +539,9 @@ inline colvarvalue const & colvar::system_force() const inline void colvar::add_bias_force(colvarvalue const &force) { + if (cvm::debug()) { + cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n"); + } fb += force; } diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 9418b52e12..f13f6993b6 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -4,7 +4,6 @@ #include "colvarparse.h" #include "colvaratoms.h" - cvm::atom::atom() { index = -1; @@ -68,8 +67,7 @@ cvm::atom_group::atom_group(std::string const &conf, char const *key_in) { key = key_in; - cvm::log("Defining atom group \""+ - std::string(key)+"\".\n"); + cvm::log("Defining atom group \"" + key + "\".\n"); init(); // real work is done by parse parse(conf); @@ -93,8 +91,9 @@ cvm::atom_group::atom_group() cvm::atom_group::~atom_group() { - if (index >= 0) { + if (is_enabled(f_ag_scalable)) { (cvm::proxy)->clear_atom_group(index); + index = -1; } if (ref_pos_group) { @@ -151,7 +150,7 @@ int cvm::atom_group::add_atom_id(int aid) int cvm::atom_group::remove_atom(cvm::atom_iter ai) { - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR); return COLVARS_ERROR; } @@ -173,10 +172,11 @@ int cvm::atom_group::remove_atom(cvm::atom_iter ai) int cvm::atom_group::init() { if (!key.size()) key = "atoms"; - atoms.clear(); - b_scalable = false; + // TODO: check with proxy whether atom forces etc are available + init_ag_requires(); + index = -1; b_center = false; @@ -216,7 +216,7 @@ void cvm::atom_group::update_total_mass() return; } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { total_mass = (cvm::proxy)->get_atom_group_mass(index); } else { total_mass = 0.0; @@ -243,7 +243,7 @@ void cvm::atom_group::update_total_charge() return; } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { total_charge = (cvm::proxy)->get_atom_group_charge(index); } else { total_charge = 0.0; @@ -280,6 +280,8 @@ int cvm::atom_group::parse(std::string const &conf) cvm::log("Initializing atom group \""+key+"\".\n"); + description = "atom group " + key; + // whether or not to include messages in the log // colvarparse::Parse_Mode mode = parse_silent; // { @@ -291,9 +293,6 @@ int cvm::atom_group::parse(std::string const &conf) int parse_error = COLVARS_OK; - // if the cvc allows it, the flag has been set to true by default - get_keyval(group_conf, "scalable", b_scalable, b_scalable); - { std::string numbers_conf = ""; size_t pos = 0; @@ -373,11 +372,7 @@ int cvm::atom_group::parse(std::string const &conf) } // Catch any errors from all the initialization steps above - if (parse_error || cvm::get_error()) return COLVARS_ERROR; - - if (b_scalable) { - index = (cvm::proxy)->init_atom_group(atoms_ids); - } + if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error()); // checks of doubly-counted atoms have been handled by add_atom() already @@ -408,6 +403,10 @@ int cvm::atom_group::parse(std::string const &conf) } } + if (is_enabled(f_ag_scalable) && !b_dummy) { + index = (cvm::proxy)->init_atom_group(atoms_ids); + } + parse_error |= parse_fitting_options(group_conf); // TODO move this to colvarparse object @@ -449,7 +448,7 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf) if (atom_indexes.size()) { atoms_ids.reserve(atoms_ids.size()+atom_indexes.size()); - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { for (size_t i = 0; i < atom_indexes.size(); i++) { add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i])); } @@ -490,7 +489,7 @@ int cvm::atom_group::add_index_group(std::string const &index_group_name) atoms_ids.reserve(atoms_ids.size()+index_groups_i->size()); - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { for (size_t i = 0; i < index_groups_i->size(); i++) { add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i])); } @@ -520,7 +519,7 @@ int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf) atoms_ids.reserve(atoms_ids.size() + (final - initial + 1)); - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { for (int anum = initial; anum <= final; anum++) { add_atom_id((cvm::proxy)->check_atom_id(anum)); } @@ -558,7 +557,7 @@ int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid, atoms_ids.reserve(atoms_ids.size() + (final - initial + 1)); - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { for (int resid = initial; resid <= final; resid++) { add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid)); } @@ -750,15 +749,24 @@ void cvm::atom_group::read_positions() int cvm::atom_group::calc_required_properties() { - if (b_dummy) return COLVARS_OK; - // TODO check if the com is needed? calc_center_of_mass(); - if (!b_scalable) { - // TODO check if calc_center_of_geometry() is needed without a fit? - calc_center_of_geometry(); + calc_center_of_geometry(); + + if (!is_enabled(f_ag_scalable)) { if (b_center || b_rotate) { + if (ref_pos_group) { + ref_pos_group->calc_center_of_geometry(); + } + calc_apply_roto_translation(); + + // update COM and COG after fitting + calc_center_of_geometry(); + calc_center_of_mass(); + if (ref_pos_group) { + ref_pos_group->calc_center_of_geometry(); + } } } @@ -767,37 +775,49 @@ int cvm::atom_group::calc_required_properties() return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } + void cvm::atom_group::calc_apply_roto_translation() { - atom_group *fit_group = ref_pos_group ? ref_pos_group : this; - if (b_center) { // center on the origin first - cvm::atom_pos const cog = fit_group->center_of_geometry(); + cvm::atom_pos const cog = ref_pos_group ? + ref_pos_group->center_of_geometry() : this->center_of_geometry(); apply_translation(-1.0 * cog); + if (ref_pos_group) { + ref_pos_group->apply_translation(-1.0 * cog); + } } if (b_rotate) { // rotate the group (around the center of geometry if b_center is // true, around the origin otherwise) - rot.calc_optimal_rotation(fit_group->positions(), ref_pos); + rot.calc_optimal_rotation(ref_pos_group ? + ref_pos_group->positions() : + this->positions(), + ref_pos); - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + cvm::atom_iter ai; + for (ai = this->begin(); ai != this->end(); ai++) { ai->pos = rot.rotate(ai->pos); } + if (ref_pos_group) { + for (ai = ref_pos_group->begin(); ai != ref_pos_group->end(); ai++) { + ai->pos = rot.rotate(ai->pos); + } + } } if (b_center) { // align with the center of geometry of ref_pos apply_translation(ref_pos_cog); - // update the center of geometry for external use - cog = ref_pos_cog; + if (ref_pos_group) { + ref_pos_group->apply_translation(ref_pos_cog); + } } - - // recalculate the center of mass - calc_center_of_mass(); + // update of COM and COG is done from the calling routine } + void cvm::atom_group::apply_translation(cvm::rvector const &t) { if (b_dummy) { @@ -805,7 +825,7 @@ void cvm::atom_group::apply_translation(cvm::rvector const &t) return; } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR); return; } @@ -815,23 +835,6 @@ void cvm::atom_group::apply_translation(cvm::rvector const &t) } } -void cvm::atom_group::apply_rotation(cvm::rotation const &rot) -{ - if (b_dummy) { - cvm::error("Error: cannot rotate the coordinates of a dummy atom group.\n", INPUT_ERROR); - return; - } - - if (b_scalable) { - cvm::error("Error: cannot rotate the coordinates of a scalable atom group.\n", INPUT_ERROR); - return; - } - - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { - ai->pos = rot.rotate(ai->pos - center_of_geometry()) + center_of_geometry(); - } -} - void cvm::atom_group::read_velocities() { @@ -852,6 +855,8 @@ void cvm::atom_group::read_velocities() } } + +// TODO make this a calc function void cvm::atom_group::read_system_forces() { if (b_dummy) return; @@ -891,7 +896,10 @@ int cvm::atom_group::calc_center_of_mass() { if (b_dummy) { com = dummy_atom_pos; - } else if (b_scalable) { + if (cvm::debug()) { + cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n"); + } + } else if (is_enabled(f_ag_scalable)) { com = (cvm::proxy)->get_atom_group_com(index); } else { com.reset(); @@ -922,7 +930,7 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad) { if (b_dummy) return; - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { scalar_com_gradient = grad; return; } @@ -993,7 +1001,7 @@ std::vector cvm::atom_group::positions() const "from a dummy atom group.\n", INPUT_ERROR); } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: atomic positions are not available " "from a scalable atom group.\n", INPUT_ERROR); } @@ -1014,7 +1022,7 @@ std::vector cvm::atom_group::positions_shifted(cvm::rvector const "from a dummy atom group.\n", INPUT_ERROR); } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: atomic positions are not available " "from a scalable atom group.\n", INPUT_ERROR); } @@ -1035,7 +1043,7 @@ std::vector cvm::atom_group::velocities() const "from a dummy atom group.\n", INPUT_ERROR); } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: atomic velocities are not available " "from a scalable atom group.\n", INPUT_ERROR); } @@ -1056,7 +1064,7 @@ std::vector cvm::atom_group::system_forces() const "from a dummy atom group.\n", INPUT_ERROR); } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { cvm::error("Error: atomic system forces are not available " "from a scalable atom group.\n", INPUT_ERROR); } @@ -1070,6 +1078,8 @@ std::vector cvm::atom_group::system_forces() const return f; } + +// TODO make this an accessor cvm::rvector cvm::atom_group::system_force() const { if (b_dummy) { @@ -1077,7 +1087,7 @@ cvm::rvector cvm::atom_group::system_force() const "from a dummy atom group.\n", INPUT_ERROR); } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { return (cvm::proxy)->get_atom_group_system_force(index); } @@ -1104,7 +1114,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force) return; } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { (cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient); return; } @@ -1156,7 +1166,7 @@ void cvm::atom_group::apply_force(cvm::rvector const &force) return; } - if (b_scalable) { + if (is_enabled(f_ag_scalable)) { (cvm::proxy)->apply_atom_group_force(index, force); } @@ -1209,3 +1219,8 @@ void cvm::atom_group::apply_forces(std::vector const &forces) } } +// Static members + +std::vector cvm::atom_group::ag_features; + + diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 746d4c2918..827601fc4e 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -5,6 +5,7 @@ #include "colvarmodule.h" #include "colvarparse.h" +#include "colvardeps.h" /// \brief Stores numeric id, mass and all mutable data for an atom, @@ -138,7 +139,7 @@ public: /// \brief Group of \link atom \endlink objects, mostly used by a /// \link cvc \endlink object to gather all atomic data class colvarmodule::atom_group - : public colvarparse + : public colvarparse, public cvm::deps { public: @@ -187,8 +188,13 @@ public: /// change atom masses after their initialization. void reset_mass(std::string &name, int i, int j); - /// \brief Whether or not the properties of this group will be computed in parallel - bool b_scalable; + /// \brief Implementation of the feature list for atom group + static std::vector ag_features; + + /// \brief Implementation of the feature list accessor for atom group + virtual std::vector &features() { + return ag_features; + } /// \brief Default constructor atom_group(); @@ -325,10 +331,6 @@ public: /// \brief Move all positions void apply_translation(cvm::rvector const &t); - /// \brief Rotate all positions around the center of geometry - void apply_rotation(cvm::rotation const &q); - - /// \brief Get the current velocities; this must be called always /// *after* read_positions(); if b_rotate is defined, the same /// rotation applied to the coordinates will be used diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index cb2037f0e5..94d74e47d4 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -10,6 +10,8 @@ colvarbias::colvarbias(std::string const &conf, char const *key) { cvm::log("Initializing a new \""+std::string(key)+"\" instance.\n"); + init_cvb_requires(); + size_t rank = 1; std::string const key_str(key); @@ -35,6 +37,8 @@ colvarbias::colvarbias(std::string const &conf, char const *key) return; } + description = "bias " + name; + // lookup the associated colvars std::vector colvars_str; if (get_keyval(conf, "colvars", colvars_str)) { @@ -46,6 +50,13 @@ colvarbias::colvarbias(std::string const &conf, char const *key) cvm::error("Error: no collective variables specified.\n"); return; } + for (size_t i=0; ienable(colvar::task_gradients); + // Removed this as nor all biases apply forces eg histogram + // cv->enable(colvar::task_gradients); if (cvm::debug()) cvm::log("Applying this bias to collective variable \""+ cv->name+"\".\n"); @@ -100,10 +112,11 @@ void colvarbias::add_colvar(std::string const &cv_name) } -cvm::real colvarbias::update() +int colvarbias::update() { + // Note: if anything is added here, it should be added also in the SMP block of calc_biases() has_data = true; - return 0.0; + return COLVARS_OK; } @@ -172,3 +185,7 @@ std::ostream & colvarbias::write_traj(std::ostream &os) << bias_energy; return os; } + +// Static members + +std::vector colvarbias::cvb_features; diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 3fe2b97833..58688a47bf 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -5,10 +5,11 @@ #include "colvar.h" #include "colvarparse.h" +#include "colvardeps.h" /// \brief Collective variable bias, base class -class colvarbias : public colvarparse { +class colvarbias : public colvarparse, public cvm::deps { public: /// Name of this bias @@ -19,7 +20,7 @@ public: /// Retrieve colvar values and calculate their biasing forces /// Return bias energy - virtual cvm::real update(); + virtual int update(); // TODO: move update_bias here (share with metadynamics) @@ -77,6 +78,14 @@ public: inline cvm::real get_energy() { return bias_energy; } + + /// \brief Implementation of the feature list for colvarbias + static std::vector cvb_features; + + /// \brief Implementation of the feature list accessor for colvarbias + virtual std::vector &features() { + return cvb_features; + } protected: /// \brief Pointers to collective variables to which the bias is diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 2d9202f4b3..cc549578ce 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -58,37 +58,22 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) cvm::error("Error: no collective variables specified for the ABF bias.\n"); } + if (update_bias) { + // Request calculation of system force (which also checks for availability) + enable(f_cvb_get_system_force); + } + if (apply_bias) { + enable(f_cvb_apply_force); + } + for (size_t i = 0; i < colvars.size(); i++) { if (colvars[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Error: ABF bias can only use scalar-type variables.\n"); } - - colvars[i]->enable(colvar::task_gradients); - - if (update_bias) { - // Request calculation of system force (which also checks for availability) - colvars[i]->enable(colvar::task_system_force); - - if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) { - // request computation of Jacobian force - // ultimately, will be done regardless of extended Lagrangian - // and colvar should then just report zero Jacobian force - colvars[i]->enable(colvar::task_Jacobian_force); - - // request Jacobian force as part as system force - // except if the user explicitly requires the "silent" Jacobian - // correction AND the colvar has a single component - if (hide_Jacobian) { - if (colvars[i]->n_components() > 1) { - cvm::log("WARNING: colvar \"" + colvars[i]->name - + "\" has multiple components; reporting its Jacobian forces\n"); - colvars[i]->enable(colvar::task_report_Jacobian_force); - } - } else { - colvars[i]->enable(colvar::task_report_Jacobian_force); - } - } + colvars[i]->enable(f_cv_grid); + if (hide_Jacobian) { + colvars[i]->enable(f_cv_hide_Jacobian); } // Here we could check for orthogonality of the Cartesian coordinates @@ -176,7 +161,7 @@ colvarbias_abf::~colvarbias_abf() /// Update the FE gradient, compute and apply biasing force /// also output data to disk if needed -cvm::real colvarbias_abf::update() +int colvarbias_abf::update() { if (cvm::debug()) cvm::log("Updating ABF bias " + this->name); @@ -283,7 +268,7 @@ cvm::real colvarbias_abf::update() cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+"."); } - return 0.0; + return COLVARS_OK; } int colvarbias_abf::replica_share() { diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index eb02d1f55e..1078ba866b 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -22,7 +22,7 @@ public: colvarbias_abf(std::string const &conf, char const *key); ~colvarbias_abf(); - cvm::real update(); + int update(); private: diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 78d7341af9..419d27cbd9 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -41,6 +41,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : current_coupling.resize(colvars.size()); coupling_rate.resize(colvars.size()); + enable(f_cvb_apply_force); for (i = 0; i < colvars.size(); i++) { colvar_centers[i].type(colvars[i]->value()); @@ -120,7 +121,7 @@ colvarbias_alb::~colvarbias_alb() { } -cvm::real colvarbias_alb::update() { +int colvarbias_alb::update() { bias_energy = 0.0; update_calls++; @@ -224,8 +225,7 @@ cvm::real colvarbias_alb::update() { } - return bias_energy; - + return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h index d8c46b499f..98de527f3a 100644 --- a/lib/colvars/colvarbias_alb.h +++ b/lib/colvars/colvarbias_alb.h @@ -13,7 +13,7 @@ public: virtual ~colvarbias_alb(); - virtual cvm::real update(); + virtual int update(); /// Read the bias configuration from a restart file virtual std::istream & read_restart(std::istream &is); diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index 21d68eeee4..bb55669614 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -10,6 +10,8 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const * : colvarbias(conf, key), grid(NULL), out_name("") { + size_t i; + get_keyval(conf, "outputFile", out_name, std::string("")); get_keyval(conf, "outputFileDX", out_name_dx, std::string("")); get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq); @@ -21,7 +23,6 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const * colvar_array_size = 0; { - size_t i; bool colvar_array = false; get_keyval(conf, "gatherVectorColvars", colvar_array, colvar_array); @@ -59,6 +60,10 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const * get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent); } + for (i = 0; i < colvars.size(); i++) { + colvars[i]->enable(f_cv_grid); + } + grid = new colvar_grid_scalar(); { @@ -86,10 +91,11 @@ colvarbias_histogram::~colvarbias_histogram() } /// Update the grid -cvm::real colvarbias_histogram::update() +int colvarbias_histogram::update() { + int error_code = COLVARS_OK; // update base class - colvarbias::update(); + error_code |= colvarbias::update(); if (cvm::debug()) { cvm::log("Updating histogram bias " + this->name); @@ -142,7 +148,7 @@ cvm::real colvarbias_histogram::update() write_output_files(); } - return 0.0; // no bias energy for histogram + return error_code | cvm::get_error(); } diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h index 90fb5441c8..60d1fdfcb2 100644 --- a/lib/colvars/colvarbias_histogram.h +++ b/lib/colvars/colvarbias_histogram.h @@ -19,7 +19,7 @@ public: colvarbias_histogram(std::string const &conf, char const *key); ~colvarbias_histogram(); - cvm::real update(); + int update(); int write_output_files(); diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 1fa804d303..46c0884e11 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -59,6 +59,9 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) comm = single_replica; } + // This implies gradients for all colvars + enable(f_cvb_apply_force); + get_keyval(conf, "useGrids", use_grids, true); if (use_grids) { @@ -68,6 +71,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) expand_grids = false; size_t i; for (i = 0; i < colvars.size(); i++) { + colvars[i]->enable(f_cv_grid); if (colvars[i]->expand_boundaries) { expand_grids = true; cvm::log("Metadynamics bias \""+this->name+"\""+ @@ -82,10 +86,6 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent); get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false); - for (i = 0; i < colvars.size(); i++) { - colvars[i]->enable(colvar::task_grid); - } - hills_energy = new colvar_grid_scalar(colvars); hills_energy_gradients = new colvar_grid_gradient(colvars); } else { @@ -253,7 +253,7 @@ colvarbias_meta::delete_hill(hill_iter &h) } -cvm::real colvarbias_meta::update() +int colvarbias_meta::update() { if (cvm::debug()) cvm::log("Updating the metadynamics bias \""+this->name+"\""+ @@ -544,7 +544,7 @@ cvm::real colvarbias_meta::update() ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); } - return bias_energy; + return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index d4c521ae62..2de1e2bade 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -36,7 +36,7 @@ public: /// Destructor virtual ~colvarbias_meta(); - virtual cvm::real update(); + virtual int update(); virtual std::istream & read_restart(std::istream &is); @@ -278,7 +278,7 @@ public: W(W_in), centers(cv.size()), widths(cv.size()), - it(cvm::it), + it(cvm::step_absolute()), replica(replica_in) { for (size_t i = 0; i < cv.size(); i++) { diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 5ac7b1318f..6812c670f7 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -18,6 +18,9 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf, colvar_centers.resize(colvars.size()); colvar_centers_raw.resize(colvars.size()); size_t i; + + enable(f_cvb_apply_force); + for (i = 0; i < colvars.size(); i++) { colvar_centers[i].type(colvars[i]->value()); colvar_centers_raw[i].type(colvars[i]->value()); @@ -53,6 +56,10 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf, size_t i; if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { + if (colvar_centers.size() != colvars.size()) { + cvm::error("Error: number of target centers does not match " + "that of collective variables.\n"); + } b_chg_centers = true; for (i = 0; i < target_centers.size(); i++) { target_centers[i].apply_constraints(); @@ -167,7 +174,7 @@ cvm::real colvarbias_restraint::energy_difference(std::string const &conf) } -cvm::real colvarbias_restraint::update() +int colvarbias_restraint::update() { bias_energy = 0.0; @@ -333,7 +340,7 @@ cvm::real colvarbias_restraint::update() cvm::log("Current forces for the restraint bias \""+ this->name+"\": "+cvm::to_str(colvar_forces)+".\n"); - return bias_energy; + return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index 3148fe0138..fb751bc10e 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -12,7 +12,7 @@ class colvarbias_restraint : public colvarbias { public: /// Retrieve colvar values and calculate their biasing forces - virtual cvm::real update(); + virtual int update(); /// Load new configuration - force constant and/or centers only virtual void change_configuration(std::string const &conf); diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 88e44ddd86..2f9ea78398 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -10,26 +10,22 @@ colvar::cvc::cvc() : sup_coeff(1.0), sup_np(1), - b_enabled(true), b_periodic(false), - b_inverse_gradients(false), - b_Jacobian_derivative(false), - b_debug_gradients(false) + b_try_scalable(true) {} colvar::cvc::cvc(std::string const &conf) : sup_coeff(1.0), sup_np(1), - b_enabled(true), b_periodic(false), - b_inverse_gradients(false), - b_Jacobian_derivative(false), - b_debug_gradients(false) + b_try_scalable(true) { if (cvm::debug()) cvm::log("Initializing cvc base object.\n"); + init_cvc_requires(); + get_keyval(conf, "name", this->name, std::string(""), parse_silent); get_keyval(conf, "componentCoeff", sup_coeff, 1.0); @@ -38,38 +34,109 @@ colvar::cvc::cvc(std::string const &conf) get_keyval(conf, "period", period, 0.0); get_keyval(conf, "wrapAround", wrap_center, 0.0); - get_keyval(conf, "debugGradients", b_debug_gradients, false, parse_silent); + { + bool b_debug_gradient; + get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent); + if (b_debug_gradient) enable(f_cvc_debug_gradient); + } + + // Attempt scalable calculations when in parallel? (By default yes, if available) + get_keyval(conf, "scalable", b_try_scalable, true); if (cvm::debug()) cvm::log("Done initializing cvc base object.\n"); } -void colvar::cvc::parse_group(std::string const &conf, - char const *group_key, - cvm::atom_group &group, - bool optional) +cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, + char const *group_key, + bool optional) { + cvm::atom_group *group = NULL; + if (key_lookup(conf, group_key)) { - // TODO turn on scalable flag for group objects in cvc init function - group.key = group_key; - if (group.parse(conf) != COLVARS_OK) { + group = new cvm::atom_group; + group->key = group_key; + + if (b_try_scalable) { + if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) { + enable(f_cvc_scalable_com); + enable(f_cvc_scalable); + group->enable(f_ag_scalable_com); + group->enable(f_ag_scalable); + } + + // TODO check for other types of parallelism here + + if (is_enabled(f_cvc_scalable)) { + cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n"); + } else { + cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n"); + } + } + + if (group->parse(conf) == COLVARS_OK) { + atom_groups.push_back(group); + } else { cvm::error("Error parsing definition for atom group \""+ std::string(group_key)+"\".\n"); - return; } } else { if (! optional) { cvm::error("Error: definition for atom group \""+ std::string(group_key)+"\" not found.\n"); - return; } } + return group; +} + + +int colvar::cvc::setup() +{ + size_t i; + description = "cvc " + name; + + for (i = 0; i < atom_groups.size(); i++) { + add_child((cvm::deps *) atom_groups[i]); + } + + if (b_try_scalable && is_available(f_cvc_scalable)) { + enable(f_cvc_scalable); + } + + return COLVARS_OK; } colvar::cvc::~cvc() -{} +{ + remove_all_children(); + for (size_t i = 0; i < atom_groups.size(); i++) { + if (atom_groups[i] != NULL) delete atom_groups[i]; + } +} + + +void colvar::cvc::read_data() +{ + size_t ig; + for (ig = 0; ig < atom_groups.size(); ig++) { + cvm::atom_group &atoms = *(atom_groups[ig]); + atoms.reset_atoms_data(); + atoms.read_positions(); + atoms.calc_required_properties(); + // each atom group will take care of its own ref_pos_group, if defined + } + +//// Don't try to get atom velocities, as no back-end currently implements it +// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) { +// for (i = 0; i < cvcs.size(); i++) { +// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { +// cvcs[i]->atom_groups[ig]->read_velocities(); +// } +// } +// } +} void colvar::cvc::calc_force_invgrads() @@ -86,15 +153,15 @@ void colvar::cvc::calc_Jacobian_derivative() } -void colvar::cvc::debug_gradients(cvm::atom_group &group) +void colvar::cvc::debug_gradients(cvm::atom_group *group) { // this function should work for any scalar variable: // the only difference will be the name of the atom group (here, "group") - if (group.b_dummy) return; + if (group->b_dummy) return; - cvm::rotation const rot_0 = group.rot; - cvm::rotation const rot_inv = group.rot.inverse(); + cvm::rotation const rot_0 = group->rot; + cvm::rotation const rot_inv = group->rot.inverse(); cvm::real x_0 = x.real_value; if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0]; @@ -103,45 +170,45 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) // it only makes sense to debug the fit gradients // when the fitting group is the same as this group - if (group.b_rotate || group.b_center) - if (group.b_fit_gradients && (group.ref_pos_group == NULL)) { - group.calc_fit_gradients(); - if (group.b_rotate) { + if (group->b_rotate || group->b_center) + if (group->b_fit_gradients && (group->ref_pos_group == NULL)) { + group->calc_fit_gradients(); + if (group->b_rotate) { // fit_gradients are in the original frame, we should print them in the rotated frame - for (size_t j = 0; j < group.fit_gradients.size(); j++) { - group.fit_gradients[j] = rot_0.rotate(group.fit_gradients[j]); + for (size_t j = 0; j < group->fit_gradients.size(); j++) { + group->fit_gradients[j] = rot_0.rotate(group->fit_gradients[j]); } } - cvm::log("fit_gradients = "+cvm::to_str(group.fit_gradients)+"\n"); - if (group.b_rotate) { - for (size_t j = 0; j < group.fit_gradients.size(); j++) { - group.fit_gradients[j] = rot_inv.rotate(group.fit_gradients[j]); + cvm::log("fit_gradients = "+cvm::to_str(group->fit_gradients)+"\n"); + if (group->b_rotate) { + for (size_t j = 0; j < group->fit_gradients.size(); j++) { + group->fit_gradients[j] = rot_inv.rotate(group->fit_gradients[j]); } } } - for (size_t ia = 0; ia < group.size(); ia++) { + for (size_t ia = 0; ia < group->size(); ia++) { // tests are best conducted in the unrotated (simulation) frame - cvm::rvector const atom_grad = group.b_rotate ? - rot_inv.rotate(group[ia].grad) : - group[ia].grad; + cvm::rvector const atom_grad = group->b_rotate ? + rot_inv.rotate((*group)[ia].grad) : + (*group)[ia].grad; for (size_t id = 0; id < 3; id++) { // (re)read original positions - group.read_positions(); + group->read_positions(); // change one coordinate - group[ia].pos[id] += cvm::debug_gradients_step_size; - group.calc_required_properties(); + (*group)[ia].pos[id] += cvm::debug_gradients_step_size; + group->calc_required_properties(); calc_value(); cvm::real x_1 = x.real_value; if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0]; cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n"); cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0, 21, 14)+"\n"); - //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ? - cvm::real const dx_pred = (group.fit_gradients.size()) ? - (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : + //cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ? + cvm::real const dx_pred = (group->fit_gradients.size()) ? + (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) : (cvm::debug_gradients_step_size * atom_grad[id]); cvm::log("dx(interp) = "+cvm::to_str(dx_pred, 21, 14)+"\n"); @@ -154,27 +221,27 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) /* * The code below is WIP */ -// if (group.ref_pos_group != NULL) { -// cvm::atom_group &ref = *group.ref_pos_group; -// group.calc_fit_gradients(); +// if (group->ref_pos_group != NULL) { +// cvm::atom_group &ref = *group->ref_pos_group; +// group->calc_fit_gradients(); // // for (size_t ia = 0; ia < ref.size(); ia++) { // // for (size_t id = 0; id < 3; id++) { // // (re)read original positions -// group.read_positions(); +// group->read_positions(); // ref.read_positions(); // // change one coordinate // ref[ia].pos[id] += cvm::debug_gradients_step_size; -// group.update(); +// group->update(); // calc_value(); // cvm::real const x_1 = x.real_value; // cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n"); // cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0, // 21, 14)+"\n"); -// //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ? -// // cvm::real const dx_pred = (group.fit_gradients.size()) ? -// // (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : +// //cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ? +// // cvm::real const dx_pred = (group->fit_gradients.size()) ? +// // (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) : // // (cvm::debug_gradients_step_size * atom_grad[id]); // cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][id]; // cvm::log("dx(interp) = "+cvm::to_str (dx_pred, @@ -190,3 +257,8 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) return; } + + +// Static members + +std::vector colvar::cvc::cvc_features; diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 2f0e4b43b7..923b2c6bd1 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -61,7 +61,7 @@ /// call to e.g. apply_force(). class colvar::cvc - : public colvarparse + : public colvarparse, public cvm::deps { public: @@ -84,10 +84,6 @@ public: /// \brief Exponent in the polynomial combination (default: 1) int sup_np; - /// \brief This defaults to true; setting it to false disables - /// update of this cvc to save compute time (useful with scriptedFunction) - bool b_enabled; - /// \brief Is this a periodic component? bool b_periodic; @@ -107,11 +103,14 @@ public: /// \brief Within the constructor, make a group parse its own /// options from the provided configuration string - void parse_group(std::string const &conf, + /// Returns reference to new group + cvm::atom_group *parse_group(std::string const &conf, char const *group_key, - cvm::atom_group &group, bool optional = false); + /// \brief After construction, set data related to dependency handling + int setup(); + /// \brief Default constructor (used when \link cvc \endlink /// objects are declared within other ones) cvc(); @@ -119,18 +118,16 @@ public: /// Destructor virtual ~cvc(); + /// \brief Implementation of the feature list for colvar + static std::vector cvc_features; - /// \brief If this flag is false (default), inverse gradients - /// (derivatives of atom coordinates with respect to x) are - /// unavailable; it should be set to true by the constructor of each - /// derived object capable of calculating them - bool b_inverse_gradients; + /// \brief Implementation of the feature list accessor for colvar + virtual std::vector &features() { + return cvc_features; + } - /// \brief If this flag is false (default), the Jacobian derivative - /// (divergence of the inverse gradients) is unavailable; it should - /// be set to true by the constructor of each derived object capable - /// of calculating it - bool b_Jacobian_derivative; + /// \brief Obtain data needed for the calculation for the backend + void read_data(); /// \brief Calculate the variable virtual void calc_value() = 0; @@ -139,11 +136,8 @@ public: /// order to apply forces virtual void calc_gradients() = 0; - /// \brief If true, calc_gradients() will call debug_gradients() for every group needed - bool b_debug_gradients; - /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component - virtual void debug_gradients(cvm::atom_group &group); + virtual void debug_gradients(cvm::atom_group *group); /// \brief Calculate the total force from the system using the /// inverse atomic gradients @@ -227,6 +221,9 @@ public: /// e.g. atomic gradients std::vector atom_groups; + /// \brief Whether or not this CVC will be computed in parallel whenever possible + bool b_try_scalable; + protected: /// \brief Cached value @@ -307,9 +304,9 @@ class colvar::distance { protected: /// First atom group - cvm::atom_group group1; + cvm::atom_group *group1; /// Second atom group - cvm::atom_group group2; + cvm::atom_group *group2; /// Vector distance, cached to be recycled cvm::rvector dist_v; /// Use absolute positions, ignoring PBCs when present @@ -389,11 +386,11 @@ class colvar::distance_z { protected: /// Main atom group - cvm::atom_group main; + cvm::atom_group *main; /// Reference atom group - cvm::atom_group ref1; + cvm::atom_group *ref1; /// Optional, second ref atom group - cvm::atom_group ref2; + cvm::atom_group *ref2; /// Use absolute positions, ignoring PBCs when present bool b_no_PBC; /// Compute system force on one site only to avoid unwanted @@ -486,9 +483,9 @@ class colvar::distance_pairs { protected: /// First atom group - cvm::atom_group group1; + cvm::atom_group *group1; /// Second atom group - cvm::atom_group group2; + cvm::atom_group *group2; /// Use absolute positions, ignoring PBCs when present bool b_no_PBC; public: @@ -515,7 +512,7 @@ class colvar::gyration { protected: /// Atoms involved - cvm::atom_group atoms; + cvm::atom_group *atoms; public: /// Constructor gyration(std::string const &conf); @@ -590,7 +587,7 @@ class colvar::eigenvector protected: /// Atom group - cvm::atom_group atoms; + cvm::atom_group * atoms; /// Reference coordinates std::vector ref_pos; @@ -632,11 +629,11 @@ class colvar::angle protected: /// Atom group - cvm::atom_group group1; + cvm::atom_group *group1; /// Atom group - cvm::atom_group group2; + cvm::atom_group *group2; /// Atom group - cvm::atom_group group3; + cvm::atom_group *group3; /// Inter site vectors cvm::rvector r21, r23; @@ -678,11 +675,11 @@ class colvar::dipole_angle protected: /// Dipole atom group - cvm::atom_group group1; + cvm::atom_group *group1; /// Atom group - cvm::atom_group group2; + cvm::atom_group *group2; /// Atom group - cvm::atom_group group3; + cvm::atom_group *group3; /// Inter site vectors cvm::rvector r21, r23; @@ -722,13 +719,13 @@ class colvar::dihedral protected: /// Atom group - cvm::atom_group group1; + cvm::atom_group *group1; /// Atom group - cvm::atom_group group2; + cvm::atom_group *group2; /// Atom group - cvm::atom_group group3; + cvm::atom_group *group3; /// Atom group - cvm::atom_group group4; + cvm::atom_group *group4; /// Inter site vectors cvm::rvector r12, r23, r34; @@ -1009,7 +1006,7 @@ class colvar::orientation protected: /// Atom group - cvm::atom_group atoms; + cvm::atom_group * atoms; /// Center of geometry of the group cvm::atom_pos atoms_cog; @@ -1155,7 +1152,7 @@ class colvar::rmsd protected: /// Atom group - cvm::atom_group atoms; + cvm::atom_group *atoms; /// Reference coordinates (for RMSD calculation only) std::vector ref_pos; @@ -1186,7 +1183,7 @@ class colvar::cartesian { protected: /// Atom group - cvm::atom_group atoms; + cvm::atom_group *atoms; /// Which Cartesian coordinates to include std::vector axes; public: diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 35b98c0e46..519288a01f 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -11,14 +11,13 @@ colvar::angle::angle(std::string const &conf) : cvc(conf) { function_type = "angle"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; - parse_group(conf, "group1", group1); - parse_group(conf, "group2", group2); - parse_group(conf, "group3", group3); - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); + + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); + group3 = parse_group(conf, "group3"); if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { cvm::log("Computing system force on group 1 only"); } @@ -27,19 +26,21 @@ colvar::angle::angle(std::string const &conf) colvar::angle::angle(cvm::atom const &a1, - cvm::atom const &a2, - cvm::atom const &a3) - : group1(std::vector(1, a1)), - group2(std::vector(1, a2)), - group3(std::vector(1, a3)) + cvm::atom const &a2, + cvm::atom const &a3) { function_type = "angle"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); b_1site_force = false; - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); + + group1 = new cvm::atom_group(std::vector(1, a1)); + group2 = new cvm::atom_group(std::vector(1, a2)); + group3 = new cvm::atom_group(std::vector(1, a3)); + atom_groups.push_back(group1); + atom_groups.push_back(group2); + atom_groups.push_back(group3); x.type(colvarvalue::type_scalar); } @@ -54,9 +55,9 @@ colvar::angle::angle() void colvar::angle::calc_value() { - cvm::atom_pos const g1_pos = group1.center_of_mass(); - cvm::atom_pos const g2_pos = group2.center_of_mass(); - cvm::atom_pos const g3_pos = group3.center_of_mass(); + cvm::atom_pos const g1_pos = group1->center_of_mass(); + cvm::atom_pos const g2_pos = group2->center_of_mass(); + cvm::atom_pos const g3_pos = group3->center_of_mass(); r21 = cvm::position_distance(g2_pos, g1_pos); r21l = r21.norm(); @@ -71,7 +72,6 @@ void colvar::angle::calc_value() void colvar::angle::calc_gradients() { - size_t i; cvm::real const cos_theta = (r21*r23)/(r21l*r23l); cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta); @@ -81,22 +81,11 @@ void colvar::angle::calc_gradients() dxdr3 = (180.0/PI) * dxdcos * (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); - for (i = 0; i < group1.size(); i++) { - group1[i].grad = (group1[i].mass/group1.total_mass) * - (dxdr1); - } + group1->set_weighted_gradient(dxdr1); + group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0)); + group3->set_weighted_gradient(dxdr3); - for (i = 0; i < group2.size(); i++) { - group2[i].grad = (group2[i].mass/group2.total_mass) * - (dxdr1 + dxdr3) * (-1.0); - } - - for (i = 0; i < group3.size(); i++) { - group3[i].grad = (group3[i].mass/group3.total_mass) * - (dxdr3); - } - - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { debug_gradients(group1); debug_gradients(group2); debug_gradients(group3); @@ -112,15 +101,15 @@ void colvar::angle::calc_force_invgrads() // when propagating changes in the angle) if (b_1site_force) { - group1.read_system_forces(); + group1->read_system_forces(); cvm::real norm_fact = 1.0 / dxdr1.norm2(); - ft.real_value = norm_fact * dxdr1 * group1.system_force(); + ft.real_value = norm_fact * dxdr1 * group1->system_force(); } else { - group1.read_system_forces(); - group3.read_system_forces(); + group1->read_system_forces(); + group3->read_system_forces(); cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2()); - ft.real_value = norm_fact * ( dxdr1 * group1.system_force() - + dxdr3 * group3.system_force()); + ft.real_value = norm_fact * ( dxdr1 * group1->system_force() + + dxdr3 * group3->system_force()); } return; } @@ -136,14 +125,14 @@ void colvar::angle::calc_Jacobian_derivative() void colvar::angle::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); - if (!group3.noforce) - group3.apply_colvar_force(force.real_value); + if (!group3->noforce) + group3->apply_colvar_force(force.real_value); } @@ -153,13 +142,10 @@ colvar::dipole_angle::dipole_angle(std::string const &conf) : cvc(conf) { function_type = "dipole_angle"; - parse_group(conf, "group1", group1); - parse_group(conf, "group2", group2); - parse_group(conf, "group3", group3); + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); + group3 = parse_group(conf, "group3"); - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { cvm::log("Computing system force on group 1 only"); } @@ -170,15 +156,16 @@ colvar::dipole_angle::dipole_angle(std::string const &conf) colvar::dipole_angle::dipole_angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3) - : group1(std::vector(1, a1)), - group2(std::vector(1, a2)), - group3(std::vector(1, a3)) { function_type = "dipole_angle"; b_1site_force = false; - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); + + group1 = new cvm::atom_group(std::vector(1, a1)); + group2 = new cvm::atom_group(std::vector(1, a2)); + group3 = new cvm::atom_group(std::vector(1, a3)); + atom_groups.push_back(group1); + atom_groups.push_back(group2); + atom_groups.push_back(group3); x.type(colvarvalue::type_scalar); } @@ -193,13 +180,13 @@ colvar::dipole_angle::dipole_angle() void colvar::dipole_angle::calc_value() { - cvm::atom_pos const g1_pos = group1.center_of_mass(); - cvm::atom_pos const g2_pos = group2.center_of_mass(); - cvm::atom_pos const g3_pos = group3.center_of_mass(); + cvm::atom_pos const g1_pos = group1->center_of_mass(); + cvm::atom_pos const g2_pos = group2->center_of_mass(); + cvm::atom_pos const g3_pos = group3->center_of_mass(); - group1.calc_dipole(g1_pos); + group1->calc_dipole(g1_pos); - r21 = group1.dipole(); + r21 = group1->dipole(); r21l = r21.norm(); r23 = cvm::position_distance(g2_pos, g3_pos); r23l = r23.norm(); @@ -225,35 +212,35 @@ void colvar::dipole_angle::calc_gradients() (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); //this auxiliar variables are to avoid numerical errors inside "for" - double aux1 = group1.total_charge/group1.total_mass; - // double aux2 = group2.total_charge/group2.total_mass; - // double aux3 = group3.total_charge/group3.total_mass; + double aux1 = group1->total_charge/group1->total_mass; + // double aux2 = group2->total_charge/group2->total_mass; + // double aux3 = group3->total_charge/group3->total_mass; size_t i; - for (i = 0; i < group1.size(); i++) { - group1[i].grad =(group1[i].charge + (-1)* group1[i].mass * aux1) * (dxdr1); + for (i = 0; i < group1->size(); i++) { + (*group1)[i].grad =((*group1)[i].charge + (-1)* (*group1)[i].mass * aux1) * (dxdr1); } - for (i = 0; i < group2.size(); i++) { - group2[i].grad = (group2[i].mass/group2.total_mass)* dxdr3 * (-1.0); + for (i = 0; i < group2->size(); i++) { + (*group2)[i].grad = ((*group2)[i].mass/group2->total_mass)* dxdr3 * (-1.0); } - for (i = 0; i < group3.size(); i++) { - group3[i].grad =(group3[i].mass/group3.total_mass) * (dxdr3); + for (i = 0; i < group3->size(); i++) { + (*group3)[i].grad =((*group3)[i].mass/group3->total_mass) * (dxdr3); } } void colvar::dipole_angle::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); - if (!group3.noforce) - group3.apply_colvar_force(force.real_value); + if (!group3->noforce) + group3->apply_colvar_force(force.real_value); } @@ -265,32 +252,26 @@ colvar::dihedral::dihedral(std::string const &conf) function_type = "dihedral"; period = 360.0; b_periodic = true; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); + if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { cvm::log("Computing system force on group 1 only"); } - parse_group(conf, "group1", group1); - parse_group(conf, "group2", group2); - parse_group(conf, "group3", group3); - parse_group(conf, "group4", group4); - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); - atom_groups.push_back(&group4); + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); + group3 = parse_group(conf, "group3"); + group4 = parse_group(conf, "group4"); x.type(colvarvalue::type_scalar); } colvar::dihedral::dihedral(cvm::atom const &a1, - cvm::atom const &a2, - cvm::atom const &a3, - cvm::atom const &a4) - : group1(std::vector(1, a1)), - group2(std::vector(1, a2)), - group3(std::vector(1, a3)), - group4(std::vector(1, a4)) + cvm::atom const &a2, + cvm::atom const &a3, + cvm::atom const &a4) { if (cvm::debug()) cvm::log("Initializing dihedral object from atom groups.\n"); @@ -298,14 +279,20 @@ colvar::dihedral::dihedral(cvm::atom const &a1, function_type = "dihedral"; period = 360.0; b_periodic = true; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); + b_1site_force = false; - atom_groups.push_back(&group1); - atom_groups.push_back(&group2); - atom_groups.push_back(&group3); - atom_groups.push_back(&group4); + group1 = new cvm::atom_group(std::vector(1, a1)); + group2 = new cvm::atom_group(std::vector(1, a2)); + group3 = new cvm::atom_group(std::vector(1, a3)); + group4 = new cvm::atom_group(std::vector(1, a4)); + atom_groups.push_back(group1); + atom_groups.push_back(group2); + atom_groups.push_back(group3); + atom_groups.push_back(group4); x.type(colvarvalue::type_scalar); @@ -319,18 +306,18 @@ colvar::dihedral::dihedral() function_type = "dihedral"; period = 360.0; b_periodic = true; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); x.type(colvarvalue::type_scalar); } void colvar::dihedral::calc_value() { - cvm::atom_pos const g1_pos = group1.center_of_mass(); - cvm::atom_pos const g2_pos = group2.center_of_mass(); - cvm::atom_pos const g3_pos = group3.center_of_mass(); - cvm::atom_pos const g4_pos = group4.center_of_mass(); + cvm::atom_pos const g1_pos = group1->center_of_mass(); + cvm::atom_pos const g2_pos = group2->center_of_mass(); + cvm::atom_pos const g3_pos = group3->center_of_mass(); + cvm::atom_pos const g4_pos = group4->center_of_mass(); // Usual sign convention: r12 = r2 - r1 r12 = cvm::position_distance(g1_pos, g2_pos); @@ -415,18 +402,10 @@ void colvar::dihedral::calc_gradients() +dsindB.y*r34.x - dsindB.x*r34.y); } - size_t i; - for (i = 0; i < group1.size(); i++) - group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1); - - for (i = 0; i < group2.size(); i++) - group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1); - - for (i = 0; i < group3.size(); i++) - group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2); - - for (i = 0; i < group4.size(); i++) - group4[i].grad = (group4[i].mass/group4.total_mass) * (f3); + group1->set_weighted_gradient(-f1); + group2->set_weighted_gradient(-f2 + f1); + group3->set_weighted_gradient(-f3 + f2); + group4->set_weighted_gradient(f3); } @@ -448,15 +427,15 @@ void colvar::dihedral::calc_force_invgrads() cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1); cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4); - group1.read_system_forces(); + group1->read_system_forces(); if ( b_1site_force ) { // This is only measuring the force on group 1 - ft.real_value = PI/180.0 * fact1 * (cross1 * group1.system_force()); + ft.real_value = PI/180.0 * fact1 * (cross1 * group1->system_force()); } else { // Default case: use groups 1 and 4 - group4.read_system_forces(); - ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1.system_force()) - + fact4 * (cross4 * group4.system_force())); + group4->read_system_forces(); + ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1->system_force()) + + fact4 * (cross4 * group4->system_force())); } } @@ -470,17 +449,17 @@ void colvar::dihedral::calc_Jacobian_derivative() void colvar::dihedral::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); - if (!group3.noforce) - group3.apply_colvar_force(force.real_value); + if (!group3->noforce) + group3->apply_colvar_force(force.real_value); - if (!group4.noforce) - group4.apply_colvar_force(force.real_value); + if (!group4->noforce) + group4->apply_colvar_force(force.real_value); } diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index 7541c8bb23..2d4a20d542 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -79,13 +79,13 @@ colvar::coordnum::coordnum(std::string const &conf) x.type(colvarvalue::type_scalar); // group1 and group2 are already initialized by distance() - if (group1.b_dummy) + if (group1->b_dummy) cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n"); // need to specify this explicitly because the distance() constructor // has set it to true - b_inverse_gradients = false; + feature_states[f_cvc_inv_gradient]->available = false; bool const b_scale = get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom())); @@ -110,7 +110,7 @@ colvar::coordnum::coordnum(std::string const &conf) cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n"); } - get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy); + get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy); } @@ -130,26 +130,26 @@ void colvar::coordnum::calc_value() // create a fake atom to hold the group2 com coordinates cvm::atom group2_com_atom; - group2_com_atom.pos = group2.center_of_mass(); + group2_com_atom.pos = group2->center_of_mass(); if (b_anisotropic) { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) x.real_value += switching_function(r0_vec, en, ed, *ai1, group2_com_atom); } else { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) x.real_value += switching_function(r0, en, ed, *ai1, group2_com_atom); } } else { if (b_anisotropic) { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { x.real_value += switching_function(r0_vec, en, ed, *ai1, *ai2); } } else { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { x.real_value += switching_function(r0, en, ed, *ai1, *ai2); } } @@ -163,29 +163,29 @@ void colvar::coordnum::calc_gradients() // create a fake atom to hold the group2 com coordinates cvm::atom group2_com_atom; - group2_com_atom.pos = group2.center_of_mass(); + group2_com_atom.pos = group2->center_of_mass(); if (b_anisotropic) { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) switching_function(r0_vec, en, ed, *ai1, group2_com_atom); } else { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) switching_function(r0, en, ed, *ai1, group2_com_atom); } - group2.set_weighted_gradient(group2_com_atom.grad); + group2->set_weighted_gradient(group2_com_atom.grad); } else { if (b_anisotropic) { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { switching_function(r0_vec, en, ed, *ai1, *ai2); } } else { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { switching_function(r0, en, ed, *ai1, *ai2); } } @@ -194,11 +194,11 @@ void colvar::coordnum::calc_gradients() void colvar::coordnum::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); } @@ -298,7 +298,8 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf) // need to specify this explicitly because the distance() constructor // has set it to true - b_inverse_gradients = false; + feature_states[f_cvc_inv_gradient]->available = false; + get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom())); get_keyval(conf, "expNumer", en, int(6) ); @@ -320,9 +321,9 @@ colvar::selfcoordnum::selfcoordnum() void colvar::selfcoordnum::calc_value() { x.real_value = 0.0; - for (size_t i = 0; i < group1.size() - 1; i++) { - for (size_t j = i + 1; j < group1.size(); j++) { - x.real_value += colvar::coordnum::switching_function(r0, en, ed, group1[i], group1[j]); + for (size_t i = 0; i < group1->size() - 1; i++) { + for (size_t j = i + 1; j < group1->size(); j++) { + x.real_value += colvar::coordnum::switching_function(r0, en, ed, (*group1)[i], (*group1)[j]); } } } @@ -330,17 +331,17 @@ void colvar::selfcoordnum::calc_value() void colvar::selfcoordnum::calc_gradients() { - for (size_t i = 0; i < group1.size() - 1; i++) { - for (size_t j = i + 1; j < group1.size(); j++) { - colvar::coordnum::switching_function(r0, en, ed, group1[i], group1[j]); + for (size_t i = 0; i < group1->size() - 1; i++) { + for (size_t j = i + 1; j < group1->size(); j++) { + colvar::coordnum::switching_function(r0, en, ed, (*group1)[i], (*group1)[j]); } } } void colvar::selfcoordnum::apply_force(colvarvalue const &force) { - if (!group1.noforce) { - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) { + group1->apply_colvar_force(force.real_value); } } diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 3531cfc234..8d15c5d67a 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -17,19 +17,19 @@ colvar::distance::distance(std::string const &conf, bool twogroups) : cvc(conf) { function_type = "distance"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); + if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { cvm::log("Computing distance using absolute positions (not minimal-image)"); } if (twogroups && get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { cvm::log("Computing system force on group 1 only"); } - parse_group(conf, "group1", group1); - atom_groups.push_back(&group1); + group1 = parse_group(conf, "group1"); if (twogroups) { - parse_group(conf, "group2", group2); - atom_groups.push_back(&group2); + group2 = parse_group(conf, "group2"); } x.type(colvarvalue::type_scalar); } @@ -39,54 +39,60 @@ colvar::distance::distance() : cvc() { function_type = "distance"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); b_1site_force = false; b_no_PBC = false; x.type(colvarvalue::type_scalar); } + void colvar::distance::calc_value() { if (b_no_PBC) { - dist_v = group2.center_of_mass() - group1.center_of_mass(); + dist_v = group2->center_of_mass() - group1->center_of_mass(); } else { - dist_v = cvm::position_distance(group1.center_of_mass(), - group2.center_of_mass()); + dist_v = cvm::position_distance(group1->center_of_mass(), + group2->center_of_mass()); } x.real_value = dist_v.norm(); } + void colvar::distance::calc_gradients() { cvm::rvector const u = dist_v.unit(); - group1.set_weighted_gradient(-1.0 * u); - group2.set_weighted_gradient( u); + group1->set_weighted_gradient(-1.0 * u); + group2->set_weighted_gradient( u); } + void colvar::distance::calc_force_invgrads() { - group1.read_system_forces(); + group1->read_system_forces(); if ( b_1site_force ) { - ft.real_value = -1.0 * (group1.system_force() * dist_v.unit()); + ft.real_value = -1.0 * (group1->system_force() * dist_v.unit()); } else { - group2.read_system_forces(); - ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit()); + group2->read_system_forces(); + ft.real_value = 0.5 * ((group2->system_force() - group1->system_force()) * dist_v.unit()); } } + void colvar::distance::calc_Jacobian_derivative() { jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0; } + void colvar::distance::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); } @@ -98,6 +104,7 @@ colvar::distance_vec::distance_vec(std::string const &conf) x.type(colvarvalue::type_3vector); } + colvar::distance_vec::distance_vec() : distance() { @@ -105,29 +112,32 @@ colvar::distance_vec::distance_vec() x.type(colvarvalue::type_3vector); } + void colvar::distance_vec::calc_value() { if (b_no_PBC) { - x.rvector_value = group2.center_of_mass() - group1.center_of_mass(); + x.rvector_value = group2->center_of_mass() - group1->center_of_mass(); } else { - x.rvector_value = cvm::position_distance(group1.center_of_mass(), - group2.center_of_mass()); + x.rvector_value = cvm::position_distance(group1->center_of_mass(), + group2->center_of_mass()); } } + void colvar::distance_vec::calc_gradients() { // gradients are not stored: a 3x3 matrix for each atom would be // needed to store just the identity matrix } + void colvar::distance_vec::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_force(-1.0 * force.rvector_value); + if (!group1->noforce) + group1->apply_force(-1.0 * force.rvector_value); - if (!group2.noforce) - group2.apply_force( force.rvector_value); + if (!group2->noforce) + group2->apply_force( force.rvector_value); } @@ -136,8 +146,9 @@ colvar::distance_z::distance_z(std::string const &conf) : cvc(conf) { function_type = "distance_z"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); x.type(colvarvalue::type_scalar); // TODO detect PBC from MD engine (in simple cases) @@ -151,15 +162,12 @@ colvar::distance_z::distance_z(std::string const &conf) return; } - parse_group(conf, "main", main); - parse_group(conf, "ref", ref1); - atom_groups.push_back(&main); - atom_groups.push_back(&ref1); + main = parse_group(conf, "main"); + ref1 = parse_group(conf, "ref"); // this group is optional - parse_group(conf, "ref2", ref2, true); + ref2 = parse_group(conf, "ref2", true); - if (ref2.size()) { - atom_groups.push_back(&ref2); + if (ref2->size()) { cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); fixed_axis = false; if (key_lookup(conf, "axis")) @@ -189,8 +197,9 @@ colvar::distance_z::distance_z(std::string const &conf) colvar::distance_z::distance_z() { function_type = "distance_z"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); x.type(colvarvalue::type_scalar); } @@ -198,21 +207,21 @@ void colvar::distance_z::calc_value() { if (fixed_axis) { if (b_no_PBC) { - dist_v = main.center_of_mass() - ref1.center_of_mass(); + dist_v = main->center_of_mass() - ref1->center_of_mass(); } else { - dist_v = cvm::position_distance(ref1.center_of_mass(), - main.center_of_mass()); + dist_v = cvm::position_distance(ref1->center_of_mass(), + main->center_of_mass()); } } else { if (b_no_PBC) { - dist_v = main.center_of_mass() - - (0.5 * (ref1.center_of_mass() + ref2.center_of_mass())); - axis = ref2.center_of_mass() - ref1.center_of_mass(); + dist_v = main->center_of_mass() - + (0.5 * (ref1->center_of_mass() + ref2->center_of_mass())); + axis = ref2->center_of_mass() - ref1->center_of_mass(); } else { - dist_v = cvm::position_distance(0.5 * (ref1.center_of_mass() + - ref2.center_of_mass()), main.center_of_mass()); - axis = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass()); + dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() + + ref2->center_of_mass()), main->center_of_mass()); + axis = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass()); } axis_norm = axis.norm(); axis = axis.unit(); @@ -223,25 +232,25 @@ void colvar::distance_z::calc_value() void colvar::distance_z::calc_gradients() { - main.set_weighted_gradient( axis ); + main->set_weighted_gradient( axis ); if (fixed_axis) { - ref1.set_weighted_gradient(-1.0 * axis); + ref1->set_weighted_gradient(-1.0 * axis); } else { if (b_no_PBC) { - ref1.set_weighted_gradient( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() - + ref1->set_weighted_gradient( 1.0 / axis_norm * (main->center_of_mass() - ref2->center_of_mass() - x.real_value * axis )); - ref2.set_weighted_gradient( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() + + ref2->set_weighted_gradient( 1.0 / axis_norm * (ref1->center_of_mass() - main->center_of_mass() + x.real_value * axis )); } else { - ref1.set_weighted_gradient( 1.0 / axis_norm * ( - cvm::position_distance(ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis )); - ref2.set_weighted_gradient( 1.0 / axis_norm * ( - cvm::position_distance(main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis )); + ref1->set_weighted_gradient( 1.0 / axis_norm * ( + cvm::position_distance(ref2->center_of_mass(), main->center_of_mass()) - x.real_value * axis )); + ref2->set_weighted_gradient( 1.0 / axis_norm * ( + cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis )); } } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients for group main:\n"); debug_gradients(main); cvm::log("Debugging gradients for group ref1:\n"); @@ -253,13 +262,13 @@ void colvar::distance_z::calc_gradients() void colvar::distance_z::calc_force_invgrads() { - main.read_system_forces(); + main->read_system_forces(); if (fixed_axis && !b_1site_force) { - ref1.read_system_forces(); - ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis); + ref1->read_system_forces(); + ft.real_value = 0.5 * ((main->system_force() - ref1->system_force()) * axis); } else { - ft.real_value = main.system_force() * axis; + ft.real_value = main->system_force() * axis; } } @@ -270,14 +279,14 @@ void colvar::distance_z::calc_Jacobian_derivative() void colvar::distance_z::apply_force(colvarvalue const &force) { - if (!ref1.noforce) - ref1.apply_colvar_force(force.real_value); + if (!ref1->noforce) + ref1->apply_colvar_force(force.real_value); - if (ref2.size() && !ref2.noforce) - ref2.apply_colvar_force(force.real_value); + if (ref2->size() && !ref2->noforce) + ref2->apply_colvar_force(force.real_value); - if (!main.noforce) - main.apply_colvar_force(force.real_value); + if (!main->noforce) + main->apply_colvar_force(force.real_value); } @@ -286,8 +295,9 @@ colvar::distance_xy::distance_xy(std::string const &conf) : distance_z(conf) { function_type = "distance_xy"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); x.type(colvarvalue::type_scalar); } @@ -295,24 +305,25 @@ colvar::distance_xy::distance_xy() : distance_z() { function_type = "distance_xy"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + provide(f_cvc_com_based); x.type(colvarvalue::type_scalar); } void colvar::distance_xy::calc_value() { if (b_no_PBC) { - dist_v = main.center_of_mass() - ref1.center_of_mass(); + dist_v = main->center_of_mass() - ref1->center_of_mass(); } else { - dist_v = cvm::position_distance(ref1.center_of_mass(), - main.center_of_mass()); + dist_v = cvm::position_distance(ref1->center_of_mass(), + main->center_of_mass()); } if (!fixed_axis) { if (b_no_PBC) { - v12 = ref2.center_of_mass() - ref1.center_of_mass(); + v12 = ref2->center_of_mass() - ref1->center_of_mass(); } else { - v12 = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass()); + v12 = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass()); } axis_norm = v12.norm(); axis = v12.unit(); @@ -333,31 +344,31 @@ void colvar::distance_xy::calc_gradients() x_inv = 1.0 / x.real_value; if (fixed_axis) { - ref1.set_weighted_gradient(-1.0 * x_inv * dist_v_ortho); - main.set_weighted_gradient( x_inv * dist_v_ortho); + ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho); + main->set_weighted_gradient( x_inv * dist_v_ortho); } else { if (b_no_PBC) { - v13 = main.center_of_mass() - ref1.center_of_mass(); + v13 = main->center_of_mass() - ref1->center_of_mass(); } else { - v13 = cvm::position_distance(ref1.center_of_mass(), main.center_of_mass()); + v13 = cvm::position_distance(ref1->center_of_mass(), main->center_of_mass()); } A = (dist_v * axis) / axis_norm; - ref1.set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho); - ref2.set_weighted_gradient( -A * x_inv * dist_v_ortho); - main.set_weighted_gradient( 1.0 * x_inv * dist_v_ortho); + ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho); + ref2->set_weighted_gradient( -A * x_inv * dist_v_ortho); + main->set_weighted_gradient( 1.0 * x_inv * dist_v_ortho); } } void colvar::distance_xy::calc_force_invgrads() { - main.read_system_forces(); + main->read_system_forces(); if (fixed_axis && !b_1site_force) { - ref1.read_system_forces(); - ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho); + ref1->read_system_forces(); + ft.real_value = 0.5 / x.real_value * ((main->system_force() - ref1->system_force()) * dist_v_ortho); } else { - ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho; + ft.real_value = 1.0 / x.real_value * main->system_force() * dist_v_ortho; } } @@ -368,14 +379,14 @@ void colvar::distance_xy::calc_Jacobian_derivative() void colvar::distance_xy::apply_force(colvarvalue const &force) { - if (!ref1.noforce) - ref1.apply_colvar_force(force.real_value); + if (!ref1->noforce) + ref1->apply_colvar_force(force.real_value); - if (ref2.size() && !ref2.noforce) - ref2.apply_colvar_force(force.real_value); + if (ref2->size() && !ref2->noforce) + ref2->apply_colvar_force(force.real_value); - if (!main.noforce) - main.apply_colvar_force(force.real_value); + if (!main->noforce) + main->apply_colvar_force(force.real_value); } @@ -399,10 +410,10 @@ colvar::distance_dir::distance_dir() void colvar::distance_dir::calc_value() { if (b_no_PBC) { - dist_v = group2.center_of_mass() - group1.center_of_mass(); + dist_v = group2->center_of_mass() - group1->center_of_mass(); } else { - dist_v = cvm::position_distance(group1.center_of_mass(), - group2.center_of_mass()); + dist_v = cvm::position_distance(group1->center_of_mass(), + group2->center_of_mass()); } x.rvector_value = dist_v.unit(); } @@ -413,6 +424,8 @@ void colvar::distance_dir::calc_gradients() // gradients are computed on the fly within apply_force() // Note: could be a problem if a future bias relies on gradient // calculations... + // TODO in new deps system: remove dependency of biasing force to gradient? + // That way we could tell apart an explicit gradient dependency } @@ -422,11 +435,11 @@ void colvar::distance_dir::apply_force(colvarvalue const &force) cvm::real const iprod = force.rvector_value * x.rvector_value; cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value; - if (!group1.noforce) - group1.apply_force(-1.0 * force_tang); + if (!group1->noforce) + group1->apply_force(-1.0 * force_tang); - if (!group2.noforce) - group2.apply_force( force_tang); + if (!group2->noforce) + group2->apply_force( force_tang); } @@ -445,8 +458,8 @@ colvar::distance_inv::distance_inv(std::string const &conf) return; } - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { if (ai1->id == ai2->id) { cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n"); return; @@ -454,8 +467,6 @@ colvar::distance_inv::distance_inv(std::string const &conf) } } - b_inverse_gradients = false; - b_Jacobian_derivative = false; x.type(colvarvalue::type_scalar); } @@ -463,8 +474,6 @@ colvar::distance_inv::distance_inv() { function_type = "distance_inv"; exponent = 6; - b_inverse_gradients = false; - b_Jacobian_derivative = false; b_1site_force = false; x.type(colvarvalue::type_scalar); } @@ -473,8 +482,8 @@ void colvar::distance_inv::calc_value() { x.real_value = 0.0; if (b_no_PBC) { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { cvm::rvector const dv = ai2->pos - ai1->pos; cvm::real const d2 = dv.norm2(); cvm::real dinv = 1.0; @@ -487,8 +496,8 @@ void colvar::distance_inv::calc_value() } } } else { - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos); cvm::real const d2 = dv.norm2(); cvm::real dinv = 1.0; @@ -502,28 +511,28 @@ void colvar::distance_inv::calc_value() } } - x.real_value *= 1.0 / cvm::real(group1.size() * group2.size()); + x.real_value *= 1.0 / cvm::real(group1->size() * group2->size()); x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent))); } void colvar::distance_inv::calc_gradients() { - cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1.size() * group2.size()); - for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { + cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size()); + for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) { ai1->grad *= dxdsum; } - for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { + for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) { ai2->grad *= dxdsum; } } void colvar::distance_inv::apply_force(colvarvalue const &force) { - if (!group1.noforce) - group1.apply_colvar_force(force.real_value); + if (!group1->noforce) + group1->apply_colvar_force(force.real_value); - if (!group2.noforce) - group2.apply_colvar_force(force.real_value); + if (!group2->noforce) + group2->apply_colvar_force(force.real_value); } @@ -532,57 +541,50 @@ colvar::distance_pairs::distance_pairs(std::string const &conf) : cvc(conf) { function_type = "distance_pairs"; - b_inverse_gradients = false; - b_Jacobian_derivative = false; if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { cvm::log("Computing distance using absolute positions (not minimal-image)"); } - parse_group(conf, "group1", group1); - atom_groups.push_back(&group1); - - parse_group(conf, "group2", group2); - atom_groups.push_back(&group2); + group1 = parse_group(conf, "group1"); + group2 = parse_group(conf, "group2"); x.type(colvarvalue::type_vector); - x.vector1d_value.resize(group1.size() * group2.size()); + x.vector1d_value.resize(group1->size() * group2->size()); } colvar::distance_pairs::distance_pairs() { function_type = "distance_pairs"; - b_inverse_gradients = false; - b_Jacobian_derivative = false; x.type(colvarvalue::type_vector); } void colvar::distance_pairs::calc_value() { - x.vector1d_value.resize(group1.size() * group2.size()); + x.vector1d_value.resize(group1->size() * group2->size()); if (b_no_PBC) { size_t i1, i2; - for (i1 = 0; i1 < group1.size(); i1++) { - for (i2 = 0; i2 < group2.size(); i2++) { - cvm::rvector const dv = group2[i2].pos - group1[i1].pos; + for (i1 = 0; i1 < group1->size(); i1++) { + for (i2 = 0; i2 < group2->size(); i2++) { + cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos; cvm::real const d = dv.norm(); - x.vector1d_value[i1*group2.size() + i2] = d; - group1[i1].grad = -1.0 * dv.unit(); - group2[i2].grad = dv.unit(); + x.vector1d_value[i1*group2->size() + i2] = d; + (*group1)[i1].grad = -1.0 * dv.unit(); + (*group2)[i2].grad = dv.unit(); } } } else { size_t i1, i2; - for (i1 = 0; i1 < group1.size(); i1++) { - for (i2 = 0; i2 < group2.size(); i2++) { - cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); + for (i1 = 0; i1 < group1->size(); i1++) { + for (i2 = 0; i2 < group2->size(); i2++) { + cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos); cvm::real const d = dv.norm(); - x.vector1d_value[i1*group2.size() + i2] = d; - group1[i1].grad = -1.0 * dv.unit(); - group2[i2].grad = dv.unit(); + x.vector1d_value[i1*group2->size() + i2] = d; + (*group1)[i1].grad = -1.0 * dv.unit(); + (*group2)[i2].grad = dv.unit(); } } } @@ -591,7 +593,7 @@ void colvar::distance_pairs::calc_value() void colvar::distance_pairs::calc_gradients() { // will be calculated on the fly in apply_force() - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(group1); } @@ -601,20 +603,20 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force) { if (b_no_PBC) { size_t i1, i2; - for (i1 = 0; i1 < group1.size(); i1++) { - for (i2 = 0; i2 < group2.size(); i2++) { - cvm::rvector const dv = group2[i2].pos - group1[i1].pos; - group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); - group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); + for (i1 = 0; i1 < group1->size(); i1++) { + for (i2 = 0; i2 < group2->size(); i2++) { + cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos; + (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit()); + (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit()); } } } else { size_t i1, i2; - for (i1 = 0; i1 < group1.size(); i1++) { - for (i2 = 0; i2 < group2.size(); i2++) { - cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); - group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); - group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); + for (i1 = 0; i1 < group1->size(); i1++) { + for (i2 = 0; i2 < group2->size(); i2++) { + cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos); + (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit()); + (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit()); } } } @@ -625,16 +627,15 @@ colvar::gyration::gyration(std::string const &conf) : cvc(conf) { function_type = "gyration"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; - parse_group(conf, "atoms", atoms); - atom_groups.push_back(&atoms); + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); + atoms = parse_group(conf, "atoms"); - if (atoms.b_user_defined_fit) { + if (atoms->b_user_defined_fit) { cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); } else { - atoms.b_center = true; - atoms.ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); + atoms->b_center = true; + atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); } x.type(colvarvalue::type_scalar); @@ -644,8 +645,8 @@ colvar::gyration::gyration(std::string const &conf) colvar::gyration::gyration() { function_type = "gyration"; - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); x.type(colvarvalue::type_scalar); } @@ -653,21 +654,21 @@ colvar::gyration::gyration() void colvar::gyration::calc_value() { x.real_value = 0.0; - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { x.real_value += (ai->pos).norm2(); } - x.real_value = std::sqrt(x.real_value / cvm::real(atoms.size())); + x.real_value = std::sqrt(x.real_value / cvm::real(atoms->size())); } void colvar::gyration::calc_gradients() { - cvm::real const drdx = 1.0/(cvm::real(atoms.size()) * x.real_value); - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + cvm::real const drdx = 1.0/(cvm::real(atoms->size()) * x.real_value); + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { ai->grad = drdx * ai->pos; } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(atoms); } @@ -676,12 +677,12 @@ void colvar::gyration::calc_gradients() void colvar::gyration::calc_force_invgrads() { - atoms.read_system_forces(); + atoms->read_system_forces(); cvm::real const dxdr = 1.0/x.real_value; ft.real_value = 0.0; - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { ft.real_value += dxdr * ai->pos * ai->system_force; } } @@ -689,14 +690,14 @@ void colvar::gyration::calc_force_invgrads() void colvar::gyration::calc_Jacobian_derivative() { - jd = x.real_value ? (3.0 * cvm::real(atoms.size()) - 4.0) / x.real_value : 0.0; + jd = x.real_value ? (3.0 * cvm::real(atoms->size()) - 4.0) / x.real_value : 0.0; } void colvar::gyration::apply_force(colvarvalue const &force) { - if (!atoms.noforce) - atoms.apply_colvar_force(force.real_value); + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); } @@ -705,8 +706,6 @@ colvar::inertia::inertia(std::string const &conf) : gyration(conf) { function_type = "inertia"; - b_inverse_gradients = false; - b_Jacobian_derivative = false; x.type(colvarvalue::type_scalar); } @@ -721,7 +720,7 @@ colvar::inertia::inertia() void colvar::inertia::calc_value() { x.real_value = 0.0; - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { x.real_value += (ai->pos).norm2(); } } @@ -729,11 +728,11 @@ void colvar::inertia::calc_value() void colvar::inertia::calc_gradients() { - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { ai->grad = 2.0 * ai->pos; } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(atoms); } @@ -742,8 +741,8 @@ void colvar::inertia::calc_gradients() void colvar::inertia::apply_force(colvarvalue const &force) { - if (!atoms.noforce) - atoms.apply_colvar_force(force.real_value); + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); } @@ -775,7 +774,7 @@ colvar::inertia_z::inertia_z() void colvar::inertia_z::calc_value() { x.real_value = 0.0; - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { cvm::real const iprod = ai->pos * axis; x.real_value += iprod * iprod; } @@ -784,11 +783,11 @@ void colvar::inertia_z::calc_value() void colvar::inertia_z::calc_gradients() { - for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { ai->grad = 2.0 * (ai->pos * axis) * axis; } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(atoms); } @@ -797,8 +796,8 @@ void colvar::inertia_z::calc_gradients() void colvar::inertia_z::apply_force(colvarvalue const &force) { - if (!atoms.noforce) - atoms.apply_colvar_force(force.real_value); + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); } @@ -806,35 +805,35 @@ void colvar::inertia_z::apply_force(colvarvalue const &force) colvar::rmsd::rmsd(std::string const &conf) : cvc(conf) { - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); function_type = "rmsd"; x.type(colvarvalue::type_scalar); - parse_group(conf, "atoms", atoms); - atom_groups.push_back(&atoms); + atoms = parse_group(conf, "atoms"); - if (atoms.b_dummy) { - cvm::error("Error: \"atoms\" cannot be a dummy atom."); + if (atoms->size() == 0) { + cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD."); return; } - if (atoms.ref_pos_group != NULL && b_Jacobian_derivative) { + bool b_Jacobian_derivative = true; + if (atoms->ref_pos_group != NULL && b_Jacobian_derivative) { cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " "Jacobian derivatives of the RMSD will not be calculated.\n"); b_Jacobian_derivative = false; } + if (b_Jacobian_derivative) provide(f_cvc_Jacobian); // the following is a simplified version of the corresponding atom group options; // we need this because the reference coordinates defined inside the atom group // may be used only for fitting, and even more so if ref_pos_group is used if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { cvm::log("Using reference positions from configuration file to calculate the variable.\n"); - if (ref_pos.size() != atoms.size()) { + if (ref_pos.size() != atoms->size()) { cvm::error("Error: the number of reference positions provided ("+ cvm::to_str(ref_pos.size())+ ") does not match the number of atoms of group \"atoms\" ("+ - cvm::to_str(atoms.size())+").\n"); + cvm::to_str(atoms->size())+").\n"); return; } } @@ -861,40 +860,40 @@ colvar::rmsd::rmsd(std::string const &conf) } } else { // if not, rely on existing atom indices for the group - atoms.create_sorted_ids(); - ref_pos.resize(atoms.size()); + atoms->create_sorted_ids(); + ref_pos.resize(atoms->size()); } - cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms.sorted_ids, + cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms->sorted_ids, ref_pos_col, ref_pos_col_value); } } - if (ref_pos.size() != atoms.size()) { + if (ref_pos.size() != atoms->size()) { cvm::error("Error: found " + cvm::to_str(ref_pos.size()) + - " reference positions; expected " + cvm::to_str(atoms.size())); + " reference positions; expected " + cvm::to_str(atoms->size())); return; } - if (atoms.b_user_defined_fit) { + if (atoms->b_user_defined_fit) { cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); } else { // Default: fit everything cvm::log("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.\n"); - atoms.b_center = true; - atoms.b_rotate = true; + atoms->b_center = true; + atoms->b_rotate = true; // default case: reference positions for calculating the rmsd are also those used // for fitting - atoms.ref_pos = ref_pos; - atoms.center_ref_pos(); + atoms->ref_pos = ref_pos; + atoms->center_ref_pos(); cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation " "will not be computed as they cancel out in the gradients."); - atoms.b_fit_gradients = false; + atoms->b_fit_gradients = false; } - if (atoms.b_rotate) { + if (atoms->b_rotate) { // TODO: finer-grained control of this would require exposing a // "request_Jacobian_derivative()" method to the colvar, and the same // from the colvar to biases @@ -902,9 +901,9 @@ colvar::rmsd::rmsd(std::string const &conf) // component - instead it should be decided in a generic way by the atom group // request the calculation of the derivatives of the rotation defined by the atom group - atoms.rot.request_group1_gradients(atoms.size()); + atoms->rot.request_group1_gradients(atoms->size()); // request derivatives of optimal rotation wrt reference coordinates for Jacobian: - atoms.rot.request_group2_gradients(atoms.size()); + atoms->rot.request_group2_gradients(atoms->size()); } } @@ -914,10 +913,10 @@ void colvar::rmsd::calc_value() // rotational-translational fit is handled by the atom group x.real_value = 0.0; - for (size_t ia = 0; ia < atoms.size(); ia++) { - x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2(); + for (size_t ia = 0; ia < atoms->size(); ia++) { + x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2(); } - x.real_value /= cvm::real(atoms.size()); // MSD + x.real_value /= cvm::real(atoms->size()); // MSD x.real_value = std::sqrt(x.real_value); } @@ -925,14 +924,14 @@ void colvar::rmsd::calc_value() void colvar::rmsd::calc_gradients() { cvm::real const drmsddx2 = (x.real_value > 0.0) ? - 0.5 / (x.real_value * cvm::real(atoms.size())) : + 0.5 / (x.real_value * cvm::real(atoms->size())) : 0.0; - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia])); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia])); } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(atoms); } @@ -941,22 +940,22 @@ void colvar::rmsd::calc_gradients() void colvar::rmsd::apply_force(colvarvalue const &force) { - if (!atoms.noforce) - atoms.apply_colvar_force(force.real_value); + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); } void colvar::rmsd::calc_force_invgrads() { - atoms.read_system_forces(); + atoms->read_system_forces(); ft.real_value = 0.0; // Note: gradient square norm is 1/N_atoms - for (size_t ia = 0; ia < atoms.size(); ia++) { - ft.real_value += atoms[ia].grad * atoms[ia].system_force; + for (size_t ia = 0; ia < atoms->size(); ia++) { + ft.real_value += (*atoms)[ia].grad * (*atoms)[ia].system_force; } - ft.real_value *= atoms.size(); + ft.real_value *= atoms->size(); } @@ -965,26 +964,26 @@ void colvar::rmsd::calc_Jacobian_derivative() // divergence of the rotated coordinates (including only derivatives of the rotation matrix) cvm::real divergence = 0.0; - if (atoms.b_rotate) { + if (atoms->b_rotate) { // gradient of the rotation matrix cvm::matrix2d grad_rot_mat(3, 3); // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; - for (size_t ia = 0; ia < atoms.size(); ia++) { + for (size_t ia = 0; ia < atoms->size(); ia++) { // Gradient of optimal quaternion wrt current Cartesian position - cvm::vector1d &dq = atoms.rot.dQ0_1[ia]; + cvm::vector1d &dq = atoms->rot.dQ0_1[ia]; - g11 = 2.0 * (atoms.rot.q)[1]*dq[1]; - g22 = 2.0 * (atoms.rot.q)[2]*dq[2]; - g33 = 2.0 * (atoms.rot.q)[3]*dq[3]; - g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0]; - g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0]; - g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0]; - g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1]; - g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1]; - g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2]; + g11 = 2.0 * (atoms->rot.q)[1]*dq[1]; + g22 = 2.0 * (atoms->rot.q)[2]*dq[2]; + g33 = 2.0 * (atoms->rot.q)[3]*dq[3]; + g01 = (atoms->rot.q)[0]*dq[1] + (atoms->rot.q)[1]*dq[0]; + g02 = (atoms->rot.q)[0]*dq[2] + (atoms->rot.q)[2]*dq[0]; + g03 = (atoms->rot.q)[0]*dq[3] + (atoms->rot.q)[3]*dq[0]; + g12 = (atoms->rot.q)[1]*dq[2] + (atoms->rot.q)[2]*dq[1]; + g13 = (atoms->rot.q)[1]*dq[3] + (atoms->rot.q)[3]*dq[1]; + g23 = (atoms->rot.q)[2]*dq[3] + (atoms->rot.q)[3]*dq[2]; // Gradient of the rotation matrix wrt current Cartesian position grad_rot_mat[0][0] = -2.0 * (g22 + g33); @@ -1009,8 +1008,7 @@ void colvar::rmsd::calc_Jacobian_derivative() } } } - - jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0; + jd.real_value = x.real_value > 0.0 ? (3.0 * atoms->size() - 4.0 - divergence) / x.real_value : 0.0; } @@ -1019,23 +1017,21 @@ void colvar::rmsd::calc_Jacobian_derivative() colvar::eigenvector::eigenvector(std::string const &conf) : cvc(conf) { - b_inverse_gradients = true; - b_Jacobian_derivative = true; + provide(f_cvc_inv_gradient); + provide(f_cvc_Jacobian); function_type = "eigenvector"; x.type(colvarvalue::type_scalar); - parse_group(conf, "atoms", atoms); - atom_groups.push_back(&atoms); - + atoms = parse_group(conf, "atoms"); { bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos); if (b_inline) { cvm::log("Using reference positions from input file.\n"); - if (ref_pos.size() != atoms.size()) { + if (ref_pos.size() != atoms->size()) { cvm::error("Error: reference positions do not " - "match the number of requested atoms.\n"); + "match the number of requested atoms->\n"); return; } } @@ -1060,43 +1056,43 @@ colvar::eigenvector::eigenvector(std::string const &conf) } } else { // if not, use atom indices - atoms.create_sorted_ids(); + atoms->create_sorted_ids(); } - ref_pos.resize(atoms.size()); - cvm::load_coords(file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value); + ref_pos.resize(atoms->size()); + cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value); } } - if (ref_pos.size() != atoms.size()) { + if (ref_pos.size() != atoms->size()) { cvm::error("Error: reference positions were not provided, or do not " - "match the number of requested atoms.\n"); + "match the number of requested atoms->\n"); return; } // save for later the geometric center of the provided positions (may not be the origin) cvm::rvector ref_pos_center(0.0, 0.0, 0.0); - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { ref_pos_center += ref_pos[i]; } - ref_pos_center *= 1.0 / atoms.size(); + ref_pos_center *= 1.0 / atoms->size(); - if (atoms.b_user_defined_fit) { + if (atoms->b_user_defined_fit) { cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); } else { // default: fit everything cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: " "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); - atoms.b_center = true; - atoms.b_rotate = true; - atoms.ref_pos = ref_pos; - atoms.center_ref_pos(); + atoms->b_center = true; + atoms->b_rotate = true; + atoms->ref_pos = ref_pos; + atoms->center_ref_pos(); // request the calculation of the derivatives of the rotation defined by the atom group - atoms.rot.request_group1_gradients(atoms.size()); + atoms->rot.request_group1_gradients(atoms->size()); // request derivatives of optimal rotation wrt reference coordinates for Jacobian: // this is only required for ABF, but we do both groups here for better caching - atoms.rot.request_group2_gradients(atoms.size()); + atoms->rot.request_group2_gradients(atoms->size()); } { @@ -1104,9 +1100,9 @@ colvar::eigenvector::eigenvector(std::string const &conf) // now load the eigenvector if (b_inline) { cvm::log("Using vector components from input file.\n"); - if (eigenvec.size() != atoms.size()) { + if (eigenvec.size() != atoms->size()) { cvm::fatal_error("Error: vector components do not " - "match the number of requested atoms.\n"); + "match the number of requested atoms->\n"); } } @@ -1129,11 +1125,11 @@ colvar::eigenvector::eigenvector(std::string const &conf) } } else { // if not, use atom indices - atoms.create_sorted_ids(); + atoms->create_sorted_ids(); } - eigenvec.resize(atoms.size()); - cvm::load_coords(file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value); + eigenvec.resize(atoms->size()); + cvm::load_coords(file_name.c_str(), eigenvec, atoms->sorted_ids, file_col, file_col_value); } } @@ -1144,10 +1140,10 @@ colvar::eigenvector::eigenvector(std::string const &conf) } cvm::atom_pos eig_center(0.0, 0.0, 0.0); - for (size_t eil = 0; eil < atoms.size(); eil++) { + for (size_t eil = 0; eil < atoms->size(); eil++) { eig_center += eigenvec[eil]; } - eig_center *= 1.0 / atoms.size(); + eig_center *= 1.0 / atoms->size(); cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n"); bool b_difference_vector = false; @@ -1155,33 +1151,33 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (b_difference_vector) { - if (atoms.b_center) { + if (atoms->b_center) { // both sets should be centered on the origin for fitting - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] -= eig_center; ref_pos[i] -= ref_pos_center; } } - if (atoms.b_rotate) { - atoms.rot.calc_optimal_rotation(eigenvec, ref_pos); - for (size_t i = 0; i < atoms.size(); i++) { - eigenvec[i] = atoms.rot.rotate(eigenvec[i]); + if (atoms->b_rotate) { + atoms->rot.calc_optimal_rotation(eigenvec, ref_pos); + for (size_t i = 0; i < atoms->size(); i++) { + eigenvec[i] = atoms->rot.rotate(eigenvec[i]); } } cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n"); - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] -= ref_pos[i]; } - if (atoms.b_center) { + if (atoms->b_center) { // bring back the ref positions to where they were - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { ref_pos[i] += ref_pos_center; } } } else { cvm::log("Centering the provided vector to zero.\n"); - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] -= eig_center; } } @@ -1190,14 +1186,14 @@ colvar::eigenvector::eigenvector(std::string const &conf) // for inverse gradients eigenvec_invnorm2 = 0.0; - for (size_t ein = 0; ein < atoms.size(); ein++) { + for (size_t ein = 0; ein < atoms->size(); ein++) { eigenvec_invnorm2 += eigenvec[ein].norm2(); } eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2; if (b_difference_vector) { cvm::log("\"differenceVector\" is on: normalizing the vector.\n"); - for (size_t i = 0; i < atoms.size(); i++) { + for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] *= eigenvec_invnorm2; } } else { @@ -1209,19 +1205,19 @@ colvar::eigenvector::eigenvector(std::string const &conf) void colvar::eigenvector::calc_value() { x.real_value = 0.0; - for (size_t i = 0; i < atoms.size(); i++) { - x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i]; + for (size_t i = 0; i < atoms->size(); i++) { + x.real_value += ((*atoms)[i].pos - ref_pos[i]) * eigenvec[i]; } } void colvar::eigenvector::calc_gradients() { - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = eigenvec[ia]; + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = eigenvec[ia]; } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging gradients:\n"); debug_gradients(atoms); } @@ -1230,19 +1226,19 @@ void colvar::eigenvector::calc_gradients() void colvar::eigenvector::apply_force(colvarvalue const &force) { - if (!atoms.noforce) - atoms.apply_colvar_force(force.real_value); + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); } void colvar::eigenvector::calc_force_invgrads() { - atoms.read_system_forces(); + atoms->read_system_forces(); ft.real_value = 0.0; - for (size_t ia = 0; ia < atoms.size(); ia++) { - ft.real_value += eigenvec_invnorm2 * atoms[ia].grad * - atoms[ia].system_force; + for (size_t ia = 0; ia < atoms->size(); ia++) { + ft.real_value += eigenvec_invnorm2 * (*atoms)[ia].grad * + (*atoms)[ia].system_force; } } @@ -1251,19 +1247,19 @@ void colvar::eigenvector::calc_Jacobian_derivative() { // gradient of the rotation matrix cvm::matrix2d grad_rot_mat(3, 3); - cvm::quaternion &quat0 = atoms.rot.q; + cvm::quaternion &quat0 = atoms->rot.q; // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; cvm::real sum = 0.0; - for (size_t ia = 0; ia < atoms.size(); ia++) { + for (size_t ia = 0; ia < atoms->size(); ia++) { // Gradient of optimal quaternion wrt current Cartesian position // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t // we can just transpose the derivatives of the direct matrix - cvm::vector1d &dq_1 = atoms.rot.dQ0_1[ia]; + cvm::vector1d &dq_1 = atoms->rot.dQ0_1[ia]; g11 = 2.0 * quat0[1]*dq_1[1]; g22 = 2.0 * quat0[2]*dq_1[2]; @@ -1301,12 +1297,9 @@ void colvar::eigenvector::calc_Jacobian_derivative() colvar::cartesian::cartesian(std::string const &conf) : cvc(conf) { - b_inverse_gradients = false; - b_Jacobian_derivative = false; function_type = "cartesian"; - parse_group(conf, "atoms", atoms); - atom_groups.push_back(&atoms); + atoms = parse_group(conf, "atoms"); bool use_x, use_y, use_z; get_keyval(conf, "useX", use_x, true); @@ -1323,7 +1316,7 @@ colvar::cartesian::cartesian(std::string const &conf) } x.type(colvarvalue::type_vector); - x.vector1d_value.resize(atoms.size() * axes.size()); + x.vector1d_value.resize(atoms->size() * axes.size()); } @@ -1331,9 +1324,9 @@ void colvar::cartesian::calc_value() { size_t const dim = axes.size(); size_t ia, j; - for (ia = 0; ia < atoms.size(); ia++) { + for (ia = 0; ia < atoms->size(); ia++) { for (j = 0; j < dim; j++) { - x.vector1d_value[dim*ia + j] = atoms[ia].pos[axes[j]]; + x.vector1d_value[dim*ia + j] = (*atoms)[ia].pos[axes[j]]; } } } @@ -1351,13 +1344,13 @@ void colvar::cartesian::apply_force(colvarvalue const &force) { size_t const dim = axes.size(); size_t ia, j; - if (!atoms.noforce) { + if (!atoms->noforce) { cvm::rvector f; - for (ia = 0; ia < atoms.size(); ia++) { + for (ia = 0; ia < atoms->size(); ia++) { for (j = 0; j < dim; j++) { f[axes[j]] = force.vector1d_value[dim*ia + j]; } - atoms[ia].apply_force(f); + (*atoms)[ia].apply_force(f); } } } diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 7e0cf75879..505f95e072 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -14,15 +14,14 @@ colvar::orientation::orientation(std::string const &conf) : cvc(conf) { function_type = "orientation"; - parse_group(conf, "atoms", atoms); - atom_groups.push_back(&atoms); + atoms = parse_group(conf, "atoms"); x.type(colvarvalue::type_quaternion); - ref_pos.reserve(atoms.size()); + ref_pos.reserve(atoms->size()); if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) { cvm::log("Using reference positions from input file.\n"); - if (ref_pos.size() != atoms.size()) { + if (ref_pos.size() != atoms->size()) { cvm::fatal_error("Error: reference positions do not " "match the number of requested atoms.\n"); } @@ -42,10 +41,10 @@ colvar::orientation::orientation(std::string const &conf) "if provided, must be non-zero.\n"); } else { // if not, use atom indices - atoms.create_sorted_ids(); + atoms->create_sorted_ids(); } - ref_pos.resize(atoms.size()); - cvm::load_coords(file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value); + ref_pos.resize(atoms->size()); + cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value); } } @@ -71,8 +70,8 @@ colvar::orientation::orientation(std::string const &conf) get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0)); // initialize rot member data - if (!atoms.noforce) { - rot.request_group2_gradients(atoms.size()); + if (!atoms->noforce) { + rot.request_group2_gradients(atoms->size()); } } @@ -88,9 +87,9 @@ colvar::orientation::orientation() void colvar::orientation::calc_value() { - atoms_cog = atoms.center_of_geometry(); + atoms_cog = atoms->center_of_geometry(); - rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); if ((rot.q).inner(ref_quat) >= 0.0) { x.quaternion_value = rot.q; @@ -113,10 +112,10 @@ void colvar::orientation::apply_force(colvarvalue const &force) { cvm::quaternion const &FQ = force.quaternion_value; - if (!atoms.noforce) { - for (size_t ia = 0; ia < atoms.size(); ia++) { + if (!atoms->noforce) { + for (size_t ia = 0; ia < atoms->size(); ia++) { for (size_t i = 0; i < 4; i++) { - atoms[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]); + (*atoms)[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]); } } } @@ -142,9 +141,9 @@ colvar::orientation_angle::orientation_angle() void colvar::orientation_angle::calc_value() { - atoms_cog = atoms.center_of_geometry(); + atoms_cog = atoms->center_of_geometry(); - rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); if ((rot.q).q0 >= 0.0) { x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0); @@ -161,10 +160,10 @@ void colvar::orientation_angle::calc_gradients() ((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) : 0.0 ); - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging orientationAngle component gradients:\n"); debug_gradients(atoms); } @@ -174,8 +173,8 @@ void colvar::orientation_angle::calc_gradients() void colvar::orientation_angle::apply_force(colvarvalue const &force) { cvm::real const &fw = force.real_value; - if (!atoms.noforce) { - atoms.apply_colvar_force(fw); + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); } } @@ -199,8 +198,8 @@ colvar::orientation_proj::orientation_proj() void colvar::orientation_proj::calc_value() { - atoms_cog = atoms.center_of_geometry(); - rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); + atoms_cog = atoms->center_of_geometry(); + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0; } @@ -208,10 +207,10 @@ void colvar::orientation_proj::calc_value() void colvar::orientation_proj::calc_gradients() { cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0; - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging orientationProj component gradients:\n"); debug_gradients(atoms); } @@ -222,8 +221,8 @@ void colvar::orientation_proj::apply_force(colvarvalue const &force) { cvm::real const &fw = force.real_value; - if (!atoms.noforce) { - atoms.apply_colvar_force(fw); + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); } } @@ -255,9 +254,9 @@ colvar::tilt::tilt() void colvar::tilt::calc_value() { - atoms_cog = atoms.center_of_geometry(); + atoms_cog = atoms->center_of_geometry(); - rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); x.real_value = rot.cos_theta(axis); } @@ -267,14 +266,14 @@ void colvar::tilt::calc_gradients() { cvm::quaternion const dxdq = rot.dcos_theta_dq(axis); - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0); for (size_t iq = 0; iq < 4; iq++) { - atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); + (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); } } - if (b_debug_gradients) { + if (is_enabled(f_cvc_debug_gradient)) { cvm::log("Debugging tilt component gradients:\n"); debug_gradients(atoms); } @@ -285,8 +284,8 @@ void colvar::tilt::apply_force(colvarvalue const &force) { cvm::real const &fw = force.real_value; - if (!atoms.noforce) { - atoms.apply_colvar_force(fw); + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); } } @@ -322,9 +321,9 @@ colvar::spin_angle::spin_angle() void colvar::spin_angle::calc_value() { - atoms_cog = atoms.center_of_geometry(); + atoms_cog = atoms->center_of_geometry(); - rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); x.real_value = rot.spin_angle(axis); this->wrap(x); @@ -335,10 +334,10 @@ void colvar::spin_angle::calc_gradients() { cvm::quaternion const dxdq = rot.dspin_angle_dq(axis); - for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0); for (size_t iq = 0; iq < 4; iq++) { - atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); + (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]); } } } @@ -348,7 +347,7 @@ void colvar::spin_angle::apply_force(colvarvalue const &force) { cvm::real const &fw = force.real_value; - if (!atoms.noforce) { - atoms.apply_colvar_force(fw); + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); } } diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index ac192ed0ee..4cf05b0ddb 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -251,15 +251,6 @@ public: size_t i; - for (i = 0; i < cv.size(); i++) { - if (!cv[i]->tasks[colvar::task_lower_boundary] || - !cv[i]->tasks[colvar::task_upper_boundary]) { - cvm::error("Tried to initialize a grid on a " - "variable with undefined boundaries.\n", INPUT_ERROR); - return COLVARS_ERROR; - } - } - if (cvm::debug()) { cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+ " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n"); diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index d5131f6ed5..9991513aca 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -15,6 +15,7 @@ #include "colvarbias_restraint.h" #include "colvarscript.h" + colvarmodule::colvarmodule(colvarproxy *proxy_in) { // pointer to the proxy object @@ -100,6 +101,7 @@ int colvarmodule::read_config_string(std::string const &config_str) return parse_config(conf); } + int colvarmodule::parse_config(std::string &conf) { // parse global options @@ -165,8 +167,9 @@ int colvarmodule::parse_global_params(std::string const &conf) // we are continuing after a reset(): default to true parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append); - parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false, - colvarparse::parse_silent); + parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false); + + parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true); if (use_scripted_forces && !proxy->force_script_defined) { cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR); @@ -218,6 +221,7 @@ int colvarmodule::parse_colvars(std::string const &conf) return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } + bool colvarmodule::check_new_bias(std::string &conf, char const *key) { if (cvm::get_error() || @@ -289,6 +293,12 @@ int colvarmodule::parse_biases(std::string const &conf) cvm::decrease_depth(); } + for (size_t i = 0; i < biases.size(); i++) { + biases[i]->enable(cvm::deps::f_cvb_active); + if (cvm::debug()) + biases[i]->print_state(); + } + if (biases.size() || use_scripted_forces) { cvm::log(cvm::line_marker); cvm::log("Collective variables biases initialized, "+ @@ -448,12 +458,9 @@ int colvarmodule::bias_share(std::string const &bias_name) } -int colvarmodule::calc() { - cvm::real total_bias_energy = 0.0; - cvm::real total_colvar_energy = 0.0; - - std::vector::iterator cvi; - std::vector::iterator bi; +int colvarmodule::calc() +{ + int error_code = COLVARS_OK; if (cvm::debug()) { cvm::log(cvm::line_marker); @@ -461,75 +468,87 @@ int colvarmodule::calc() { cvm::to_str(cvm::step_absolute())+"\n"); } - // calculate collective variables and their gradients - if (cvm::debug()) - cvm::log("Calculating collective variables.\n"); - cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->calc(); - if (cvm::get_error()) { - return COLVARS_ERROR; - } - } - cvm::decrease_depth(); - - // update the biases and communicate their forces to the collective - // variables - if (cvm::debug() && biases.size()) - cvm::log("Updating collective variable biases.\n"); - cvm::increase_depth(); - for (bi = biases.begin(); bi != biases.end(); bi++) { - total_bias_energy += (*bi)->update(); - if (cvm::get_error()) { - return COLVARS_ERROR; - } - } - cvm::decrease_depth(); - - // sum the forces from all biases for each collective variable - if (cvm::debug() && biases.size()) - cvm::log("Collecting forces from all biases.\n"); - cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + error_code |= calc_colvars(); + // set biasing forces to zero before biases are calculated and summed over + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); cvi++) { (*cvi)->reset_bias_force(); } - for (bi = biases.begin(); bi != biases.end(); bi++) { - (*bi)->communicate_forces(); - if (cvm::get_error()) { - return COLVARS_ERROR; - } - } - - // Run user force script, if provided, - // potentially adding scripted forces to the colvars - if (use_scripted_forces) { - int res; - res = proxy->run_force_callback(); - if (res == COLVARS_NOT_IMPLEMENTED) { - cvm::error("Colvar forces scripts are not implemented."); - return COLVARS_ERROR; - } - if (res != COLVARS_OK) { - cvm::error("Error running user colvar forces script"); - return COLVARS_ERROR; - } - } - - cvm::decrease_depth(); + error_code |= calc_biases(); + error_code |= update_colvar_forces(); if (cvm::b_analysis) { - // perform runtime analysis of colvars and biases - if (cvm::debug() && biases.size()) - cvm::log("Perform runtime analyses.\n"); + error_code |= analyze(); + } + + // write trajectory files, if needed + if (cv_traj_freq && cv_traj_name.size()) { + error_code |= write_traj_files(); + } + + // write restart files, if needed + if (restart_out_freq && restart_out_name.size()) { + error_code |= write_restart_files(); + } + + return error_code; +} + + +int colvarmodule::calc_colvars() +{ + if (cvm::debug()) + cvm::log("Calculating collective variables.\n"); + // calculate collective variables and their gradients + + int error_code = COLVARS_OK; + std::vector::iterator cvi; + + // if SMP support is available, split up the work + if (proxy->smp_enabled() == COLVARS_OK) { + + // first, calculate how much work (currently, how many active CVCs) each colvar has + + colvars_smp.resize(0); + colvars_smp_items.resize(0); + + colvars_smp.reserve(colvars.size()); + colvars_smp_items.reserve(colvars.size()); + + // set up a vector containing all components + size_t num_colvar_items = 0; cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->analyze(); - if (cvm::get_error()) { - return COLVARS_ERROR; + + error_code |= (*cvi)->update_cvc_flags(); + + size_t num_items = (*cvi)->num_active_cvcs(); + colvars_smp.reserve(colvars_smp.size() + num_items); + colvars_smp_items.reserve(colvars_smp_items.size() + num_items); + for (size_t icvc = 0; icvc < num_items; icvc++) { + colvars_smp.push_back(*cvi); + colvars_smp_items.push_back(icvc); } + + num_colvar_items += num_items; } - for (bi = biases.begin(); bi != biases.end(); bi++) { - (*bi)->analyze(); + cvm::decrease_depth(); + + // calculate colvar components in parallel + error_code |= proxy->smp_colvars_loop(); + + cvm::increase_depth(); + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + error_code |= (*cvi)->collect_cvc_data(); + } + cvm::decrease_depth(); + + } else { + + // calculate colvars one at a time + cvm::increase_depth(); + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + error_code |= (*cvi)->calc(); if (cvm::get_error()) { return COLVARS_ERROR; } @@ -537,6 +556,81 @@ int colvarmodule::calc() { cvm::decrease_depth(); } + return error_code | cvm::get_error(); +} + + +int colvarmodule::calc_biases() +{ + // update the biases and communicate their forces to the collective + // variables + if (cvm::debug() && biases.size()) + cvm::log("Updating collective variable biases.\n"); + + std::vector::iterator bi; + int error_code = COLVARS_OK; + + // if SMP support is available, split up the work + if (proxy->smp_enabled() == COLVARS_OK) { + + if (use_scripted_forces && !scripting_after_biases) { + // calculate biases and scripted forces in parallel + error_code |= proxy->smp_biases_script_loop(); + } else { + // calculate biases in parallel + error_code |= proxy->smp_biases_loop(); + } + + } else { + + if (use_scripted_forces && !scripting_after_biases) { + error_code |= calc_scripted_forces(); + } + + cvm::increase_depth(); + for (bi = biases.begin(); bi != biases.end(); bi++) { + error_code |= (*bi)->update(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } + } + cvm::decrease_depth(); + } + + cvm::real total_bias_energy = 0.0; + for (bi = biases.begin(); bi != biases.end(); bi++) { + total_bias_energy += (*bi)->get_energy(); + } + + proxy->add_energy(total_bias_energy); + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + + +int colvarmodule::update_colvar_forces() +{ + int error_code = COLVARS_OK; + + std::vector::iterator cvi; + std::vector::iterator bi; + + // sum the forces from all biases for each collective variable + if (cvm::debug() && biases.size()) + cvm::log("Collecting forces from all biases.\n"); + cvm::increase_depth(); + for (bi = biases.begin(); bi != biases.end(); bi++) { + (*bi)->communicate_forces(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } + } + cvm::decrease_depth(); + + if (use_scripted_forces && scripting_after_biases) { + error_code |= calc_scripted_forces(); + } + + cvm::real total_colvar_energy = 0.0; // sum up the forces for each colvar, including wall forces // and integrate any internal // equation of motion (extended system) @@ -545,13 +639,13 @@ int colvarmodule::calc() { "of colvars (if they have any).\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - total_colvar_energy += (*cvi)->update(); + total_colvar_energy += (*cvi)->update_forces_energy(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); - proxy->add_energy(total_bias_energy + total_colvar_energy); + proxy->add_energy(total_colvar_energy); // make collective variables communicate their forces to their // coupled degrees of freedom (i.e. atoms) @@ -559,7 +653,7 @@ int colvarmodule::calc() { cvm::log("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - if ((*cvi)->tasks[colvar::task_gradients]) { + if ((*cvi)->is_enabled(cvm::deps::f_cv_gradient)) { (*cvi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; @@ -568,46 +662,69 @@ int colvarmodule::calc() { } cvm::decrease_depth(); - // write restart file, if needed - if (restart_out_freq && restart_out_name.size()) { - if ( (cvm::step_relative() > 0) && - ((cvm::step_absolute() % restart_out_freq) == 0) ) { - cvm::log("Writing the state file \""+ - restart_out_name+"\".\n"); - proxy->backup_file(restart_out_name.c_str()); - restart_out_os.open(restart_out_name.c_str()); - if (!restart_out_os.is_open() || !write_restart(restart_out_os)) - cvm::error("Error: in writing restart file.\n"); - restart_out_os.close(); - } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + + +int colvarmodule::calc_scripted_forces() +{ + // Run user force script, if provided, + // potentially adding scripted forces to the colvars + int res; + res = proxy->run_force_callback(); + if (res == COLVARS_NOT_IMPLEMENTED) { + cvm::error("Colvar forces scripts are not implemented."); + return COLVARS_NOT_IMPLEMENTED; + } + if (res != COLVARS_OK) { + cvm::error("Error running user colvar forces script"); + return COLVARS_ERROR; + } + return COLVARS_OK; +} + + +int colvarmodule::write_restart_files() +{ + if ( (cvm::step_relative() > 0) && + ((cvm::step_absolute() % restart_out_freq) == 0) ) { + cvm::log("Writing the state file \""+ + restart_out_name+"\".\n"); + proxy->backup_file(restart_out_name.c_str()); + restart_out_os.open(restart_out_name.c_str()); + if (!restart_out_os.is_open() || !write_restart(restart_out_os)) + cvm::error("Error: in writing restart file.\n"); + restart_out_os.close(); } - // write trajectory file, if needed - if (cv_traj_freq && cv_traj_name.size()) { + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} - if (!cv_traj_os.is_open()) { - open_traj_file(cv_traj_name); - } - // write labels in the traj file every 1000 lines and at first timestep - if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { - write_traj_label(cv_traj_os); - } +int colvarmodule::write_traj_files() +{ + if (!cv_traj_os.is_open()) { + open_traj_file(cv_traj_name); + } - if ((cvm::step_absolute() % cv_traj_freq) == 0) { - write_traj(cv_traj_os); - } + // write labels in the traj file every 1000 lines and at first timestep + if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { + write_traj_label(cv_traj_os); + } - if (restart_out_freq && cv_traj_os.is_open()) { - // flush the trajectory file if we are at the restart frequency - if ( (cvm::step_relative() > 0) && - ((cvm::step_absolute() % restart_out_freq) == 0) ) { - cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+ - cv_traj_name+"\".\n"); - cv_traj_os.flush(); - } + if ((cvm::step_absolute() % cv_traj_freq) == 0) { + write_traj(cv_traj_os); + } + + if (restart_out_freq && cv_traj_os.is_open()) { + // flush the trajectory file if we are at the restart frequency + if ( (cvm::step_relative() > 0) && + ((cvm::step_absolute() % restart_out_freq) == 0) ) { + cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+ + cv_traj_name+"\".\n"); + cv_traj_os.flush(); } - } // end if (cv_traj_freq) + } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -643,6 +760,7 @@ int colvarmodule::analyze() return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } + int colvarmodule::setup() { // loop over all components of all colvars to reset masses of all groups @@ -653,25 +771,23 @@ int colvarmodule::setup() return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } + colvarmodule::~colvarmodule() { - reset(); - delete parse; - proxy = NULL; + if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) || + (proxy->smp_thread_id() == 0)) { + reset(); + delete parse; + proxy = NULL; + } } + int colvarmodule::reset() { parse->reset(); cvm::log("Resetting the Collective Variables Module.\n"); - // Iterate backwards because we are deleting the elements as we go - for (std::vector::reverse_iterator cvi = colvars.rbegin(); - cvi != colvars.rend(); - cvi++) { - delete *cvi; // the colvar destructor updates the colvars array - } - colvars.clear(); // Iterate backwards because we are deleting the elements as we go for (std::vector::reverse_iterator bi = biases.rbegin(); @@ -681,6 +797,14 @@ int colvarmodule::reset() } biases.clear(); + // Iterate backwards because we are deleting the elements as we go + for (std::vector::reverse_iterator cvi = colvars.rbegin(); + cvi != colvars.rend(); + cvi++) { + delete *cvi; // the colvar destructor updates the colvars array + } + colvars.clear(); + index_groups.clear(); index_group_names.clear(); @@ -983,7 +1107,7 @@ int colvarmodule::open_traj_file(std::string const &file_name) int colvarmodule::close_traj_file() { - if (cv_traj_os.good()) { + if (cv_traj_os.is_open()) { cv_traj_os.close(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -1048,28 +1172,71 @@ std::ostream & colvarmodule::write_traj(std::ostream &os) void cvm::log(std::string const &message) { - if (depth > 0) - proxy->log((std::string(2*depth, ' '))+message); + size_t const d = depth(); + if (d > 0) + proxy->log((std::string(2*d, ' '))+message); else proxy->log(message); } + void cvm::increase_depth() { - depth++; + (depth())++; } + void cvm::decrease_depth() { - if (depth) depth--; + if (depth() > 0) { + (depth())--; + } } + +size_t & cvm::depth() +{ + // NOTE: do not call log() or error() here, to avoid recursion + size_t const nt = proxy->smp_num_threads(); + if (proxy->smp_enabled() == COLVARS_OK) { + if (depth_v.size() != nt) { + // update array of depths + proxy->smp_lock(); + if (depth_v.size() > 0) { depth_s = depth_v[0]; } + depth_v.clear(); + depth_v.assign(nt, depth_s); + proxy->smp_unlock(); + } + return depth_v[proxy->smp_thread_id()]; + } + return depth_s; +} + + +void colvarmodule::set_error_bits(int code) +{ + proxy->smp_lock(); + errorCode |= code; + errorCode |= COLVARS_ERROR; + proxy->smp_unlock(); +} + + +void colvarmodule::clear_error() +{ + proxy->smp_lock(); + errorCode = COLVARS_OK; + proxy->smp_unlock(); +} + + void cvm::error(std::string const &message, int code) { set_error_bits(code); proxy->error(message); } + void cvm::fatal_error(std::string const &message) { // TODO once all non-fatal errors have been set to be handled by error(), @@ -1078,6 +1245,7 @@ void cvm::fatal_error(std::string const &message) proxy->fatal_error(message); } + void cvm::exit(std::string const &message) { proxy->exit(message); @@ -1232,11 +1400,13 @@ long colvarmodule::it = 0; long colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; -size_t colvarmodule::depth = 0; +size_t colvarmodule::depth_s = 0; +std::vector colvarmodule::depth_v(0); bool colvarmodule::b_analysis = false; std::list colvarmodule::index_group_names; std::list > colvarmodule::index_groups; bool colvarmodule::use_scripted_forces = false; +bool colvarmodule::scripting_after_biases = true; // file name prefixes std::string colvarmodule::output_prefix = ""; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 32ec8fe4ad..eac071805c 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2016-03-08" +#define COLVARS_VERSION "2016-04-14" #endif #ifndef COLVARS_DEBUG @@ -74,6 +74,9 @@ private: public: + /// Base class to handle mutual dependencies of most objects + class deps; + friend class colvarproxy; // TODO colvarscript should be unaware of colvarmodule's internals friend class colvarscript; @@ -106,20 +109,21 @@ public: /// Module-wide error state /// see constants at the top of this file +protected: + static int errorCode; - static inline void set_error_bits(int code) - { - errorCode |= code; - errorCode |= COLVARS_ERROR; - } + +public: + + static void set_error_bits(int code); + static inline int get_error() { return errorCode; } - static inline void clear_error() - { - errorCode = 0; - } + + static void clear_error(); + /// Current step number static long it; @@ -167,6 +171,12 @@ public: } */ + /// Collective variables to be calculated on different threads; + /// colvars with multple items (e.g. multiple active CVCs) are duplicated + std::vector colvars_smp; + /// Indexes of the items to calculate for each colvar + std::vector colvars_smp_items; + /// Array of collective variable biases static std::vector biases; /// \brief Number of ABF biases initialized (in normal conditions @@ -269,6 +279,10 @@ public: /// Write explanatory labels in the trajectory file std::ostream & write_traj_label(std::ostream &os); + /// Write all trajectory files + int write_traj_files(); + /// Write all restart files + int write_restart_files(); /// Write all FINAL output files int write_output_files(); /// Backup a file before writing it @@ -300,11 +314,21 @@ public: //// Share among replicas. int bias_share(std::string const &bias_name); - /// Calculate collective variables and biases + /// Main worker function int calc(); + /// Calculate collective variables + int calc_colvars(); + + /// Calculate biases + int calc_biases(); + + /// Integrate bias and restraint forces, send colvar forces to atoms + int update_colvar_forces(); + /// Perform analysis int analyze(); + /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) int read_traj(char const *traj_filename, @@ -480,18 +504,18 @@ protected: /// Output restart file colvarmodule::ofstream restart_out_os; - /// \brief Counter for the current depth in the object hierarchy (useg e.g. in output - static size_t depth; +protected: - /// Use scripted colvars forces? - static bool use_scripted_forces; + /// Counter for the current depth in the object hierarchy (useg e.g. in output) + static size_t depth_s; + + /// Thread-specific depth + static std::vector depth_v; public: - /// \brief Pointer to the proxy object, used to retrieve atomic data - /// from the hosting program; it is static in order to be accessible - /// from static functions in the colvarmodule class - static colvarproxy *proxy; + /// Get the current object depth in the hierarchy + static size_t & depth(); /// Increase the depth (number of indentations in the output) static void increase_depth(); @@ -499,7 +523,24 @@ public: /// Decrease the depth (number of indentations in the output) static void decrease_depth(); - static inline bool scripted_forces() { return use_scripted_forces; } + static inline bool scripted_forces() + { + return use_scripted_forces; + } + + /// Use scripted colvars forces? + static bool use_scripted_forces; + + /// Wait for all biases before calculating scripted forces? + static bool scripting_after_biases; + + /// Calculate the energy and forces of scripted biases + int calc_scripted_forces(); + + /// \brief Pointer to the proxy object, used to retrieve atomic data + /// from the hosting program; it is static in order to be accessible + /// from static functions in the colvarmodule class + static colvarproxy *proxy; }; diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index d18ec3715a..585a40cdfc 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -113,6 +113,62 @@ protected: public: + // ***************** SHARED-MEMORY PARALLELIZATION ***************** + + /// Whether or not threaded parallelization is available + virtual int smp_enabled() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Distribute calculation of colvars (and their components) across threads + virtual int smp_colvars_loop() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Distribute calculation of biases across threads + virtual int smp_biases_loop() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread + virtual int smp_biases_script_loop() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Index of this thread + virtual int smp_thread_id() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Number of threads sharing this address space + virtual int smp_num_threads() + { + return COLVARS_NOT_IMPLEMENTED; + } + + /// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing + virtual int smp_lock() + { + return COLVARS_OK; + } + + /// Attempt to lock the proxy's shared data + virtual int smp_trylock() + { + return COLVARS_OK; + } + + /// Release the lock + virtual int smp_unlock() + { + return COLVARS_OK; + } + // **************** MULTIPLE REPLICAS COMMUNICATION **************** // Replica exchange commands: @@ -466,9 +522,9 @@ protected: public: /// \brief Whether this proxy implementation has capability for scalable groups - virtual bool has_scalable_groups() const + virtual int scalable_group_coms() { - return false; + return COLVARS_NOT_IMPLEMENTED; } /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet @@ -476,6 +532,7 @@ public: inline int add_atom_group_slot(int atom_group_id) { atom_groups_ids.push_back(atom_group_id); + atom_groups_ncopies.push_back(1); atom_groups_masses.push_back(1.0); atom_groups_charges.push_back(0.0); atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0)); @@ -496,18 +553,18 @@ public: /// \brief Used by the atom_group class destructor virtual void clear_atom_group(int index) { + if (cvm::debug()) { + log("Trying to remove/disable atom group number "+cvm::to_str(index)+"\n"); + } + if (((size_t) index) >= atom_groups_ids.size()) { cvm::error("Error: trying to disable an atom group that was not previously requested.\n", INPUT_ERROR); } - atom_groups_ids.erase(atom_groups_ids.begin()+index); - atom_groups_masses.erase(atom_groups_masses.begin()+index); - atom_groups_charges.erase(atom_groups_charges.begin()+index); - atom_groups_coms.erase(atom_groups_coms.begin()+index); - atom_groups_total_forces.erase(atom_groups_total_forces.begin()+index); - atom_groups_applied_forces.erase(atom_groups_applied_forces.begin()+index); - atom_groups_new_colvar_forces.erase(atom_groups_new_colvar_forces.begin()+index); + if (atom_groups_ncopies[index] > 0) { + atom_groups_ncopies[index] -= 1; + } } /// Get the numeric ID of the given atom group (for the MD program) diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index efc6731ef7..bd1c8878fb 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -221,7 +221,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { if (subcmd == "update") { cv->calc(); - cv->update(); + cv->update_forces_energy(); result = (cv->value()).to_simple_string(); return COLVARS_OK; } @@ -243,6 +243,16 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { return COLVARS_OK; } + if (subcmd == "getappliedforce") { + result = (cv->bias_force()).to_simple_string(); + return COLVARS_OK; + } + + if (subcmd == "getsystemforce") { + result = (cv->system_force()).to_simple_string(); + return COLVARS_OK; + } + if (subcmd == "addforce") { if (argc < 4) { result = "addforce: missing parameter: force value\n" + help_string();