diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 5b776defca36759dadaaabebcc6df409000d52cd..37201275fe4b4efa0da77afbb72e6c55a127c042 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index d60a9e90f774179fb1fa2bcf6d8b61452786048f..e8c7e88324a9e88299c4c7c1834ce4c7b21e8e90 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1,5 +1,5 @@
-// -*- c++ -*-
+// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
@@ -7,6 +7,7 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
@@ -23,20 +24,31 @@ bool compare(colvar::cvc *i, colvar::cvc *j) {
}
-colvar::colvar(std::string const &conf)
- : colvarparse(conf)
+colvar::colvar()
+{
+ // Initialize static array once and for all
+ init_cv_requires();
+}
+
+
+int colvar::init(std::string const &conf)
{
cvm::log("Initializing a new collective variable.\n");
+ colvarparse::init(conf);
+
int error_code = COLVARS_OK;
+ colvarmodule *cv = cvm::main();
+
get_keyval(conf, "name", this->name,
- (std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));
+ (std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));
- if (cvm::colvar_by_name(this->name) != NULL) {
+ if ((cvm::colvar_by_name(this->name) != NULL) &&
+ (cvm::colvar_by_name(this->name) != this)) {
cvm::error("Error: this colvar cannot have the same name, \""+this->name+
"\", as another colvar.\n",
INPUT_ERROR);
- return;
+ return INPUT_ERROR;
}
// Initialize dependency members
@@ -44,14 +56,13 @@ colvar::colvar(std::string const &conf)
this->description = "colvar " + this->name;
- // Initialize static array once and for all
- init_cv_requires();
-
kinetic_energy = 0.0;
potential_energy = 0.0;
error_code |= init_components(conf);
- if (error_code != COLVARS_OK) return;
+ if (error_code != COLVARS_OK) {
+ return cvm::get_error();
+ }
size_t i;
@@ -59,8 +70,6 @@ colvar::colvar(std::string const &conf)
if (get_keyval(conf, "scriptedFunction", scripted_function,
"", colvarparse::parse_silent)) {
- // Make feature available only on user request
- provide(f_cv_scripted);
enable(f_cv_scripted);
cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
@@ -76,8 +85,8 @@ colvar::colvar(std::string const &conf)
}
}
if (x.type() == colvarvalue::type_notset) {
- cvm::error("Could not parse scripted colvar type.");
- return;
+ cvm::error("Could not parse scripted colvar type.", INPUT_ERROR);
+ return INPUT_ERROR;
}
cvm::log(std::string("Expecting colvar value of type ")
@@ -86,8 +95,9 @@ colvar::colvar(std::string const &conf)
if (x.type() == colvarvalue::type_vector) {
int size;
if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
- cvm::error("Error: no size specified for vector scripted function.");
- return;
+ cvm::error("Error: no size specified for vector scripted function.",
+ INPUT_ERROR);
+ return INPUT_ERROR;
}
x.vector1d_value.resize(size);
}
@@ -154,7 +164,7 @@ colvar::colvar(std::string const &conf)
}
}
}
- feature_states[f_cv_linear]->enabled = lin;
+ set_enabled(f_cv_linear, lin);
}
// Colvar is homogeneous if:
@@ -168,7 +178,7 @@ colvar::colvar(std::string const &conf)
homogeneous = false;
}
}
- feature_states[f_cv_homogeneous]->enabled = homogeneous;
+ set_enabled(f_cv_homogeneous, homogeneous);
}
// Colvar is deemed periodic if:
@@ -188,7 +198,7 @@ colvar::colvar(std::string const &conf)
"Make sure that you know what you are doing!");
}
}
- feature_states[f_cv_periodic]->enabled = b_periodic;
+ set_enabled(f_cv_periodic, b_periodic);
}
// check that cvcs are compatible
@@ -202,7 +212,7 @@ colvar::colvar(std::string const &conf)
"by using components of different types. "
"You must use the same type in order to "
"sum them together.\n", INPUT_ERROR);
- return;
+ return INPUT_ERROR;
}
}
@@ -215,44 +225,110 @@ colvar::colvar(std::string const &conf)
f.type(value());
f_accumulated.type(value());
+ x_old.type(value());
+ v_fdiff.type(value());
+ v_reported.type(value());
+ fj.type(value());
+ ft.type(value());
+ ft_reported.type(value());
+ f_old.type(value());
+ f_old.reset();
+
+ x_restart.type(value());
+ after_restart = false;
+
reset_bias_force();
+ // TODO use here information from the CVCs' own natural boundaries
+ error_code |= init_grid_parameters(conf);
+
+ get_keyval(conf, "timeStepFactor", time_step_factor, 1);
+
+ error_code |= init_extended_Lagrangian(conf);
+ error_code |= init_output_flags(conf);
+
+ // Start in active state by default
+ enable(f_cv_active);
+ // Make sure dependency side-effects are correct
+ refresh_deps();
+
+ if (cvm::b_analysis)
+ parse_analysis(conf);
+
+ if (cvm::debug())
+ cvm::log("Done initializing collective variable \""+this->name+"\".\n");
+
+ return error_code;
+}
+
+
+int colvar::init_grid_parameters(std::string const &conf)
+{
+ colvarmodule *cv = cvm::main();
+
get_keyval(conf, "width", width, 1.0);
if (width <= 0.0) {
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
- return;
+ return INPUT_ERROR;
}
lower_boundary.type(value());
- lower_wall.type(value());
upper_boundary.type(value());
upper_wall.type(value());
- feature_states[f_cv_scalar]->enabled = (value().type() == colvarvalue::type_scalar);
+ set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
if (is_enabled(f_cv_scalar)) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
- provide(f_cv_lower_boundary);
enable(f_cv_lower_boundary);
}
+ std::string lw_conf, uw_conf;
- get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0);
- if (lower_wall_k > 0.0) {
+ if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) {
+ cvm::log("Warning: lowerWallConstant and lowerWall are deprecated, "
+ "please define a harmonicWalls bias instead.\n");
+ lower_wall.type(value());
get_keyval(conf, "lowerWall", lower_wall, lower_boundary);
- enable(f_cv_lower_wall);
+ lw_conf = std::string("\n\
+ lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
+ lowerWalls "+cvm::to_str(lower_wall)+"\n");
}
if (get_keyval(conf, "upperBoundary", upper_boundary, 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) {
+ if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) {
+ cvm::log("Warning: upperWallConstant and upperWall are deprecated, "
+ "please define a harmonicWalls bias instead.\n");
+ upper_wall.type(value());
get_keyval(conf, "upperWall", upper_wall, upper_boundary);
- enable(f_cv_upper_wall);
+ uw_conf = std::string("\n\
+ upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
+ upperWalls "+cvm::to_str(upper_wall)+"\n");
+ }
+
+ if (lw_conf.size() && uw_conf.size()) {
+ if (lower_wall >= upper_wall) {
+ cvm::error("Error: the upper wall, "+
+ cvm::to_str(upper_wall)+
+ ", is not higher than the lower wall, "+
+ cvm::to_str(lower_wall)+".\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ }
+
+ if (lw_conf.size() || uw_conf.size()) {
+ cvm::log("Generating a new harmonicWalls bias for compatibility purposes.\n");
+ std::string const walls_conf("\n\
+harmonicWalls {\n\
+ name "+this->name+"w\n\
+ colvars "+this->name+"\n"+lw_conf+uw_conf+
+ "}\n");
+ cv->append_new_config(walls_conf);
}
}
@@ -271,16 +347,7 @@ colvar::colvar(std::string const &conf)
", is not higher than the lower boundary, "+
cvm::to_str(lower_boundary)+".\n",
INPUT_ERROR);
- }
- }
-
- 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)+
- ", is not higher than the lower wall, "+
- cvm::to_str(lower_wall)+".\n",
- INPUT_ERROR);
+ return INPUT_ERROR;
}
}
@@ -289,83 +356,90 @@ colvar::colvar(std::string const &conf)
cvm::error("Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.\n",
INPUT_ERROR);
+ return INPUT_ERROR;
}
if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
cvm::error("Error: inconsistent configuration "
"(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).\n",
INPUT_ERROR);
+ return INPUT_ERROR;
}
- get_keyval(conf, "timeStepFactor", time_step_factor, 1);
+ return COLVARS_OK;
+}
- {
- bool b_extended_Lagrangian;
- get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
- if (b_extended_Lagrangian) {
- cvm::real temp, tolerance, period;
+int colvar::init_extended_Lagrangian(std::string const &conf)
+{
+ bool b_extended_Lagrangian;
+ get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
- cvm::log("Enabling the extended Lagrangian term for colvar \""+
- this->name+"\".\n");
+ if (b_extended_Lagrangian) {
+ cvm::real temp, tolerance, period;
- // Make feature available only on user request
- provide(f_cv_extended_Lagrangian);
- enable(f_cv_extended_Lagrangian);
- provide(f_cv_Langevin);
-
- // The extended mass will apply forces
- enable(f_cv_gradient);
-
- xr.type(value());
- vr.type(value());
- fr.type(value());
-
- const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
- if (temp <= 0.0) {
- if (found)
- cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
- else
- cvm::error("Error: a positive temperature must be provided, either "
- "by enabling a thermostat, or through \"extendedTemp\".\n",
- INPUT_ERROR);
- }
+ cvm::log("Enabling the extended Lagrangian term for colvar \""+
+ this->name+"\".\n");
- get_keyval(conf, "extendedFluctuation", tolerance);
- if (tolerance <= 0.0) {
- cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
- }
- ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
- cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
+ enable(f_cv_extended_Lagrangian);
- get_keyval(conf, "extendedTimeConstant", period, 200.0);
- if (period <= 0.0) {
- cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
- }
- ext_mass = (cvm::boltzmann() * temp * period * period)
- / (4.0 * PI * PI * tolerance * tolerance);
- cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");
-
- {
- bool b_output_energy;
- get_keyval(conf, "outputEnergy", b_output_energy, false);
- if (b_output_energy) {
- enable(f_cv_output_energy);
- }
- }
+ xr.type(value());
+ vr.type(value());
+ fr.type(value());
- get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
- if (ext_gamma < 0.0) {
- cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
- }
- if (ext_gamma != 0.0) {
- 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());
+ const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
+ if (temp <= 0.0) {
+ if (found)
+ cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
+ else
+ cvm::error("Error: a positive temperature must be provided, either "
+ "by enabling a thermostat, or through \"extendedTemp\".\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+
+ get_keyval(conf, "extendedFluctuation", tolerance);
+ if (tolerance <= 0.0) {
+ cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
+ cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
+
+ get_keyval(conf, "extendedTimeConstant", period, 200.0);
+ if (period <= 0.0) {
+ cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
+ }
+ ext_mass = (cvm::boltzmann() * temp * period * period)
+ / (4.0 * PI * PI * tolerance * tolerance);
+ cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");
+
+ {
+ bool b_output_energy;
+ get_keyval(conf, "outputEnergy", b_output_energy, false);
+ if (b_output_energy) {
+ enable(f_cv_output_energy);
}
}
+
+ get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
+ if (ext_gamma < 0.0) {
+ cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ if (ext_gamma != 0.0) {
+ 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());
+ }
}
+ return COLVARS_OK;
+}
+
+
+int colvar::init_output_flags(std::string const &conf)
+{
{
bool b_output_value;
get_keyval(conf, "outputValue", b_output_value, true);
@@ -387,7 +461,7 @@ colvar::colvar(std::string const &conf)
if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
"The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
- return;
+ return INPUT_ERROR;
}
}
@@ -395,28 +469,7 @@ colvar::colvar(std::string const &conf)
get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
- // 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());
- f_old.type(value());
- f_old.reset();
-
- x_restart.type(value());
- after_restart = false;
-
- if (cvm::b_analysis)
- parse_analysis(conf);
-
- if (cvm::debug())
- cvm::log("Done initializing collective variable \""+this->name+"\".\n");
+ return COLVARS_OK;
}
@@ -637,7 +690,7 @@ int colvar::parse_analysis(std::string const &conf)
std::string runave_outfile;
get_keyval(conf, "runAveOutputFile", runave_outfile,
- std::string(cvm::output_prefix+"."+
+ std::string(cvm::output_prefix()+"."+
this->name+".runave.traj"));
size_t const this_cv_width = x.output_width(cvm::cv_width);
@@ -693,7 +746,7 @@ int colvar::parse_analysis(std::string const &conf)
get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
get_keyval(conf, "corrFuncOutputFile", acf_outfile,
- std::string(cvm::output_prefix+"."+this->name+
+ std::string(cvm::output_prefix()+"."+this->name+
".corrfunc.dat"));
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@@ -730,11 +783,12 @@ colvar::~colvar()
}
// remove reference to this colvar from the CVM
- for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
- cvi != cvm::colvars.end();
+ colvarmodule *cv = cvm::main();
+ for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
+ cvi != cv->variables()->end();
++cvi) {
if ( *cvi == this) {
- cvm::colvars.erase(cvi);
+ cv->variables()->erase(cvi);
break;
}
}
@@ -892,7 +946,6 @@ int colvar::collect_cvc_values()
cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
if (after_restart) {
- after_restart = false;
if (cvm::proxy->simulation_running()) {
cvm::real const jump2 = dist2(x, x_restart) / (width*width);
if (jump2 > 0.25) {
@@ -1122,8 +1175,7 @@ int colvar::calc_colvar_properties()
// initialize the restraint center in the first step to the value
// just calculated from the cvcs
- // TODO: put it in the restart information
- if (cvm::step_relative() == 0) {
+ if (cvm::step_relative() == 0 && !after_restart) {
xr = x;
vr.reset(); // (already 0; added for clarity)
}
@@ -1148,6 +1200,8 @@ int colvar::calc_colvar_properties()
ft_reported = ft;
}
+ // At the end of the first update after a restart, we can reset the flag
+ after_restart = false;
return COLVARS_OK;
}
@@ -1173,51 +1227,17 @@ cvm::real colvar::update_forces_energy()
f -= fj;
}
- // Wall force
- colvarvalue fw(x);
- fw.reset();
-
- if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
-
- if (cvm::debug())
- cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
-
- // For a periodic colvar, both walls may be applicable at the same time
- // in which case we pick the closer one
- 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);
- if (grad < 0.0) {
- fw = -0.5 * lower_wall_k * grad;
- if (cvm::debug())
- cvm::log("Applying a lower wall force("+
- cvm::to_str(fw)+") to \""+this->name+"\".\n");
- }
-
- } else {
-
- cvm::real const grad = this->dist2_lgrad(x, upper_wall);
- if (grad > 0.0) {
- fw = -0.5 * upper_wall_k * grad;
- if (cvm::debug())
- cvm::log("Applying an upper wall force("+
- cvm::to_str(fw)+") to \""+this->name+"\".\n");
- }
- }
- }
-
// At this point f is the force f from external biases that will be applied to the
// extended variable if there is one
if (is_enabled(f_cv_extended_Lagrangian)) {
if (cvm::debug()) {
- cvm::log("Updating extended-Lagrangian degrees of freedom.\n");
+ cvm::log("Updating extended-Lagrangian degree of freedom.\n");
}
cvm::real dt = cvm::dt();
- colvarvalue f_ext(fr.type());
+ colvarvalue f_ext(fr.type()); // force acting on the extended variable
f_ext.reset();
// the total force is applied to the fictitious mass, while the
@@ -1226,7 +1246,7 @@ cvm::real colvar::update_forces_energy()
// f_ext: total force on extended variable (including harmonic spring)
// f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates
- // ie. spring force + wall force
+ // ie. spring force (fb_actual will be added just below)
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
@@ -1259,8 +1279,6 @@ cvm::real colvar::update_forces_energy()
if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
}
- // TODO remove the wall force
- f += fw;
// Now adding the force on the actual colvar (for those biases who
// bypass the extended Lagrangian mass)
f += fb_actual;
@@ -1286,8 +1304,10 @@ cvm::real colvar::update_forces_energy()
void colvar::communicate_forces()
{
size_t i;
- if (cvm::debug())
+ if (cvm::debug()) {
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
+ cvm::log("Force to be applied: " + cvm::to_str(f_accumulated) + "\n");
+ }
if (is_enabled(f_cv_scripted)) {
std::vector<cvm::matrix2d<cvm::real> > func_grads;
@@ -1364,7 +1384,7 @@ int colvar::update_cvc_flags()
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];
+ cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
if (cvcs[i]->is_enabled()) {
n_active_cvcs++;
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 2cf3d2dac512b7492452940393067e8d6f7088ef..0cbda450b8af878ba611f1689cc8d85628e8e0e5 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -227,11 +227,23 @@ public:
/// Constructor
- colvar(std::string const &conf);
+ colvar();
+
+ /// Main init function
+ int init(std::string const &conf);
/// Parse the CVC configuration and allocate their data
int init_components(std::string const &conf);
+ /// Init defaults for grid options
+ int init_grid_parameters(std::string const &conf);
+
+ /// Init extended Lagrangian parameters
+ int init_extended_Lagrangian(std::string const &conf);
+
+ /// Init output flags
+ int init_output_flags(std::string const &conf);
+
private:
/// Parse the CVC configuration for all components of a certain type
template<typename def_class_name> int init_components_type(std::string const &conf,
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 48c16e887af95c9e623aa3b8b525b54594e10ca2..32cfadf3b6c84aae0af422fb63f6bcad0d870158 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -98,7 +98,7 @@ cvm::atom_group::atom_group()
cvm::atom_group::~atom_group()
{
- if (is_enabled(f_ag_scalable)) {
+ if (is_enabled(f_ag_scalable) && !b_dummy) {
(cvm::proxy)->clear_atom_group(index);
index = -1;
}
@@ -418,7 +418,7 @@ int cvm::atom_group::parse(std::string const &conf)
// We need to know the fitting options to decide whether the group is scalable
parse_error |= parse_fitting_options(group_conf);
- if (is_available(f_ag_scalable_com) && !b_rotate) {
+ if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) {
enable(f_ag_scalable_com);
enable(f_ag_scalable);
}
@@ -500,14 +500,16 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
int cvm::atom_group::add_index_group(std::string const &index_group_name)
{
- std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
- std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
- for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) {
+ colvarmodule *cv = cvm::main();
+
+ std::list<std::string>::iterator names_i = cv->index_group_names.begin();
+ std::list<std::vector<int> >::iterator index_groups_i = cv->index_groups.begin();
+ for ( ; names_i != cv->index_group_names.end() ; ++names_i, ++index_groups_i) {
if (*names_i == index_group_name)
break;
}
- if (names_i == cvm::index_group_names.end()) {
+ if (names_i == cv->index_group_names.end()) {
cvm::error("Error: could not find index group "+
index_group_name+" among those provided by the index file.\n",
INPUT_ERROR);
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index fdd2b6254ca12dd8b450244f72c75f8b570f8a2f..3779c82aa300b86814a8f3ab23eaf1b161c34946 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -19,20 +19,6 @@ colvarbias::colvarbias(char const *key)
rank = 1;
- if (bias_type == std::string("abf")) {
- rank = cvm::n_abf_biases+1;
- }
- if (bias_type == std::string("harmonic") ||
- bias_type == std::string("linear")) {
- rank = cvm::n_rest_biases+1;
- }
- if (bias_type == std::string("histogram")) {
- rank = cvm::n_histo_biases+1;
- }
- if (bias_type == std::string("metadynamics")) {
- rank = cvm::n_meta_biases+1;
- }
-
has_data = false;
b_output_energy = false;
reset();
@@ -48,7 +34,11 @@ int colvarbias::init(std::string const &conf)
colvarparse::init(conf);
if (name.size() == 0) {
+
+ // first initialization
+
cvm::log("Initializing a new \""+bias_type+"\" instance.\n");
+ rank = cvm::num_biases_type(bias_type);
get_keyval(conf, "name", name, bias_type+cvm::to_str(rank));
{
@@ -69,7 +59,7 @@ int colvarbias::init(std::string const &conf)
// lookup the associated colvars
std::vector<std::string> colvar_names;
if (get_keyval(conf, "colvars", colvar_names)) {
- if (colvars.size()) {
+ if (num_variables()) {
cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n",
INPUT_ERROR);
return INPUT_ERROR;
@@ -80,7 +70,7 @@ int colvarbias::init(std::string const &conf)
}
}
- if (!colvars.size()) {
+ if (!num_variables()) {
cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
return INPUT_ERROR;
}
@@ -89,6 +79,8 @@ int colvarbias::init(std::string const &conf)
cvm::log("Reinitializing bias \""+name+"\".\n");
}
+ output_prefix = cvm::output_prefix();
+
get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
return COLVARS_OK;
@@ -98,7 +90,7 @@ int colvarbias::init(std::string const &conf)
int colvarbias::reset()
{
bias_energy = 0.0;
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
colvar_forces[i].reset();
}
return COLVARS_OK;
@@ -132,12 +124,13 @@ int colvarbias::clear()
}
}
+ colvarmodule *cv = cvm::main();
// ...and from the colvars module
- for (std::vector<colvarbias *>::iterator bi = cvm::biases.begin();
- bi != cvm::biases.end();
+ for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
+ bi != cv->biases.end();
++bi) {
if ( *bi == this) {
- cvm::biases.erase(bi);
+ cv->biases.erase(bi);
break;
}
}
@@ -185,21 +178,29 @@ int colvarbias::add_colvar(std::string const &cv_name)
int colvarbias::update()
{
- // Note: if anything is added here, it should be added also in the SMP block of calc_biases()
- // TODO move here debug msg of bias update
+ if (cvm::debug()) {
+ cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n");
+ }
+
has_data = true;
+
+ bias_energy = 0.0;
+ for (size_t ir = 0; ir < num_variables(); ir++) {
+ colvar_forces[ir].reset();
+ }
+
return COLVARS_OK;
}
void colvarbias::communicate_forces()
{
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+
- colvars[i]->name+"\".\n");
+ variables(i)->name+"\".\n");
}
- colvars[i]->add_bias_force(colvar_forces[i]);
+ variables(i)->add_bias_force(colvar_forces[i]);
}
}
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index 12397dcb8fa0ba8aa2356fd776d1a8ccb16b4546..6d5776d3db87047c78ff15c637d6703ceeb77300 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -32,17 +32,39 @@ public:
/// Add a new collective variable to this bias
int add_colvar(std::string const &cv_name);
- /// Add a new collective variable to this bias
- size_t number_of_colvars() const
+ /// How many variables are defined for this bias
+ inline size_t num_variables() const
{
return colvars.size();
}
+ /// Access the variables vector
+ inline std::vector<colvar *> *variables()
+ {
+ return &colvars;
+ }
+
+ /// Access the i-th variable
+ inline colvar * variables(int i) const
+ {
+ return colvars[i];
+ }
+
/// Retrieve colvar values and calculate their biasing forces
/// Return bias energy
virtual int update();
- // TODO: move update_bias here (share with metadynamics)
+ /// \brief Compute the energy of the bias with alternative values of the
+ /// collective variables (suitable for bias exchange)
+ virtual int calc_energy(std::vector<colvarvalue> const &values =
+ std::vector<colvarvalue>(0))
+ {
+ cvm::error("Error: calc_energy() not implemented.\n", COLVARS_NOT_IMPLEMENTED);
+ return COLVARS_NOT_IMPLEMENTED;
+ }
+
+ /// Send forces to the collective variables
+ virtual void communicate_forces();
/// Load new configuration - force constant and/or centers only
virtual int change_configuration(std::string const &conf);
@@ -51,10 +73,13 @@ public:
virtual cvm::real energy_difference(std::string const &conf);
/// Give the total number of bins for a given bias.
+ // FIXME this is currently 1D only
virtual int bin_num();
/// Calculate the bin index for a given bias.
+ // FIXME this is currently 1D only
virtual int current_bin();
//// Give the count at a given bin index.
+ // FIXME this is currently 1D only
virtual int bin_count(int bin_index);
//// Share information between replicas, whatever it may be.
virtual int replica_share();
@@ -62,9 +87,6 @@ public:
/// Perform analysis tasks
virtual void analyze() {}
- /// Send forces to the collective variables
- virtual void communicate_forces();
-
/// \brief Constructor
colvarbias(char const *key);
@@ -135,6 +157,9 @@ public:
return COLVARS_OK;
}
+ /// Use this prefix for all output files
+ std::string output_prefix;
+
/// If this bias is communicating with other replicas through files, send it to them
virtual int write_state_to_replicas()
{
@@ -162,7 +187,7 @@ protected:
/// through each colvar object
std::vector<colvar *> colvars;
- /// \brief Current forces from this bias to the colvars
+ /// \brief Current forces from this bias to the variables
std::vector<colvarvalue> colvar_forces;
/// \brief Current energy of this bias (colvar_forces should be obtained by deriving this)
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index ccfa401697693bb97a2d403c6cc06da3f2f32d9e..d039004f09d47f3b216c62a3e5d5304e9cb89678 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -30,9 +30,8 @@ int colvarbias_abf::init(std::string const &conf)
{
colvarbias::init(conf);
- provide(f_cvb_scalar_variables);
enable(f_cvb_scalar_variables);
- provide(f_cvb_history_dependent);
+ enable(f_cvb_calc_pmf);
// TODO relax this in case of VMD plugin
if (cvm::temperature() == 0.0)
@@ -221,9 +220,6 @@ colvarbias_abf::~colvarbias_abf()
delete [] system_force;
system_force = NULL;
}
-
- if (cvm::n_abf_biases > 0)
- cvm::n_abf_biases -= 1;
}
@@ -319,11 +315,11 @@ int colvarbias_abf::update()
}
// update the output prefix; TODO: move later to setup_output() function
- if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) {
- // This is the only ABF bias
- output_prefix = cvm::output_prefix;
+ if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
+ // This is the only bias computing PMFs
+ output_prefix = cvm::output_prefix();
} else {
- output_prefix = cvm::output_prefix + "." + this->name;
+ output_prefix = cvm::output_prefix() + "." + this->name;
}
if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp
index 9e414d4db49637087dc954ff8239a5465287664e..d096ac3daf15a508e076b38e2dc1e18368f321b7 100644
--- a/lib/colvars/colvarbias_alb.cpp
+++ b/lib/colvars/colvarbias_alb.cpp
@@ -40,11 +40,8 @@ int colvarbias_alb::init(std::string const &conf)
{
colvarbias::init(conf);
- provide(f_cvb_scalar_variables);
enable(f_cvb_scalar_variables);
- provide(f_cvb_history_dependent);
-
size_t i;
// get the initial restraint centers
@@ -138,8 +135,6 @@ int colvarbias_alb::init(std::string const &conf)
colvarbias_alb::~colvarbias_alb()
{
- if (cvm::n_rest_biases > 0)
- cvm::n_rest_biases -= 1;
}
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index 900ad213d51d6d0301063d6edba45561e285a0e7..502a7455b1ed4ab99ec9c6d553e621ec4a4ff7d7 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -24,10 +24,7 @@ int colvarbias_histogram::init(std::string const &conf)
{
colvarbias::init(conf);
- provide(f_cvb_scalar_variables);
enable(f_cvb_scalar_variables);
-
- provide(f_cvb_history_dependent);
enable(f_cvb_history_dependent);
size_t i;
@@ -104,9 +101,6 @@ colvarbias_histogram::~colvarbias_histogram()
delete grid;
grid = NULL;
}
-
- if (cvm::n_histo_biases > 0)
- cvm::n_histo_biases -= 1;
}
@@ -127,14 +121,14 @@ int colvarbias_histogram::update()
// At the first timestep, we need to assign out_name since
// output_prefix is unset during the constructor
if (cvm::step_relative() == 0) {
- out_name = cvm::output_prefix + "." + this->name + ".dat";
+ out_name = cvm::output_prefix() + "." + this->name + ".dat";
cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"");
}
}
if (out_name_dx.size() == 0) {
if (cvm::step_relative() == 0) {
- out_name_dx = cvm::output_prefix + "." + this->name + ".dx";
+ out_name_dx = cvm::output_prefix() + "." + this->name + ".dx";
cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\"");
}
}
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index 01a04d1a016dbb116f707296cfd54bb99ad86f83..b0acfe974a188a179de18da36d1c537e70f77995 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -43,7 +43,7 @@ int colvarbias_meta::init(std::string const &conf)
{
colvarbias::init(conf);
- provide(f_cvb_history_dependent);
+ enable(f_cvb_calc_pmf);
get_keyval(conf, "hillWeight", hill_weight, 0.0);
if (hill_weight > 0.0) {
@@ -59,9 +59,9 @@ int colvarbias_meta::init(std::string const &conf)
get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0);
cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
- for (size_t i = 0; i < colvars.size(); i++) {
- cvm::log(colvars[i]->name+std::string(": ")+
- cvm::to_str(0.5 * colvars[i]->width * hill_width));
+ for (size_t i = 0; i < num_variables(); i++) {
+ cvm::log(variables(i)->name+std::string(": ")+
+ cvm::to_str(0.5 * variables(i)->width * hill_width));
}
{
@@ -73,8 +73,10 @@ int colvarbias_meta::init(std::string const &conf)
comm = single_replica;
}
- // This implies gradients for all colvars
- enable(f_cvb_apply_force);
+ // in all cases, the first replica is this bias itself
+ if (replicas.size() == 0) {
+ replicas.push_back(this);
+ }
get_keyval(conf, "useGrids", use_grids, true);
@@ -84,14 +86,14 @@ int colvarbias_meta::init(std::string const &conf)
expand_grids = false;
size_t i;
- for (i = 0; i < colvars.size(); i++) {
- colvars[i]->enable(f_cv_grid);
- if (colvars[i]->expand_boundaries) {
+ for (i = 0; i < num_variables(); i++) {
+ variables(i)->enable(f_cv_grid);
+ if (variables(i)->expand_boundaries) {
expand_grids = true;
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": Will expand grids when the colvar \""+
- colvars[i]->name+"\" approaches its boundaries.\n");
+ variables(i)->name+"\" approaches its boundaries.\n");
}
}
@@ -100,7 +102,7 @@ int colvarbias_meta::init(std::string const &conf)
get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
- "please use \"keepFreeEnergyFiles\" instead.");
+ "please use \"keepFreeEnergyFile\" instead.");
}
get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
@@ -154,6 +156,20 @@ int colvarbias_meta::init(std::string const &conf)
get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false);
+ init_well_tempered_params(conf);
+ init_ebmeta_params(conf);
+
+ if (cvm::debug())
+ cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
+
+ save_delimiters = false;
+ return COLVARS_OK;
+}
+
+
+int colvarbias_meta::init_well_tempered_params(std::string const &conf)
+{
// for well-tempered metadynamics
get_keyval(conf, "wellTempered", well_tempered, false);
get_keyval(conf, "biasTemperature", bias_temperature, -1.0);
@@ -164,8 +180,12 @@ int colvarbias_meta::init(std::string const &conf)
cvm::log("Well-tempered metadynamics is used.\n");
cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
}
+ return COLVARS_OK;
+}
+int colvarbias_meta::init_ebmeta_params(std::string const &conf)
+{
// for ebmeta
target_dist = NULL;
get_keyval(conf, "ebMeta", ebmeta, false);
@@ -203,11 +223,6 @@ int colvarbias_meta::init(std::string const &conf)
get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0);
}
- if (cvm::debug())
- cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
-
- save_delimiters = false;
return COLVARS_OK;
}
@@ -234,9 +249,6 @@ colvarbias_meta::~colvarbias_meta()
delete target_dist;
target_dist = NULL;
}
-
- if (cvm::n_meta_biases > 0)
- cvm::n_meta_biases -= 1;
}
@@ -314,23 +326,45 @@ colvarbias_meta::delete_hill(hill_iter &h)
int colvarbias_meta::update()
{
- if (cvm::debug())
- cvm::log("Updating the metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
+ int error_code = COLVARS_OK;
+
+ // update base class
+ error_code |= colvarbias::update();
+
+ // update grid definition, if needed
+ error_code |= update_grid_params();
+ // add new biasing energy/forces
+ error_code |= update_bias();
+ // update grid content to reflect new bias
+ error_code |= update_grid_data();
+
+ if (comm != single_replica &&
+ (cvm::step_absolute() % replica_update_freq) == 0) {
+ // sync with the other replicas (if needed)
+ error_code |= replica_share();
+ }
+
+ error_code |= calc_energy();
+ error_code |= calc_forces();
+ return error_code;
+}
+
+
+int colvarbias_meta::update_grid_params()
+{
if (use_grids) {
std::vector<int> curr_bin = hills_energy->get_colvars_index();
+ if (cvm::debug()) {
+ cvm::log("Metadynamics bias \""+this->name+"\""+
+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+ ": current coordinates on the grid: "+
+ cvm::to_str(curr_bin)+".\n");
+ }
if (expand_grids) {
-
// first of all, expand the grids, if specified
- if (cvm::debug())
- cvm::log("Metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- ": current coordinates on the grid: "+
- cvm::to_str(curr_bin)+".\n");
-
bool changed_grids = false;
int const min_buffer =
(3 * (size_t) std::floor(hill_width)) + 1;
@@ -339,9 +373,9 @@ int colvarbias_meta::update()
std::vector<colvarvalue> new_lower_boundaries(hills_energy->lower_boundaries);
std::vector<colvarvalue> new_upper_boundaries(hills_energy->upper_boundaries);
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
- if (! colvars[i]->expand_boundaries)
+ if (! variables(i)->expand_boundaries)
continue;
cvm::real &new_lb = new_lower_boundaries[i].real_value;
@@ -349,10 +383,10 @@ int colvarbias_meta::update()
int &new_size = new_sizes[i];
bool changed_lb = false, changed_ub = false;
- if (!colvars[i]->hard_lower_boundary)
+ if (!variables(i)->hard_lower_boundary)
if (curr_bin[i] < min_buffer) {
int const extra_points = (min_buffer - curr_bin[i]);
- new_lb -= extra_points * colvars[i]->width;
+ new_lb -= extra_points * variables(i)->width;
new_size += extra_points;
// changed offset in this direction => the pointer needs to
// be changed, too
@@ -362,21 +396,21 @@ int colvarbias_meta::update()
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new lower boundary for colvar \""+
- colvars[i]->name+"\", at "+
+ variables(i)->name+"\", at "+
cvm::to_str(new_lower_boundaries[i])+".\n");
}
- if (!colvars[i]->hard_upper_boundary)
+ if (!variables(i)->hard_upper_boundary)
if (curr_bin[i] > new_size - min_buffer - 1) {
int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
- new_ub += extra_points * colvars[i]->width;
+ new_ub += extra_points * variables(i)->width;
new_size += extra_points;
changed_ub = true;
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new upper boundary for colvar \""+
- colvars[i]->name+"\", at "+
+ variables(i)->name+"\", at "+
cvm::to_str(new_upper_boundaries[i])+".\n");
}
@@ -401,7 +435,7 @@ int colvarbias_meta::update()
new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
- new_hills_energy_gradients->setup(new_sizes, 0.0, colvars.size());
+ new_hills_energy_gradients->setup(new_sizes, 0.0, num_variables());
new_hills_energy->map_grid(*hills_energy);
new_hills_energy_gradients->map_grid(*hills_energy_gradients);
@@ -418,25 +452,32 @@ int colvarbias_meta::update()
}
}
}
+ return COLVARS_OK;
+}
+
+int colvarbias_meta::update_bias()
+{
// add a new hill if the required time interval has passed
- if ((cvm::step_absolute() % new_hill_freq) == 0) {
+ if ((cvm::step_absolute() % new_hill_freq) == 0 &&
+ is_enabled(f_cvb_history_dependent)) {
- if (cvm::debug())
+ if (cvm::debug()) {
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n");
+ }
cvm::real hills_scale=1.0;
if (ebmeta) {
- hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
- if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
- cvm::real const hills_lambda =
- (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
- (cvm::real(ebmeta_equil_steps));
- hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
- }
+ hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
+ if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
+ cvm::real const hills_lambda =
+ (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
+ (cvm::real(ebmeta_equil_steps));
+ hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
+ }
}
if (well_tempered) {
@@ -471,160 +512,165 @@ int colvarbias_meta::update()
}
}
- // sync with the other replicas (if needed)
- if (comm != single_replica) {
+ return COLVARS_OK;
+}
- // reread the replicas registry
- if ((cvm::step_absolute() % replica_update_freq) == 0) {
- update_replicas_registry();
- // empty the output buffer
- if (replica_hills_os.is_open())
- replica_hills_os.flush();
- read_replica_files();
+int colvarbias_meta::update_grid_data()
+{
+ if ((cvm::step_absolute() % grids_freq) == 0) {
+ // map the most recent gaussians to the grids
+ project_hills(new_hills_begin, hills.end(),
+ hills_energy, hills_energy_gradients);
+ new_hills_begin = hills.end();
+
+ // TODO: we may want to condense all into one replicas array,
+ // including "this" as the first element
+ if (comm == multiple_replicas) {
+ for (size_t ir = 0; ir < replicas.size(); ir++) {
+ replicas[ir]->project_hills(replicas[ir]->new_hills_begin,
+ replicas[ir]->hills.end(),
+ replicas[ir]->hills_energy,
+ replicas[ir]->hills_energy_gradients);
+ replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
+ }
}
}
- // calculate the biasing energy and forces
- bias_energy = 0.0;
- for (size_t ir = 0; ir < colvars.size(); ir++) {
- colvar_forces[ir].reset();
+ return COLVARS_OK;
+}
+
+
+int colvarbias_meta::calc_energy(std::vector<colvarvalue> const &values)
+{
+ size_t ir = 0;
+
+ for (ir = 0; ir < replicas.size(); ir++) {
+ replicas[ir]->bias_energy = 0.0;
}
- if (comm == multiple_replicas)
- for (size_t ir = 0; ir < replicas.size(); ir++) {
- replicas[ir]->bias_energy = 0.0;
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- replicas[ir]->colvar_forces[ic].reset();
- }
- }
- if (use_grids) {
+ std::vector<int> const curr_bin = values.size() ?
+ hills_energy->get_colvars_index(values) :
+ hills_energy->get_colvars_index();
- // get the forces from the grid
-
- if ((cvm::step_absolute() % grids_freq) == 0) {
- // map the most recent gaussians to the grids
- project_hills(new_hills_begin, hills.end(),
- hills_energy, hills_energy_gradients);
- new_hills_begin = hills.end();
-
- // TODO: we may want to condense all into one replicas array,
- // including "this" as the first element
- if (comm == multiple_replicas) {
- for (size_t ir = 0; ir < replicas.size(); ir++) {
- replicas[ir]->project_hills(replicas[ir]->new_hills_begin,
- replicas[ir]->hills.end(),
- replicas[ir]->hills_energy,
- replicas[ir]->hills_energy_gradients);
- replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
- }
+ if (hills_energy->index_ok(curr_bin)) {
+ // index is within the grid: get the energy from there
+ for (ir = 0; ir < replicas.size(); ir++) {
+
+ bias_energy += replicas[ir]->hills_energy->value(curr_bin);
+ if (cvm::debug()) {
+ cvm::log("Metadynamics bias \""+this->name+"\""+
+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+ ": current coordinates on the grid: "+
+ cvm::to_str(curr_bin)+".\n");
+ cvm::log("Grid energy = "+cvm::to_str(bias_energy)+".\n");
}
}
+ } else {
+ // off the grid: compute analytically only the hills at the grid's edges
+ for (ir = 0; ir < replicas.size(); ir++) {
+ calc_hills(replicas[ir]->hills_off_grid.begin(),
+ replicas[ir]->hills_off_grid.end(),
+ bias_energy,
+ values);
+ }
+ }
- std::vector<int> curr_bin = hills_energy->get_colvars_index();
- if (cvm::debug())
- cvm::log("Metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- ": current coordinates on the grid: "+
- cvm::to_str(curr_bin)+".\n");
+ // now include the hills that have not been binned yet (starting
+ // from new_hills_begin)
- if (hills_energy->index_ok(curr_bin)) {
+ for (ir = 0; ir < replicas.size(); ir++) {
+ calc_hills(replicas[ir]->new_hills_begin,
+ replicas[ir]->hills.end(),
+ bias_energy);
+ if (cvm::debug()) {
+ cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n");
+ }
+ }
- // within the grid: add the energy and the forces from there
+ return COLVARS_OK;
+}
- bias_energy += hills_energy->value(curr_bin);
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- cvm::real const *f = &(hills_energy_gradients->value(curr_bin));
- colvar_forces[ic].real_value += -1.0 * f[ic];
- // the gradients are stored, not the forces
- }
- if (comm == multiple_replicas)
- for (size_t ir = 0; ir < replicas.size(); ir++) {
- bias_energy += replicas[ir]->hills_energy->value(curr_bin);
- cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- colvar_forces[ic].real_value += -1.0 * f[ic];
- }
- }
- } else {
+int colvarbias_meta::calc_forces(std::vector<colvarvalue> const &values)
+{
+ size_t ir = 0, ic = 0;
+ for (ir = 0; ir < replicas.size(); ir++) {
+ for (ic = 0; ic < num_variables(); ic++) {
+ replicas[ir]->colvar_forces[ic].reset();
+ }
+ }
- // off the grid: compute analytically only the hills at the grid's edges
+ std::vector<int> const curr_bin = values.size() ?
+ hills_energy->get_colvars_index(values) :
+ hills_energy->get_colvars_index();
- calc_hills(hills_off_grid.begin(), hills_off_grid.end(), bias_energy);
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- calc_hills_force(ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
+ if (hills_energy->index_ok(curr_bin)) {
+ for (ir = 0; ir < replicas.size(); ir++) {
+ cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
+ for (ic = 0; ic < num_variables(); ic++) {
+ // the gradients are stored, not the forces
+ colvar_forces[ic].real_value += -1.0 * f[ic];
+ }
+ }
+ } else {
+ // off the grid: compute analytically only the hills at the grid's edges
+ for (ir = 0; ir < replicas.size(); ir++) {
+ for (ic = 0; ic < num_variables(); ic++) {
+ calc_hills_force(ic,
+ replicas[ir]->hills_off_grid.begin(),
+ replicas[ir]->hills_off_grid.end(),
+ colvar_forces,
+ values);
}
-
- if (comm == multiple_replicas)
- for (size_t ir = 0; ir < replicas.size(); ir++) {
- calc_hills(replicas[ir]->hills_off_grid.begin(),
- replicas[ir]->hills_off_grid.end(),
- bias_energy);
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- calc_hills_force(ic,
- replicas[ir]->hills_off_grid.begin(),
- replicas[ir]->hills_off_grid.end(),
- colvar_forces);
- }
- }
}
}
// now include the hills that have not been binned yet (starting
// from new_hills_begin)
- calc_hills(new_hills_begin, hills.end(), bias_energy);
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- calc_hills_force(ic, new_hills_begin, hills.end(), colvar_forces);
- }
-
- if (cvm::debug())
- cvm::log("Hills energy = "+cvm::to_str(bias_energy)+
- ", hills forces = "+cvm::to_str(colvar_forces)+".\n");
-
- if (cvm::debug())
+ if (cvm::debug()) {
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding the forces from the other replicas.\n");
+ }
- if (comm == multiple_replicas)
- for (size_t ir = 0; ir < replicas.size(); ir++) {
- calc_hills(replicas[ir]->new_hills_begin,
- replicas[ir]->hills.end(),
- bias_energy);
- for (size_t ic = 0; ic < colvars.size(); ic++) {
- calc_hills_force(ic,
- replicas[ir]->new_hills_begin,
- replicas[ir]->hills.end(),
- colvar_forces);
+ for (ir = 0; ir < replicas.size(); ir++) {
+ for (ic = 0; ic < num_variables(); ic++) {
+ calc_hills_force(ic,
+ replicas[ir]->new_hills_begin,
+ replicas[ir]->hills.end(),
+ colvar_forces,
+ values);
+ if (cvm::debug()) {
+ cvm::log("Hills forces = "+cvm::to_str(colvar_forces)+".\n");
}
- if (cvm::debug())
- cvm::log("Hills energy = "+cvm::to_str(bias_energy)+
- ", hills forces = "+cvm::to_str(colvar_forces)+".\n");
}
+ }
return COLVARS_OK;
}
+
void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
cvm::real &energy,
std::vector<colvarvalue> const &colvar_values)
{
- std::vector<colvarvalue> curr_values(colvars.size());
- for (size_t i = 0; i < colvars.size(); i++) {
- curr_values[i].type(colvars[i]->value());
+ size_t i = 0;
+ std::vector<colvarvalue> curr_values(num_variables());
+ for (i = 0; i < num_variables(); i++) {
+ curr_values[i].type(variables(i)->value());
}
if (colvar_values.size()) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
curr_values[i] = colvar_values[i];
}
} else {
- for (size_t i = 0; i < colvars.size(); i++) {
- curr_values[i] = colvars[i]->value();
+ for (i = 0; i < num_variables(); i++) {
+ curr_values[i] = variables(i)->value();
}
}
@@ -632,11 +678,11 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first,
// compute the gaussian exponent
cvm::real cv_sqdev = 0.0;
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvarvalue const &x = curr_values[i];
colvarvalue const ¢er = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[i];
- cv_sqdev += (colvars[i]->dist2(x, center)) / (half_width*half_width);
+ cv_sqdev += (variables(i)->dist2(x, center)) / (half_width*half_width);
}
// compute the gaussian
@@ -658,14 +704,14 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
std::vector<colvarvalue> const &values)
{
// Retrieve the value of the colvar
- colvarvalue const x(values.size() ? values[i] : colvars[i]->value());
+ colvarvalue const x(values.size() ? values[i] : variables(i)->value());
// do the type check only once (all colvarvalues in the hills series
// were already saved with their types matching those in the
// colvars)
hill_iter h;
- switch (colvars[i]->value().type()) {
+ switch (variables(i)->value().type()) {
case colvarvalue::type_scalar:
for (h = h_first; h != h_last; h++) {
@@ -674,7 +720,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].real_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
- (colvars[i]->dist2_lgrad(x, center)).real_value );
+ (variables(i)->dist2_lgrad(x, center)).real_value );
}
break;
@@ -687,7 +733,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].rvector_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
- (colvars[i]->dist2_lgrad(x, center)).rvector_value );
+ (variables(i)->dist2_lgrad(x, center)).rvector_value );
}
break;
@@ -699,7 +745,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].quaternion_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
- (colvars[i]->dist2_lgrad(x, center)).quaternion_value );
+ (variables(i)->dist2_lgrad(x, center)).quaternion_value );
}
break;
@@ -710,7 +756,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].vector1d_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
- (colvars[i]->dist2_lgrad(x, center)).vector1d_value );
+ (variables(i)->dist2_lgrad(x, center)).vector1d_value );
}
break;
@@ -739,13 +785,13 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
// TODO: improve it by looping over a small subgrid instead of the whole grid
- std::vector<colvarvalue> colvar_values(colvars.size());
- std::vector<cvm::real> colvar_forces_scalar(colvars.size());
+ std::vector<colvarvalue> colvar_values(num_variables());
+ std::vector<cvm::real> colvar_forces_scalar(num_variables());
std::vector<int> he_ix = he->new_index();
std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
cvm::real hills_energy_here = 0.0;
- std::vector<colvarvalue> hills_forces_here(colvars.size(), 0.0);
+ std::vector<colvarvalue> hills_forces_here(num_variables(), 0.0);
size_t count = 0;
size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
@@ -757,7 +803,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
(he->index_ok(he_ix)) && (hg->index_ok(hg_ix));
count++) {
size_t i;
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
}
@@ -766,7 +812,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
calc_hills(h_first, h_last, hills_energy_here, colvar_values);
he->acc_value(he_ix, hills_energy_here);
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
hills_forces_here[i].reset();
calc_hills_force(i, h_first, h_last, hills_forces_here, colvar_values);
colvar_forces_scalar[i] = hills_forces_here[i].real_value;
@@ -795,7 +841,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
for ( ; (he->index_ok(he_ix)); ) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
}
@@ -851,6 +897,21 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first
// **********************************************************************
+int colvarbias_meta::replica_share()
+{
+ // sync with the other replicas (if needed)
+ if (comm == multiple_replicas) {
+ // reread the replicas registry
+ update_replicas_registry();
+ // empty the output buffer
+ if (replica_hills_os.is_open())
+ replica_hills_os.flush();
+ read_replica_files();
+ }
+ return COLVARS_OK;
+}
+
+
void colvarbias_meta::update_replicas_registry()
{
if (cvm::debug())
@@ -975,7 +1036,6 @@ void colvarbias_meta::update_replicas_registry()
}
}
-
if (cvm::debug())
cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
cvm::to_str(replicas.size())+" elements.\n");
@@ -984,7 +1044,8 @@ void colvarbias_meta::update_replicas_registry()
void colvarbias_meta::read_replica_files()
{
- for (size_t ir = 0; ir < replicas.size(); ir++) {
+ // Note: we start from the 2nd replica.
+ for (size_t ir = 1; ir < replicas.size(); ir++) {
if (! (replicas[ir])->replica_state_file_in_sync) {
// if a new state file is being read, the hills file is also new
@@ -1352,9 +1413,9 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
cvm::real h_weight;
get_keyval(data, "weight", h_weight, hill_weight, parse_silent);
- std::vector<colvarvalue> h_centers(colvars.size());
- for (size_t i = 0; i < colvars.size(); i++) {
- h_centers[i].type(colvars[i]->value());
+ std::vector<colvarvalue> h_centers(num_variables());
+ for (size_t i = 0; i < num_variables(); i++) {
+ h_centers[i].type(variables(i)->value());
}
{
// it is safer to read colvarvalue objects one at a time;
@@ -1362,14 +1423,14 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
std::string centers_input;
key_lookup(data, "centers", centers_input);
std::istringstream centers_is(centers_input);
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
centers_is >> h_centers[i];
}
}
- std::vector<cvm::real> h_widths(colvars.size());
+ std::vector<cvm::real> h_widths(num_variables());
get_keyval(data, "widths", h_widths,
- std::vector<cvm::real> (colvars.size(), (std::sqrt(2.0 * PI) / 2.0)),
+ std::vector<cvm::real>(num_variables(), (std::sqrt(2.0 * PI) / 2.0)),
parse_silent);
std::string h_replica = "";
@@ -1406,6 +1467,13 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
int colvarbias_meta::setup_output()
{
+ output_prefix = cvm::output_prefix();
+ if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) {
+ // if this is not the only free energy integrator, append
+ // this bias's name, to distinguish it from the output of the other
+ // biases producing a .pmf file
+ output_prefix += ("."+this->name);
+ }
if (comm == multiple_replicas) {
@@ -1421,10 +1489,10 @@ int colvarbias_meta::setup_output()
// those by another replica
replica_hills_file =
(std::string(pwd)+std::string(PATHSEP)+
- cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills");
+ cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills");
replica_state_file =
(std::string(pwd)+std::string(PATHSEP)+
- cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state");
+ cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state");
delete[] pwd;
// now register this replica
@@ -1492,13 +1560,14 @@ int colvarbias_meta::setup_output()
}
if (b_hills_traj) {
- std::string const traj_file_name(cvm::output_prefix+
+ std::string const traj_file_name(cvm::output_prefix()+
".colvars."+this->name+
( (comm != single_replica) ?
("."+replica_id) :
("") )+
".hills.traj");
if (!hills_traj_os.is_open()) {
+ cvm::backup_file(traj_file_name.c_str());
hills_traj_os.open(traj_file_name.c_str());
}
if (!hills_traj_os.is_open())
@@ -1585,16 +1654,6 @@ void colvarbias_meta::write_pmf()
colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy);
pmf->setup();
- std::string fes_file_name_prefix(cvm::output_prefix);
-
- if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) {
- // if this is not the only free energy integrator, append
- // this bias's name, to distinguish it from the output of the other
- // biases producing a .pmf file
- // TODO: fix for ABF with updateBias == no
- fes_file_name_prefix += ("."+this->name);
- }
-
if ((comm == single_replica) || (dump_replica_fes)) {
// output the PMF from this instance or replica
pmf->reset();
@@ -1607,7 +1666,7 @@ void colvarbias_meta::write_pmf()
pmf->multiply_constant(well_temper_scale);
}
{
- std::string const fes_file_name(fes_file_name_prefix +
+ std::string const fes_file_name(this->output_prefix +
((comm != single_replica) ? ".partial" : "") +
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
@@ -1632,7 +1691,7 @@ void colvarbias_meta::write_pmf()
cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
pmf->multiply_constant(well_temper_scale);
}
- std::string const fes_file_name(fes_file_name_prefix +
+ std::string const fes_file_name(this->output_prefix +
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h
index a88a34ba00a7749895390cd81e5179687ebb5599..01921eaf640855f40fabcd8147b4815e7579cd89 100644
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@@ -36,8 +36,20 @@ public:
colvarbias_meta(char const *key);
virtual int init(std::string const &conf);
+ virtual int init_well_tempered_params(std::string const &conf);
+ virtual int init_ebmeta_params(std::string const &conf);
virtual ~colvarbias_meta();
+
virtual int update();
+ virtual int update_grid_params();
+ virtual int update_bias();
+ virtual int update_grid_data();
+ virtual int replica_share();
+
+ virtual int calc_energy(std::vector<colvarvalue> const &values =
+ std::vector<colvarvalue>(0));
+ virtual int calc_forces(std::vector<colvarvalue> const &values =
+ std::vector<colvarvalue>(0));
virtual std::string const get_state_params() const;
virtual int set_state_params(std::string const &state_conf);
@@ -102,18 +114,18 @@ protected:
/// \brief Calculate the values of the hills, incrementing
/// bias_energy
virtual void calc_hills(hill_iter h_first,
- hill_iter h_last,
- cvm::real &energy,
- std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
+ hill_iter h_last,
+ cvm::real &energy,
+ std::vector<colvarvalue> const &values = std::vector<colvarvalue>(0));
/// \brief Calculate the forces acting on the i-th colvar,
/// incrementing colvar_forces[i]; must be called after calc_hills
/// each time the values of the colvars are changed
virtual void calc_hills_force(size_t const &i,
- hill_iter h_first,
- hill_iter h_last,
- std::vector<colvarvalue> &forces,
- std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
+ hill_iter h_first,
+ hill_iter h_last,
+ std::vector<colvarvalue> &forces,
+ std::vector<colvarvalue> const &values = std::vector<colvarvalue>(0));
/// Height of new hills
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index 84630984e5085fb9b07e96ca9ad8da9480608b62..159d9eae6451bbfc9a307d8507d64888188c5bc7 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -33,17 +33,15 @@ int colvarbias_restraint::init(std::string const &conf)
int colvarbias_restraint::update()
{
- bias_energy = 0.0;
-
- if (cvm::debug())
- cvm::log("Updating the restraint bias \""+this->name+"\".\n");
+ // Update base class (bias_energy and colvar_forces are zeroed there)
+ colvarbias::update();
// Force and energy calculation
- for (size_t i = 0; i < colvars.size(); i++) {
- colvar_forces[i].type(colvars[i]->value());
+ for (size_t i = 0; i < num_variables(); i++) {
+ bias_energy += restraint_potential(i);
+ colvar_forces[i].type(variables(i)->value());
colvar_forces[i].is_derivative();
colvar_forces[i] = restraint_force(i);
- bias_energy += restraint_potential(i);
}
if (cvm::debug())
@@ -59,8 +57,6 @@ int colvarbias_restraint::update()
colvarbias_restraint::~colvarbias_restraint()
{
- if (cvm::n_rest_biases > 0)
- cvm::n_rest_biases -= 1;
}
@@ -102,18 +98,18 @@ int colvarbias_restraint_centers::init(std::string const &conf)
bool null_centers = (colvar_centers.size() == 0);
if (null_centers) {
// try to initialize the restraint centers for the first time
- colvar_centers.resize(colvars.size());
- colvar_centers_raw.resize(colvars.size());
- for (i = 0; i < colvars.size(); i++) {
- colvar_centers[i].type(colvars[i]->value());
+ colvar_centers.resize(num_variables());
+ colvar_centers_raw.resize(num_variables());
+ for (i = 0; i < num_variables(); i++) {
+ colvar_centers[i].type(variables(i)->value());
colvar_centers[i].reset();
- colvar_centers_raw[i].type(colvars[i]->value());
+ colvar_centers_raw[i].type(variables(i)->value());
colvar_centers_raw[i].reset();
}
}
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
if (cvm::debug()) {
cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n");
}
@@ -129,7 +125,7 @@ int colvarbias_restraint_centers::init(std::string const &conf)
return INPUT_ERROR;
}
- if (colvar_centers.size() != colvars.size()) {
+ if (colvar_centers.size() != num_variables()) {
cvm::error("Error: number of centers does not match "
"that of collective variables.\n", INPUT_ERROR);
return INPUT_ERROR;
@@ -142,10 +138,10 @@ int colvarbias_restraint_centers::init(std::string const &conf)
int colvarbias_restraint_centers::change_configuration(std::string const &conf)
{
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
- for (size_t i = 0; i < colvars.size(); i++) {
- colvar_centers[i].type(colvars[i]->value());
+ for (size_t i = 0; i < num_variables(); i++) {
+ colvar_centers[i].type(variables(i)->value());
colvar_centers[i].apply_constraints();
- colvar_centers_raw[i].type(colvars[i]->value());
+ colvar_centers_raw[i].type(variables(i)->value());
colvar_centers_raw[i] = colvar_centers[i];
}
}
@@ -269,7 +265,7 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf)
size_t i;
if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
- if (colvar_centers.size() != colvars.size()) {
+ if (colvar_centers.size() != num_variables()) {
cvm::error("Error: number of target centers does not match "
"that of collective variables.\n");
}
@@ -308,9 +304,9 @@ int colvarbias_restraint_centers_moving::update()
// at each simulation step (or stage, if applicable)
// (take current stage into account: it can be non-zero
// if we are restarting a staged calculation)
- centers_incr.resize(colvars.size());
- for (size_t i = 0; i < colvars.size(); i++) {
- centers_incr[i].type(colvars[i]->value());
+ centers_incr.resize(num_variables());
+ for (size_t i = 0; i < num_variables(); i++) {
+ centers_incr[i].type(variables(i)->value());
centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) /
cvm::real( target_nstages ? (target_nstages - stage) :
(target_nsteps - cvm::step_absolute()));
@@ -326,10 +322,10 @@ int colvarbias_restraint_centers_moving::update()
&& (cvm::step_absolute() % target_nsteps) == 0
&& stage < target_nstages) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
colvar_centers_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[i];
- colvars[i]->wrap(colvar_centers[i]);
+ variables(i)->wrap(colvar_centers[i]);
colvar_centers[i].apply_constraints();
}
stage++;
@@ -341,10 +337,10 @@ int colvarbias_restraint_centers_moving::update()
} else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) {
// move the restraint centers in the direction of the targets
// (slow growth)
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
colvar_centers_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[i];
- colvars[i]->wrap(colvar_centers[i]);
+ variables(i)->wrap(colvar_centers[i]);
colvar_centers[i].apply_constraints();
}
}
@@ -363,7 +359,7 @@ int colvarbias_restraint_centers_moving::update_acc_work()
{
if (b_output_acc_work) {
if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
// project forces on the calculated increments at this step
acc_work += colvar_forces[i] * centers_incr[i];
}
@@ -381,14 +377,14 @@ std::string const colvarbias_restraint_centers_moving::get_state_params() const
if (b_chg_centers) {
size_t i;
os << "centers ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers[i];
}
os << "\n";
os << "centers_raw ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers_raw[i];
@@ -429,10 +425,10 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con
std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os)
{
if (b_output_centers) {
- for (size_t i = 0; i < colvars.size(); i++) {
- size_t const this_cv_width = (colvars[i]->value()).output_width(cvm::cv_width);
+ for (size_t i = 0; i < num_variables(); i++) {
+ size_t const this_cv_width = (variables(i)->value()).output_width(cvm::cv_width);
os << " x0_"
- << cvm::wrap_string(colvars[i]->name, this_cv_width-3);
+ << cvm::wrap_string(variables(i)->name, this_cv_width-3);
}
}
@@ -448,7 +444,7 @@ std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostrea
std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
{
if (b_output_centers) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers[i];
@@ -539,9 +535,9 @@ int colvarbias_restraint_k_moving::update()
}
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow(lambda, force_k_exp);
- cvm::log("Restraint " + this->name + ", stage " +
- cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda));
- cvm::log("Setting force constant to " + cvm::to_str(force_k));
+ cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
+ + " : lambda = " + cvm::to_str(lambda)
+ + ", k = " + cvm::to_str(force_k));
}
// TI calculation: estimate free energy derivative
@@ -557,7 +553,7 @@ int colvarbias_restraint_k_moving::update()
// Square distance normalized by square colvar width
cvm::real dist_sq = 0.0;
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
dist_sq += d_restraint_potential_dk(i);
}
@@ -569,7 +565,8 @@ int colvarbias_restraint_k_moving::update()
if (cvm::step_absolute() % target_nsteps == 0 &&
cvm::step_absolute() > 0) {
- cvm::log("Lambda= " + cvm::to_str(lambda) + " dA/dLambda= "
+ cvm::log("Restraint " + this->name + " Lambda= "
+ + cvm::to_str(lambda) + " dA/dLambda= "
+ cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps)));
// ...and move on to the next one
@@ -584,9 +581,9 @@ int colvarbias_restraint_k_moving::update()
}
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow(lambda, force_k_exp);
- cvm::log("Restraint " + this->name + ", stage " +
- cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda));
- cvm::log("Setting force constant to " + cvm::to_str(force_k));
+ cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
+ + " : lambda = " + cvm::to_str(lambda)
+ + ", k = " + cvm::to_str(force_k));
}
}
@@ -721,11 +718,11 @@ int colvarbias_restraint_harmonic::init(std::string const &conf)
colvarbias_restraint_centers_moving::init(conf);
colvarbias_restraint_k_moving::init(conf);
- for (size_t i = 0; i < colvars.size(); i++) {
- if (colvars[i]->width != 1.0)
- cvm::log("The force constant for colvar \""+colvars[i]->name+
+ for (size_t i = 0; i < num_variables(); i++) {
+ if (variables(i)->width != 1.0)
+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+
- cvm::to_str(force_k / (colvars[i]->width * colvars[i]->width))+
+ cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+
" according to the specified width.\n");
}
@@ -751,22 +748,22 @@ int colvarbias_restraint_harmonic::update()
cvm::real colvarbias_restraint_harmonic::restraint_potential(size_t i) const
{
- return 0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
- colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]);
+ return 0.5 * force_k / (variables(i)->width * variables(i)->width) *
+ variables(i)->dist2(variables(i)->value(), colvar_centers[i]);
}
colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const
{
- return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
- colvars[i]->dist2_lgrad(colvars[i]->value(), colvar_centers[i]);
+ return -0.5 * force_k / (variables(i)->width * variables(i)->width) *
+ variables(i)->dist2_lgrad(variables(i)->value(), colvar_centers[i]);
}
cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const
{
- return 0.5 / (colvars[i]->width * colvars[i]->width) *
- colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]);
+ return 0.5 / (variables(i)->width * variables(i)->width) *
+ variables(i)->dist2(variables(i)->value(), colvar_centers[i]);
}
@@ -840,6 +837,8 @@ colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char co
colvarbias_restraint_moving(key),
colvarbias_restraint_k_moving(key)
{
+ lower_wall_k = 0.0;
+ upper_wall_k = 0.0;
}
@@ -849,7 +848,11 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
colvarbias_restraint_moving::init(conf);
colvarbias_restraint_k_moving::init(conf);
- provide(f_cvb_scalar_variables);
+ get_keyval(conf, "lowerWallConstant", lower_wall_k,
+ (lower_wall_k > 0.0) ? lower_wall_k : force_k);
+ get_keyval(conf, "upperWallConstant", upper_wall_k,
+ (upper_wall_k > 0.0) ? upper_wall_k : force_k);
+
enable(f_cvb_scalar_variables);
size_t i;
@@ -857,9 +860,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
bool b_null_lower_walls = false;
if (lower_walls.size() == 0) {
b_null_lower_walls = true;
- lower_walls.resize(number_of_colvars());
- for (i = 0; i < colvars.size(); i++) {
- lower_walls[i].type(colvars[i]->value());
+ lower_walls.resize(num_variables());
+ for (i = 0; i < num_variables(); i++) {
+ lower_walls[i].type(variables(i)->value());
lower_walls[i].reset();
}
}
@@ -872,9 +875,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
bool b_null_upper_walls = false;
if (upper_walls.size() == 0) {
b_null_upper_walls = true;
- upper_walls.resize(number_of_colvars());
- for (i = 0; i < colvars.size(); i++) {
- upper_walls[i].type(colvars[i]->value());
+ upper_walls.resize(num_variables());
+ for (i = 0; i < num_variables(); i++) {
+ upper_walls[i].type(variables(i)->value());
upper_walls[i].reset();
}
}
@@ -890,17 +893,17 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
}
if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) {
- for (i = 0; i < colvars.size(); i++) {
- if (colvars[i]->is_enabled(f_cv_periodic)) {
+ for (i = 0; i < num_variables(); i++) {
+ if (variables(i)->is_enabled(f_cv_periodic)) {
cvm::error("Error: at least one variable is periodic, "
- "both walls must be provided .\n", INPUT_ERROR);
+ "both walls must be provided.\n", INPUT_ERROR);
return INPUT_ERROR;
}
}
}
if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) {
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
if (lower_walls[i] >= upper_walls[i]) {
cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+
@@ -909,13 +912,24 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
INPUT_ERROR);
}
}
+ if (lower_wall_k * upper_wall_k == 0.0) {
+ cvm::error("Error: lowerWallConstant and upperWallConstant, "
+ "when defined, must both be positive.\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ force_k = lower_wall_k * upper_wall_k;
+ // transform the two constants to relative values
+ // (allow changing both via force_k)
+ lower_wall_k /= force_k;
+ upper_wall_k /= force_k;
}
- for (i = 0; i < colvars.size(); i++) {
- if (colvars[i]->width != 1.0)
- cvm::log("The force constant for colvar \""+colvars[i]->name+
+ for (i = 0; i < num_variables(); i++) {
+ if (variables(i)->width != 1.0)
+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+
- cvm::to_str(force_k / (colvars[i]->width * colvars[i]->width))+
+ cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+
" according to the specified width.\n");
}
@@ -935,20 +949,20 @@ int colvarbias_restraint_harmonic_walls::update()
void colvarbias_restraint_harmonic_walls::communicate_forces()
{
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+
- colvars[i]->name+"\".\n");
+ variables(i)->name+"\".\n");
}
- colvars[i]->add_bias_force_actual_value(colvar_forces[i]);
+ variables(i)->add_bias_force_actual_value(colvar_forces[i]);
}
}
cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
{
- colvar *cv = colvars[i];
- colvarvalue const &cvv = colvars[i]->actual_value();
+ colvar *cv = variables(i);
+ colvarvalue const &cvv = variables(i)->actual_value();
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
@@ -958,21 +972,21 @@ cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
cvm::real const upper_wall_dist2 = cv->dist2(cvv, upper_walls[i]);
if (lower_wall_dist2 < upper_wall_dist2) {
cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]);
- if (grad < 0.0) { return grad; }
+ if (grad < 0.0) { return 0.5 * grad; }
} else {
cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
- if (grad > 0.0) { return grad; }
+ if (grad > 0.0) { return 0.5 * grad; }
}
return 0.0;
}
if (lower_walls.size() > 0) {
cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]);
- if (grad < 0.0) { return grad; }
+ if (grad < 0.0) { return 0.5 * grad; }
}
if (upper_walls.size() > 0) {
cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
- if (grad > 0.0) { return grad; }
+ if (grad > 0.0) { return 0.5 * grad; }
}
return 0.0;
}
@@ -981,7 +995,8 @@ cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) const
{
cvm::real const dist = colvar_distance(i);
- return 0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
+ cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
+ return 0.5 * force_k * scale / (variables(i)->width * variables(i)->width) *
dist * dist;
}
@@ -989,15 +1004,16 @@ cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) con
colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const
{
cvm::real const dist = colvar_distance(i);
- return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
- dist;
+ cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
+ return - force_k * scale / (variables(i)->width * variables(i)->width) * dist;
}
cvm::real colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i) const
{
cvm::real const dist = colvar_distance(i);
- return 0.5 / (colvars[i]->width * colvars[i]->width) *
+ cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
+ return 0.5 * scale / (variables(i)->width * variables(i)->width) *
dist * dist;
}
@@ -1054,16 +1070,16 @@ int colvarbias_restraint_linear::init(std::string const &conf)
colvarbias_restraint_centers_moving::init(conf);
colvarbias_restraint_k_moving::init(conf);
- for (size_t i = 0; i < colvars.size(); i++) {
- if (colvars[i]->is_enabled(f_cv_periodic)) {
+ for (size_t i = 0; i < num_variables(); i++) {
+ if (variables(i)->is_enabled(f_cv_periodic)) {
cvm::error("Error: linear biases cannot be applied to periodic variables.\n",
INPUT_ERROR);
return INPUT_ERROR;
}
- if (colvars[i]->width != 1.0)
- cvm::log("The force constant for colvar \""+colvars[i]->name+
+ if (variables(i)->width != 1.0)
+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+
- cvm::to_str(force_k / colvars[i]->width)+
+ cvm::to_str(force_k / variables(i)->width)+
" according to the specified width.\n");
}
@@ -1113,19 +1129,19 @@ cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf
cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
{
- return force_k / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]);
+ return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
}
colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
{
- return -1.0 * force_k / colvars[i]->width;
+ return -1.0 * force_k / variables(i)->width;
}
cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const
{
- return 1.0 / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]);
+ return 1.0 / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
}
@@ -1279,16 +1295,16 @@ int colvarbias_restraint_histogram::update()
size_t vector_size = 0;
size_t icv;
- for (icv = 0; icv < colvars.size(); icv++) {
- vector_size += colvars[icv]->value().size();
+ for (icv = 0; icv < num_variables(); icv++) {
+ vector_size += variables(icv)->value().size();
}
cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
// calculate the histogram
p.reset();
- for (icv = 0; icv < colvars.size(); icv++) {
- colvarvalue const &cv = colvars[icv]->value();
+ for (icv = 0; icv < num_variables(); icv++) {
+ colvarvalue const &cv = variables(icv)->value();
if (cv.type() == colvarvalue::type_scalar) {
cvm::real const cv_value = cv.real_value;
size_t igrid;
@@ -1309,7 +1325,9 @@ int colvarbias_restraint_histogram::update()
}
}
} else {
- // TODO
+ cvm::error("Error: unsupported type for variable "+variables(icv)->name+".\n",
+ COLVARS_NOT_IMPLEMENTED);
+ return COLVARS_NOT_IMPLEMENTED;
}
}
@@ -1320,8 +1338,8 @@ int colvarbias_restraint_histogram::update()
bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
// calculate the forces
- for (icv = 0; icv < colvars.size(); icv++) {
- colvarvalue const &cv = colvars[icv]->value();
+ for (icv = 0; icv < num_variables(); icv++) {
+ colvarvalue const &cv = variables(icv)->value();
colvarvalue &cv_force = colvar_forces[icv];
cv_force.type(cv);
cv_force.reset();
@@ -1363,10 +1381,10 @@ int colvarbias_restraint_histogram::update()
std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
{
if (b_write_histogram) {
- std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat");
+ std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat");
std::ofstream os(file_name.c_str());
- os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width)
- << " " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3)
+ os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
+ << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3)
<< ")\n";
size_t igrid;
for (igrid = 0; igrid < p.size(); igrid++) {
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index cbdd9c44d163a3c4684eafafaeb07cfa0a9b0495..98b967abdb1ba6a50754e2f26c3d9600e742557a 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -260,6 +260,12 @@ protected:
/// \brief Location of the upper walls
std::vector<colvarvalue> upper_walls;
+ /// \brief If both walls are defined, use this k for the lower
+ cvm::real lower_wall_k;
+
+ /// \brief If both walls are defined, use this k for the upper
+ cvm::real upper_wall_k;
+
virtual cvm::real colvar_distance(size_t i) const;
virtual cvm::real restraint_potential(size_t i) const;
virtual colvarvalue const restraint_force(size_t i) const;
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index d99da56386f05970d4d90f30d1cee27811390aa0..786bc032d2ab11e5eb466c1ef3b43dda01809560 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -48,8 +48,6 @@ colvar::cvc::cvc(std::string const &conf)
get_keyval(conf, "period", period, 0.0);
get_keyval(conf, "wrapAround", wrap_center, 0.0);
- // All cvcs implement this
- provide(f_cvc_debug_gradient);
get_keyval_feature((colvarparse *)this, conf, "debugGradients",
f_cvc_debug_gradient, false, parse_silent);
@@ -63,6 +61,8 @@ colvar::cvc::cvc(std::string const &conf)
int colvar::cvc::init_total_force_params(std::string const &conf)
{
+ if (cvm::get_error()) return COLVARS_ERROR;
+
if (get_keyval_feature(this, conf, "oneSiteSystemForce",
f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
@@ -72,6 +72,19 @@ int colvar::cvc::init_total_force_params(std::string const &conf)
f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
cvm::log("Computing total force on group 1 only");
}
+
+ if (! is_enabled(f_cvc_one_site_total_force)) {
+ // check whether any of the other atom groups is dummy
+ std::vector<cvm::atom_group *>::iterator agi = atom_groups.begin();
+ agi++;
+ for ( ; agi != atom_groups.end(); agi++) {
+ if ((*agi)->b_dummy) {
+ provide(f_cvc_inv_gradient, false);
+ provide(f_cvc_Jacobian, false);
+ }
+ }
+ }
+
return COLVARS_OK;
}
@@ -87,8 +100,7 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
group->key = group_key;
if (b_try_scalable) {
- // TODO rewrite this logic in terms of dependencies
- if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
+ if (is_available(f_cvc_scalable_com) && is_enabled(f_cvc_com_based)) {
enable(f_cvc_scalable_com);
enable(f_cvc_scalable);
// The CVC makes the feature available;
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index a230cae8a9a04bf9b42e4b769f8e693062b8926d..ec215cbad107485d84f534d6c0eeeaff7b552e96 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -343,6 +343,15 @@ public:
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
+ /// Redefined to override the distance ones
+ virtual cvm::real dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const;
+ /// Redefined to override the distance ones
+ virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const;
+ /// Redefined to override the distance ones
+ virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const;
};
diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp
index 4c3390a8bd1de9bf0fb37288d114d6cb8177daca..0204f3b4b1b50d83c58fe4d27be05c015cae1f10 100644
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@@ -21,7 +21,7 @@ colvar::angle::angle(std::string const &conf)
function_type = "angle";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
@@ -40,7 +40,7 @@ colvar::angle::angle(cvm::atom const &a1,
function_type = "angle";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@@ -259,7 +259,7 @@ colvar::dihedral::dihedral(std::string const &conf)
b_periodic = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
@@ -285,7 +285,7 @@ colvar::dihedral::dihedral(cvm::atom const &a1,
b_periodic = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
b_1site_force = false;
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 7481dd360ba170e327a31977d21c33d47f4507a1..f46270246f5abbfd304568c3d7a0dca6317efbd3 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -23,7 +23,7 @@ colvar::distance::distance(std::string const &conf)
function_type = "distance";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
@@ -44,7 +44,7 @@ colvar::distance::distance()
function_type = "distance";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
b_no_PBC = false;
x.type(colvarvalue::type_scalar);
}
@@ -106,7 +106,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
: distance(conf)
{
function_type = "distance_vec";
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_3vector);
}
@@ -115,7 +115,7 @@ colvar::distance_vec::distance_vec()
: distance()
{
function_type = "distance_vec";
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_3vector);
}
@@ -176,7 +176,7 @@ colvar::distance_z::distance_z(std::string const &conf)
function_type = "distance_z";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar);
// TODO detect PBC from MD engine (in simple cases)
@@ -228,7 +228,7 @@ colvar::distance_z::distance_z()
function_type = "distance_z";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar);
}
@@ -372,7 +372,7 @@ colvar::distance_xy::distance_xy(std::string const &conf)
function_type = "distance_xy";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar);
}
@@ -383,7 +383,7 @@ colvar::distance_xy::distance_xy()
function_type = "distance_xy";
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar);
}
@@ -479,7 +479,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
: distance(conf)
{
function_type = "distance_dir";
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_unit3vector);
}
@@ -488,7 +488,7 @@ colvar::distance_dir::distance_dir()
: distance()
{
function_type = "distance_dir";
- provide(f_cvc_com_based);
+ enable(f_cvc_com_based);
x.type(colvarvalue::type_unit3vector);
}
@@ -529,6 +529,27 @@ void colvar::distance_dir::apply_force(colvarvalue const &force)
}
+cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return (x1.rvector_value - x2.rvector_value).norm2();
+}
+
+
+colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
+}
+
+
+colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
+}
+
+
colvar::distance_inv::distance_inv(std::string const &conf)
: distance(conf)
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index e3ccd4ecf5815fcabc0bf7b337a06bf32340c1e7..8252f77e6270d309fdf0db3bad08644ab1c1745b 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -14,9 +14,6 @@
colvardeps::~colvardeps() {
size_t i;
- for (i=0; i<feature_states.size(); i++) {
- if (feature_states[i] != NULL) delete feature_states[i];
- }
// Do not delete features if it's static
// for (i=0; i<features.size(); i++) {
// if (features[i] != NULL) delete features[i];
@@ -34,16 +31,34 @@ colvardeps::~colvardeps() {
}
-void colvardeps::provide(int feature_id) {
- feature_states[feature_id]->available = true;
+void colvardeps::provide(int feature_id, bool truefalse) {
+ feature_states[feature_id].available = truefalse;
+}
+
+
+void colvardeps::set_enabled(int feature_id, bool truefalse) {
+// if (!is_static(feature_id)) {
+// cvm::error("Cannot set feature " + features()[feature_id]->description + " statically in " + description + ".\n");
+// return;
+// }
+ if (truefalse) {
+ // Resolve dependencies too
+ enable(feature_id);
+ } else {
+ feature_states[feature_id].enabled = false;
+ }
}
bool colvardeps::get_keyval_feature(colvarparse *cvp,
- std::string const &conf, char const *key,
- int feature_id, bool const &def_value,
- colvarparse::Parse_Mode const parse_mode)
+ std::string const &conf, char const *key,
+ int feature_id, bool const &def_value,
+ colvarparse::Parse_Mode const parse_mode)
{
+ if (!is_user(feature_id)) {
+ cvm::error("Cannot set feature " + features()[feature_id]->description + " from user input in " + description + ".\n");
+ return false;
+ }
bool value;
bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode);
if (value) enable(feature_id);
@@ -52,19 +67,19 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
int colvardeps::enable(int feature_id,
- bool dry_run /* default: false */,
- // dry_run: fail silently, do not enable if available
- // flag is passed recursively to deps of this feature
- bool toplevel /* default: true */)
- // toplevel: false if this is called as part of a chain of dependency resolution
- // this is used to diagnose failed dependencies by displaying the full stack
- // only the toplevel dependency will throw a fatal error
+ bool dry_run /* default: false */,
+ // dry_run: fail silently, do not enable if available
+ // flag is passed recursively to deps of this feature
+ bool toplevel /* default: true */)
+// toplevel: false if this is called as part of a chain of dependency resolution
+// this is used to diagnose failed dependencies by displaying the full stack
+// only the toplevel dependency will throw a fatal error
{
int res;
size_t i, j;
bool ok;
feature *f = features()[feature_id];
- feature_state *fs = feature_states[feature_id];
+ feature_state *fs = &feature_states[feature_id];
if (cvm::debug()) {
cvm::log("DEPS: " + description +
@@ -88,6 +103,14 @@ int colvardeps::enable(int feature_id,
return COLVARS_ERROR;
}
+ if (!toplevel && !is_dynamic(feature_id)) {
+ if (!dry_run) {
+ cvm::log("Non-dynamic feature : \"" + f->description
+ + "\" in " + description + " may not be enabled as a dependency.\n");
+ }
+ return COLVARS_ERROR;
+ }
+
// 1) enforce exclusions
for (i=0; i<f->requires_exclude.size(); i++) {
feature *g = features()[f->requires_exclude[i]];
@@ -168,9 +191,9 @@ int colvardeps::enable(int feature_id,
if (res != COLVARS_OK) {
if (!dry_run) {
cvm::log("...required by \"" + f->description + "\" in " + description);
- }
- if (toplevel) {
- cvm::error("Error: Failed dependency in " + description + ".");
+ if (toplevel) {
+ cvm::error("Error: Failed dependency in " + description + ".");
+ }
}
return res;
}
@@ -194,9 +217,12 @@ int colvardeps::enable(int feature_id,
// // we need refs to parents to walk up the deps tree!
// // or refresh
// }
+void colvardeps::init_feature(int feature_id, const char *description, feature_type type) {
+ features()[feature_id]->description = description;
+ features()[feature_id]->type = type;
+}
- // Shorthand macros for describing dependencies
-#define f_description(f, d) features()[f]->description = d
+// Shorthand macros for describing dependencies
#define f_req_self(f, g) features()[f]->requires_self.push_back(g)
// This macro ensures that exclusions are symmetric
#define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \
@@ -216,35 +242,31 @@ void colvardeps::init_cvb_requires() {
for (i = 0; i < f_cvb_ntot; i++) {
features().push_back(new feature);
}
- }
- f_description(f_cvb_active, "active");
- f_req_children(f_cvb_active, f_cv_active);
+ init_feature(f_cvb_active, "active", f_type_dynamic);
+ f_req_children(f_cvb_active, f_cv_active);
+
+ init_feature(f_cvb_apply_force, "apply force", f_type_user);
+ f_req_children(f_cvb_apply_force, f_cv_gradient);
- f_description(f_cvb_apply_force, "apply force");
- f_req_children(f_cvb_apply_force, f_cv_gradient);
+ init_feature(f_cvb_get_total_force, "obtain total force");
+ f_req_children(f_cvb_get_total_force, f_cv_total_force);
- f_description(f_cvb_get_total_force, "obtain total force");
- f_req_children(f_cvb_get_total_force, f_cv_total_force);
+ init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
- f_description(f_cvb_history_dependent, "history-dependent");
+ init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
+ f_req_children(f_cvb_scalar_variables, f_cv_scalar);
- f_description(f_cvb_scalar_variables, "require scalar variables");
- f_req_children(f_cvb_scalar_variables, f_cv_scalar);
+ init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
+ }
// Initialize feature_states for each instance
feature_states.reserve(f_cvb_ntot);
for (i = 0; i < f_cvb_ntot; i++) {
- feature_states.push_back(new feature_state(true, false));
+ feature_states.push_back(feature_state(true, false));
// Most features are available, so we set them so
// and list exceptions below
}
-
- // some biases are not history-dependent
- feature_states[f_cvb_history_dependent]->available = false;
-
- // by default, biases should work with vector variables, too
- feature_states[f_cvb_scalar_variables]->available = false;
}
@@ -255,117 +277,111 @@ void colvardeps::init_cv_requires() {
features().push_back(new feature);
}
- f_description(f_cv_active, "active");
+ init_feature(f_cv_active, "active", f_type_dynamic);
f_req_children(f_cv_active, f_cvc_active);
// Colvars must be either a linear combination, or scalar (and polynomial) or scripted
f_req_alt3(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted);
- f_description(f_cv_gradient, "gradient");
+ init_feature(f_cv_gradient, "gradient", f_type_dynamic);
f_req_children(f_cv_gradient, f_cvc_gradient);
- f_description(f_cv_collect_gradient, "collect gradient");
+ init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
f_req_self(f_cv_collect_gradient, f_cv_gradient);
f_req_self(f_cv_collect_gradient, f_cv_scalar);
- f_description(f_cv_fdiff_velocity, "fdiff_velocity");
+ init_feature(f_cv_fdiff_velocity, "fdiff_velocity", f_type_dynamic);
// System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
- f_description(f_cv_total_force, "total force");
+ init_feature(f_cv_total_force, "total force", f_type_dynamic);
f_req_alt2(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
// Deps for explicit total force calculation
- f_description(f_cv_total_force_calc, "total force calculation");
+ init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic);
f_req_self(f_cv_total_force_calc, f_cv_scalar);
f_req_self(f_cv_total_force_calc, f_cv_linear);
f_req_children(f_cv_total_force_calc, f_cvc_inv_gradient);
f_req_self(f_cv_total_force_calc, f_cv_Jacobian);
- f_description(f_cv_Jacobian, "Jacobian derivative");
+ init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic);
f_req_self(f_cv_Jacobian, f_cv_scalar);
f_req_self(f_cv_Jacobian, f_cv_linear);
f_req_children(f_cv_Jacobian, f_cvc_Jacobian);
- f_description(f_cv_hide_Jacobian, "hide Jacobian force");
+ init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user);
f_req_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
- f_description(f_cv_extended_Lagrangian, "extended Lagrangian");
+ init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user);
+ f_req_self(f_cv_extended_Lagrangian, f_cv_scalar);
+ f_req_self(f_cv_extended_Lagrangian, f_cv_gradient);
- f_description(f_cv_Langevin, "Langevin dynamics");
+ init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user);
f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian);
- f_description(f_cv_linear, "linear");
+ init_feature(f_cv_linear, "linear", f_type_static);
- f_description(f_cv_scalar, "scalar");
+ init_feature(f_cv_scalar, "scalar", f_type_static);
- f_description(f_cv_output_energy, "output energy");
+ init_feature(f_cv_output_energy, "output energy", f_type_user);
- f_description(f_cv_output_value, "output value");
+ init_feature(f_cv_output_value, "output value", f_type_user);
- f_description(f_cv_output_velocity, "output velocity");
+ init_feature(f_cv_output_velocity, "output velocity", f_type_user);
f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity);
- f_description(f_cv_output_applied_force, "output applied force");
+ init_feature(f_cv_output_applied_force, "output applied force", f_type_user);
- f_description(f_cv_output_total_force, "output total force");
+ init_feature(f_cv_output_total_force, "output total force", f_type_user);
f_req_self(f_cv_output_total_force, f_cv_total_force);
- f_description(f_cv_subtract_applied_force, "subtract applied force from total force");
+ init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
- f_description(f_cv_lower_boundary, "lower boundary");
+ init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
f_req_self(f_cv_lower_boundary, f_cv_scalar);
- f_description(f_cv_upper_boundary, "upper boundary");
+ init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
f_req_self(f_cv_upper_boundary, f_cv_scalar);
- f_description(f_cv_grid, "grid");
+ init_feature(f_cv_grid, "grid", f_type_user);
f_req_self(f_cv_grid, f_cv_lower_boundary);
f_req_self(f_cv_grid, f_cv_upper_boundary);
- f_description(f_cv_lower_wall, "lower wall");
- f_req_self(f_cv_lower_wall, f_cv_lower_boundary);
- f_req_self(f_cv_lower_wall, f_cv_gradient);
-
- f_description(f_cv_upper_wall, "upper wall");
- f_req_self(f_cv_upper_wall, f_cv_upper_boundary);
- f_req_self(f_cv_upper_wall, f_cv_gradient);
+ init_feature(f_cv_runave, "running average", f_type_user);
- f_description(f_cv_runave, "running average");
+ init_feature(f_cv_corrfunc, "correlation function", f_type_user);
- f_description(f_cv_corrfunc, "correlation function");
-
- // The features below are set programmatically
- f_description(f_cv_scripted, "scripted");
- f_description(f_cv_periodic, "periodic");
+ init_feature(f_cv_scripted, "scripted", f_type_static);
+ init_feature(f_cv_periodic, "periodic", f_type_static);
f_req_self(f_cv_periodic, f_cv_homogeneous);
- f_description(f_cv_scalar, "scalar");
- f_description(f_cv_linear, "linear");
- f_description(f_cv_homogeneous, "homogeneous");
+ init_feature(f_cv_scalar, "scalar", f_type_static);
+ init_feature(f_cv_linear, "linear", f_type_static);
+ init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
}
// Initialize feature_states for each instance
feature_states.reserve(f_cv_ntot);
for (i = 0; i < f_cv_ntot; i++) {
- feature_states.push_back(new feature_state(true, false));
+ feature_states.push_back(feature_state(true, false));
// Most features are available, so we set them so
// and list exceptions below
}
- // properties that may NOT be enabled as a dependency
- int unavailable_deps[] = {
- f_cv_lower_boundary,
- f_cv_upper_boundary,
- f_cv_extended_Lagrangian,
- f_cv_Langevin,
- f_cv_scripted,
- f_cv_periodic,
- f_cv_scalar,
- f_cv_linear,
- f_cv_homogeneous
- };
- for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) {
- feature_states[unavailable_deps[i]]->available = false;
- }
+// // properties that may NOT be enabled as a dependency
+// // This will be deprecated by feature types
+// int unavailable_deps[] = {
+// f_cv_lower_boundary,
+// f_cv_upper_boundary,
+// f_cv_extended_Lagrangian,
+// f_cv_Langevin,
+// f_cv_scripted,
+// f_cv_periodic,
+// f_cv_scalar,
+// f_cv_linear,
+// f_cv_homogeneous
+// };
+// for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) {
+// feature_states[unavailable_deps[i]].available = false;
+// }
}
@@ -377,34 +393,34 @@ void colvardeps::init_cvc_requires() {
features().push_back(new feature);
}
- f_description(f_cvc_active, "active");
+ init_feature(f_cvc_active, "active", f_type_dynamic);
// The dependency below may become useful if we use dynamic atom groups
// f_req_children(f_cvc_active, f_ag_active);
- f_description(f_cvc_scalar, "scalar");
+ init_feature(f_cvc_scalar, "scalar", f_type_static);
- f_description(f_cvc_gradient, "gradient");
+ init_feature(f_cvc_gradient, "gradient", f_type_dynamic);
- f_description(f_cvc_inv_gradient, "inverse gradient");
+ init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
- f_description(f_cvc_debug_gradient, "debug gradient");
+ init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
- f_description(f_cvc_Jacobian, "Jacobian derivative");
+ init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
- f_description(f_cvc_com_based, "depends on group centers of mass");
+ init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static);
// Compute total force on first site only to avoid unwanted
// coupling to other colvars (see e.g. Ciccotti et al., 2005)
- f_description(f_cvc_one_site_total_force, "compute total collective force only from one group center");
+ init_feature(f_cvc_one_site_total_force, "compute total collective force only from one group center", f_type_user);
f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
- f_description(f_cvc_scalable, "scalable calculation");
+ init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
f_req_self(f_cvc_scalable, f_cvc_scalable_com);
- f_description(f_cvc_scalable_com, "scalable calculation of centers of mass");
+ init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static);
f_req_self(f_cvc_scalable_com, f_cvc_com_based);
@@ -414,23 +430,25 @@ void colvardeps::init_cvc_requires() {
}
// Initialize feature_states for each instance
- // default as unavailable, not enabled
+ // default as available, not enabled
+ // except dynamic features which default as unavailable
feature_states.reserve(f_cvc_ntot);
for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
- feature_states.push_back(new feature_state(false, false));
+ bool avail = is_dynamic(i) ? false : true;
+ feature_states.push_back(feature_state(avail, false));
}
// Features that are implemented by all cvcs by default
// Each cvc specifies what other features are available
- feature_states[f_cvc_active]->available = true;
- feature_states[f_cvc_gradient]->available = true;
+ feature_states[f_cvc_active].available = true;
+ feature_states[f_cvc_gradient].available = true;
// Features that are implemented by default if their requirements are
- feature_states[f_cvc_one_site_total_force]->available = true;
+ feature_states[f_cvc_one_site_total_force].available = true;
// Features That are implemented only for certain simulation engine configurations
- feature_states[f_cvc_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
- feature_states[f_cvc_scalable]->available = feature_states[f_cvc_scalable_com]->available;
+ feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
+ feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available;
}
@@ -442,21 +460,21 @@ void colvardeps::init_ag_requires() {
features().push_back(new feature);
}
- f_description(f_ag_active, "active");
- f_description(f_ag_center, "translational fit");
- f_description(f_ag_rotate, "rotational fit");
- f_description(f_ag_fitting_group, "reference positions group");
- f_description(f_ag_fit_gradient_group, "fit gradient for main group");
- f_description(f_ag_fit_gradient_ref, "fit gradient for reference group");
- f_description(f_ag_atom_forces, "atomic forces");
+ init_feature(f_ag_active, "active", f_type_dynamic);
+ init_feature(f_ag_center, "translational fit", f_type_static);
+ init_feature(f_ag_rotate, "rotational fit", f_type_static);
+ init_feature(f_ag_fitting_group, "reference positions group", f_type_static);
+ init_feature(f_ag_fit_gradient_group, "fit gradient for main group", f_type_static);
+ init_feature(f_ag_fit_gradient_ref, "fit gradient for reference group", f_type_static);
+ init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
// parallel calculation implies that we have at least a scalable center of mass,
// but f_ag_scalable is kept as a separate feature to deal with future dependencies
- f_description(f_ag_scalable, "scalable group calculation");
- f_description(f_ag_scalable_com, "scalable group center of mass calculation");
+ init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
+ init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static);
f_req_self(f_ag_scalable, f_ag_scalable_com);
-// f_description(f_ag_min_msd_fit, "minimum MSD fit")
+// init_feature(f_ag_min_msd_fit, "minimum MSD fit")
// f_req_self(f_ag_min_msd_fit, f_ag_center)
// f_req_self(f_ag_min_msd_fit, f_ag_rotate)
// f_req_exclude(f_ag_min_msd_fit, f_ag_fitting_group)
@@ -466,15 +484,15 @@ void colvardeps::init_ag_requires() {
// default as unavailable, not enabled
feature_states.reserve(f_ag_ntot);
for (i = 0; i < colvardeps::f_ag_ntot; i++) {
- feature_states.push_back(new feature_state(false, false));
+ feature_states.push_back(feature_state(false, false));
}
// Features that are implemented (or not) by all atom groups
- feature_states[f_ag_active]->available = true;
+ feature_states[f_ag_active].available = true;
// f_ag_scalable_com is provided by the CVC iff it is COM-based
- feature_states[f_ag_scalable_com]->available = false;
+ feature_states[f_ag_scalable_com].available = false;
// TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
- feature_states[f_ag_scalable]->available = true;
+ feature_states[f_ag_scalable].available = true;
}
@@ -482,7 +500,7 @@ void colvardeps::print_state() {
size_t i;
cvm::log("Enabled features of " + description);
for (i = 0; i < feature_states.size(); i++) {
- if (feature_states[i]->enabled)
+ if (feature_states[i].enabled)
cvm::log("- " + features()[i]->description);
}
for (i=0; i<children.size(); i++) {
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index 4ef27ded85a5b2a1bbcaaad54933a9294e11c51a..fd07cb645715eac8a8e367d654b6fbd06735a633 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -13,17 +13,16 @@
#include "colvarmodule.h"
#include "colvarparse.h"
-/// Parent class for a member object of a bias, cv or cvc etc. containing dependencies
-/// (features) and handling dependency resolution
-
-// Some features like colvar::f_linear have no dependencies, require() doesn't enable anything but fails if unavailable
-// Policy: those features are unavailable at all times
-// Other features are under user control
-// They are unavailable unless requested by the user, then they may be enabled subject to
-// satisfied dependencies
-
-// It seems important to have available default to false (for safety) and enabled to false (for efficiency)
-
+/// \brief Parent class for a member object of a bias, cv or cvc etc. containing features and
+/// their dependencies, and handling dependency resolution
+///
+/// There are 3 kinds of features:
+/// 1. Dynamic features are under the control of the dependency resolution
+/// system. They may be enabled or disabled depending on dependencies.
+/// 2. User features may be enabled based on user input (they may trigger a failure upon dependency resolution, though)
+/// 3. Static features are static properties of the object, determined
+/// programatically at initialization time.
+///
class colvardeps {
public:
@@ -39,12 +38,7 @@ public:
feature_state(bool a, bool e)
: available(a), enabled(e) {}
- /// Available means: supported, subject to dependencies as listed,
- /// MAY BE ENABLED AS A RESULT OF DEPENDENCY SOLVING
- /// Remains false for passive flags that are set based on other object properties,
- /// eg. f_cv_linear
- /// Is set to true upon user request for features that are implemented by the user
- /// or under his/her direct control, e.g. f_cv_scripted or f_cv_extended_Lagrangian
+ /// Feature may be enabled, subject to possible dependencies
bool available;
/// Currently enabled - this flag is subject to change dynamically
/// TODO consider implications for dependency solving: anyone who disables
@@ -53,8 +47,22 @@ public:
// bool enabledOnce; // this should trigger an update when object is evaluated
};
- /// List of the state of all features
- std::vector<feature_state *> feature_states;
+
+private:
+ /// List of the states of all features
+ std::vector<feature_state> feature_states;
+
+ /// Enum of possible feature types
+ enum feature_type {
+ f_type_not_set,
+ f_type_dynamic,
+ f_type_user,
+ f_type_static
+ };
+
+public:
+ /// Pair a numerical feature ID with a description and type
+ void init_feature(int feature_id, const char *description, feature_type type = f_type_not_set);
/// Describes a feature and its dependecies
/// used in a static array within each subclass
@@ -83,8 +91,18 @@ public:
// features that this feature requires in children
std::vector<int> requires_children;
+
+ inline bool is_dynamic() { return type == f_type_dynamic; }
+ inline bool is_static() { return type == f_type_static; }
+ inline bool is_user() { return type == f_type_user; }
+ /// Type of this feature, from the enum feature_type
+ feature_type type;
};
+ inline bool is_dynamic(int id) { return features()[id]->type == f_type_dynamic; }
+ inline bool is_static(int id) { return features()[id]->type == f_type_static; }
+ inline bool is_user(int id) { return features()[id]->type == f_type_user; }
+
// Accessor to array of all features with deps, static in most derived classes
// Subclasses with dynamic dependency trees may override this
// with a non-static array
@@ -100,9 +118,8 @@ public:
/// (useful for cvcs because their children are member objects)
void remove_all_children();
-
-
private:
+
// pointers to objects this object depends on
// list should be maintained by any code that modifies the object
// this could be secured by making lists of colvars / cvcs / atom groups private and modified through accessor functions
@@ -130,16 +147,28 @@ public:
// Checks whether given feature is enabled
// Defaults to querying f_*_active
inline bool is_enabled(int f = f_cv_active) const {
- return feature_states[f]->enabled;
+ return feature_states[f].enabled;
}
// Checks whether given feature is available
// Defaults to querying f_*_active
inline bool is_available(int f = f_cv_active) const {
- return feature_states[f]->available;
+ return feature_states[f].available;
}
- void provide(int feature_id); // set the feature's flag to available in local object
+ /// Set the feature's available flag, without checking
+ /// To be used for dynamic properties
+ /// dependencies will be checked by enable()
+ void provide(int feature_id, bool truefalse = true);
+
+ /// Set the feature's enabled flag, without dependency check or resolution
+ /// To be used for static properties only
+ /// Checking for availability is up to the caller
+ void set_enabled(int feature_id, bool truefalse = true);
+
+protected:
+
+
/// Parse a keyword and enable a feature accordingly
bool get_keyval_feature(colvarparse *cvp,
@@ -147,6 +176,8 @@ public:
int feature_id, bool const &def_value,
colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
+public:
+
int enable(int f, bool dry_run = false, bool toplevel = true); // enable a feature and recursively solve its dependencies
// dry_run is set to true to recursively test if a feature is available, without enabling it
// int disable(int f);
@@ -165,6 +196,7 @@ public:
f_cvb_get_total_force, // requires total forces
f_cvb_history_dependent, // depends on simulation history
f_cvb_scalar_variables, // requires scalar colvars
+ f_cvb_calc_pmf, // whether this bias will compute a PMF
f_cvb_ntot
};
@@ -216,12 +248,6 @@ public:
/// be used by the biases or in analysis (needs lower and upper
/// boundary)
f_cv_grid,
- /// \brief Apply a restraining potential (|x-xb|^2) to the colvar
- /// when it goes below the lower wall
- f_cv_lower_wall,
- /// \brief Apply a restraining potential (|x-xb|^2) to the colvar
- /// when it goes above the upper wall
- f_cv_upper_wall,
/// \brief Compute running average
f_cv_runave,
/// \brief Compute time correlation function
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index b9a435152b2cc56b50987b08adb7cf2bd75e4362..10cd3c0e4737905dabb664bcaa0d159608fd16f7 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -36,6 +36,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
"variable module twice.\n");
return;
}
+
+ depth_s = 0;
+
cvm::log(cvm::line_marker);
cvm::log("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n");
@@ -71,6 +74,48 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
}
+colvarmodule * colvarmodule::main()
+{
+ return proxy->colvars;
+}
+
+
+std::vector<colvar *> *colvarmodule::variables()
+{
+ return &colvars;
+}
+
+
+std::vector<colvar *> *colvarmodule::variables_active()
+{
+ return &colvars_active;
+}
+
+
+std::vector<colvar *> *colvarmodule::variables_active_smp()
+{
+ return &colvars_smp;
+}
+
+
+std::vector<int> *colvarmodule::variables_active_smp_items()
+{
+ return &colvars_smp_items;
+}
+
+
+std::vector<colvarbias *> *colvarmodule::biases_active()
+{
+ return &(biases_active_);
+}
+
+
+size_t colvarmodule::size() const
+{
+ return colvars.size() + biases.size();
+}
+
+
int colvarmodule::read_config_file(char const *config_filename)
{
cvm::log(cvm::line_marker);
@@ -118,8 +163,29 @@ int colvarmodule::read_config_string(std::string const &config_str)
}
+std::istream & colvarmodule::getline(std::istream &is, std::string &line)
+{
+ std::string l;
+ if (std::getline(is, l)) {
+ size_t const sz = l.size();
+ if (sz > 0) {
+ if (l[sz-1] == '\r' ) {
+ line = l.substr(0, sz-1);
+ } else {
+ line = l;
+ }
+ } else {
+ line.clear();
+ }
+ }
+ return is;
+}
+
+
int colvarmodule::parse_config(std::string &conf)
{
+ extra_conf.clear();
+
// parse global options
if (catch_input_errors(parse_global_params(conf))) {
return get_error();
@@ -140,6 +206,15 @@ int colvarmodule::parse_config(std::string &conf)
return get_error();
}
+ if (extra_conf.size()) {
+ catch_input_errors(parse_global_params(extra_conf));
+ catch_input_errors(parse_colvars(extra_conf));
+ catch_input_errors(parse_biases(extra_conf));
+ parse->check_keywords(extra_conf, "colvarmodule");
+ extra_conf.clear();
+ if (get_error() != COLVARS_OK) return get_error();
+ }
+
cvm::log(cvm::line_marker);
cvm::log("Collective variables module (re)initialized.\n");
cvm::log(cvm::line_marker);
@@ -156,11 +231,20 @@ int colvarmodule::parse_config(std::string &conf)
}
+int colvarmodule::append_new_config(std::string const &new_conf)
+{
+ extra_conf += new_conf;
+ return COLVARS_OK;
+}
+
+
int colvarmodule::parse_global_params(std::string const &conf)
{
+ colvarmodule *cvm = cvm::main();
+
std::string index_file_name;
if (parse->get_keyval(conf, "indexFile", index_file_name)) {
- read_index_file(index_file_name.c_str());
+ cvm->read_index_file(index_file_name.c_str());
}
if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
@@ -216,8 +300,8 @@ int colvarmodule::parse_colvars(std::string const &conf)
if (colvar_conf.size()) {
cvm::log(cvm::line_marker);
cvm::increase_depth();
- colvars.push_back(new colvar(colvar_conf));
- if (cvm::get_error() ||
+ colvars.push_back(new colvar());
+ if (((colvars.back())->init(colvar_conf) != COLVARS_OK) ||
((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
cvm::log("Error while constructing colvar number " +
cvm::to_str(colvars.size()) + " : deleting.");
@@ -262,8 +346,7 @@ bool colvarmodule::check_new_bias(std::string &conf, char const *key)
template <class bias_type>
int colvarmodule::parse_biases_type(std::string const &conf,
- char const *keyword,
- size_t &bias_count)
+ char const *keyword)
{
std::string bias_conf = "";
size_t conf_saved_pos = 0;
@@ -277,7 +360,6 @@ int colvarmodule::parse_biases_type(std::string const &conf,
return COLVARS_ERROR;
}
cvm::decrease_depth();
- bias_count++;
} else {
cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n",
INPUT_ERROR);
@@ -295,28 +377,28 @@ int colvarmodule::parse_biases(std::string const &conf)
cvm::log("Initializing the collective variables biases.\n");
/// initialize ABF instances
- parse_biases_type<colvarbias_abf>(conf, "abf", n_abf_biases);
+ parse_biases_type<colvarbias_abf>(conf, "abf");
/// initialize adaptive linear biases
- parse_biases_type<colvarbias_alb>(conf, "ALB", n_rest_biases);
+ parse_biases_type<colvarbias_alb>(conf, "ALB");
/// initialize harmonic restraints
- parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);
+ parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic");
/// initialize harmonic walls restraints
- parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls", n_rest_biases);
+ parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls");
/// initialize histograms
- parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
+ parse_biases_type<colvarbias_histogram>(conf, "histogram");
/// initialize histogram restraints
- parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint", n_rest_biases);
+ parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint");
/// initialize linear restraints
- parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);
+ parse_biases_type<colvarbias_restraint_linear>(conf, "linear");
/// initialize metadynamics instances
- parse_biases_type<colvarbias_meta>(conf, "metadynamics", n_meta_biases);
+ parse_biases_type<colvarbias_meta>(conf, "metadynamics");
if (use_scripted_forces) {
cvm::log(cvm::line_marker);
@@ -363,6 +445,36 @@ int colvarmodule::parse_biases(std::string const &conf)
}
+int colvarmodule::num_biases_feature(int feature_id)
+{
+ colvarmodule *cv = cvm::main();
+ size_t n = 0;
+ for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
+ bi != cv->biases.end();
+ bi++) {
+ if ((*bi)->is_enabled(feature_id)) {
+ n++;
+ }
+ }
+ return n;
+}
+
+
+int colvarmodule::num_biases_type(std::string const &type)
+{
+ colvarmodule *cv = cvm::main();
+ size_t n = 0;
+ for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
+ bi != cv->biases.end();
+ bi++) {
+ if ((*bi)->bias_type == type) {
+ n++;
+ }
+ }
+ return n;
+}
+
+
int colvarmodule::catch_input_errors(int result)
{
if (result != COLVARS_OK || get_error()) {
@@ -376,8 +488,9 @@ int colvarmodule::catch_input_errors(int result)
colvarbias * colvarmodule::bias_by_name(std::string const &name) {
- for (std::vector<colvarbias *>::iterator bi = biases.begin();
- bi != biases.end();
+ colvarmodule *cv = cvm::main();
+ for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
+ bi != cv->biases.end();
bi++) {
if ((*bi)->name == name) {
return (*bi);
@@ -388,8 +501,9 @@ colvarbias * colvarmodule::bias_by_name(std::string const &name) {
colvar *colvarmodule::colvar_by_name(std::string const &name) {
- for (std::vector<colvar *>::iterator cvi = colvars.begin();
- cvi != colvars.end();
+ colvarmodule *cv = cvm::main();
+ for (std::vector<colvar *>::iterator cvi = cv->colvars.begin();
+ cvi != cv->colvars.end();
cvi++) {
if ((*cvi)->name == name) {
return (*cvi);
@@ -555,9 +669,15 @@ int colvarmodule::calc_colvars()
int error_code = COLVARS_OK;
std::vector<colvar *>::iterator cvi;
- // Determine which colvars are active at this time step
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
- (*cvi)->feature_states[colvardeps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 0);
+ // Determine which colvars are active at this iteration
+ variables_active()->resize(0);
+ variables_active()->reserve(variables()->size());
+ for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
+ // This is a dynamic feature - the next call should be to enable()
+ // or disable() when dynamic dependency resolution is fully implemented
+ (*cvi)->set_enabled(colvardeps::f_cv_active,
+ step_absolute() % (*cvi)->get_time_step_factor() == 0);
+ variables_active()->push_back(*cvi);
}
// if SMP support is available, split up the work
@@ -565,25 +685,24 @@ int colvarmodule::calc_colvars()
// first, calculate how much work (currently, how many active CVCs) each colvar has
- colvars_smp.resize(0);
- colvars_smp_items.resize(0);
+ variables_active_smp()->resize(0);
+ variables_active_smp_items()->resize(0);
- colvars_smp.reserve(colvars.size());
- colvars_smp_items.reserve(colvars.size());
+ variables_active_smp()->reserve(variables_active()->size());
+ variables_active_smp_items()->reserve(variables_active()->size());
// set up a vector containing all components
cvm::increase_depth();
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+ for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
- if (!(*cvi)->is_enabled()) continue;
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);
+ variables_active_smp()->reserve(variables_active_smp()->size() + num_items);
+ variables_active_smp_items()->reserve(variables_active_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);
+ variables_active_smp()->push_back(*cvi);
+ variables_active_smp_items()->push_back(icvc);
}
}
cvm::decrease_depth();
@@ -592,8 +711,7 @@ int colvarmodule::calc_colvars()
error_code |= proxy->smp_colvars_loop();
cvm::increase_depth();
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
- if (!(*cvi)->is_enabled()) continue;
+ for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
error_code |= (*cvi)->collect_cvc_data();
}
cvm::decrease_depth();
@@ -602,8 +720,7 @@ int colvarmodule::calc_colvars()
// calculate colvars one at a time
cvm::increase_depth();
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
- if (!(*cvi)->is_enabled()) continue;
+ for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
error_code |= (*cvi)->calc();
if (cvm::get_error()) {
return COLVARS_ERROR;
@@ -627,6 +744,18 @@ int colvarmodule::calc_biases()
std::vector<colvarbias *>::iterator bi;
int error_code = COLVARS_OK;
+ // Total bias energy is reset before calling scripted biases
+ total_bias_energy = 0.0;
+
+ // update the list of active biases
+ biases_active()->resize(0);
+ biases_active()->reserve(biases.size());
+ for (bi = biases.begin(); bi != biases.end(); bi++) {
+ if ((*bi)->is_enabled()) {
+ biases_active()->push_back(*bi);
+ }
+ }
+
// if SMP support is available, split up the work
if (proxy->smp_enabled() == COLVARS_OK) {
@@ -645,7 +774,7 @@ int colvarmodule::calc_biases()
}
cvm::increase_depth();
- for (bi = biases.begin(); bi != biases.end(); bi++) {
+ for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
error_code |= (*bi)->update();
if (cvm::get_error()) {
return error_code;
@@ -654,12 +783,10 @@ int colvarmodule::calc_biases()
cvm::decrease_depth();
}
- cvm::real total_bias_energy = 0.0;
- for (bi = biases.begin(); bi != biases.end(); bi++) {
+ for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
total_bias_energy += (*bi)->get_energy();
}
- proxy->add_energy(total_bias_energy);
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@@ -675,7 +802,7 @@ int colvarmodule::update_colvar_forces()
if (cvm::debug() && biases.size())
cvm::log("Collecting forces from all biases.\n");
cvm::increase_depth();
- for (bi = biases.begin(); bi != biases.end(); bi++) {
+ for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
(*bi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
@@ -687,6 +814,11 @@ int colvarmodule::update_colvar_forces()
error_code |= calc_scripted_forces();
}
+ // Now we have collected energies from both built-in and scripted biases
+ if (cvm::debug())
+ cvm::log("Adding total bias energy: " + cvm::to_str(total_bias_energy) + "\n");
+ proxy->add_energy(total_bias_energy);
+
cvm::real total_colvar_energy = 0.0;
// sum up the forces for each colvar, including wall forces
// and integrate any internal
@@ -695,7 +827,7 @@ int colvarmodule::update_colvar_forces()
cvm::log("Updating the internal degrees of freedom "
"of colvars (if they have any).\n");
cvm::increase_depth();
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+ for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
// Here we call even inactive colvars, so they accumulate biasing forces
// as well as update their extended-system dynamics
total_colvar_energy += (*cvi)->update_forces_energy();
@@ -704,6 +836,8 @@ int colvarmodule::update_colvar_forces()
}
}
cvm::decrease_depth();
+ if (cvm::debug())
+ cvm::log("Adding total colvar energy: " + cvm::to_str(total_colvar_energy) + "\n");
proxy->add_energy(total_colvar_energy);
// make collective variables communicate their forces to their
@@ -711,9 +845,8 @@ int colvarmodule::update_colvar_forces()
if (cvm::debug())
cvm::log("Communicating forces from the colvars to the atoms.\n");
cvm::increase_depth();
- for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+ for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) {
- if (!(*cvi)->is_enabled()) continue;
(*cvi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
@@ -800,8 +933,8 @@ int colvarmodule::analyze()
cvm::log("Performing analysis.\n");
// perform colvar-specific analysis
- for (std::vector<colvar *>::iterator cvi = colvars.begin();
- cvi != colvars.end();
+ for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
+ cvi != variables_active()->end();
cvi++) {
cvm::increase_depth();
(*cvi)->analyze();
@@ -825,8 +958,8 @@ int colvarmodule::setup()
{
if (this->size() == 0) return cvm::get_error();
// loop over all components of all colvars to reset masses of all groups
- for (std::vector<colvar *>::iterator cvi = colvars.begin();
- cvi != colvars.end(); cvi++) {
+ for (std::vector<colvar *>::iterator cvi = variables()->begin();
+ cvi != variables()->end(); cvi++) {
(*cvi)->setup();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@@ -934,18 +1067,18 @@ int colvarmodule::setup_output()
cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
}
- output_prefix = proxy->output_prefix();
- if (output_prefix.size()) {
+ output_prefix() = proxy->output_prefix();
+ if (output_prefix().size()) {
cvm::log("The final output state file will be \""+
- (output_prefix.size() ?
- std::string(output_prefix+".colvars.state") :
+ (output_prefix().size() ?
+ std::string(output_prefix()+".colvars.state") :
std::string("colvars.state"))+"\".\n");
// cvm::log (cvm::line_marker);
}
cv_traj_name =
- (output_prefix.size() ?
- std::string(output_prefix+".colvars.traj") :
+ (output_prefix().size() ?
+ std::string(output_prefix()+".colvars.traj") :
std::string(""));
if (cv_traj_freq && cv_traj_name.size()) {
@@ -1077,14 +1210,14 @@ continue the previous simulation.\n\n");
cvm::log(cvm::line_marker);
// update this ahead of time in this special case
- output_prefix = proxy->input_prefix();
- cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n");
+ output_prefix() = proxy->input_prefix();
+ cvm::log("All output files will now be saved with the prefix \""+output_prefix()+".tmp.*\".\n");
cvm::log(cvm::line_marker);
cvm::log("Please review the important warning above. After that, you may rename:\n\
-\""+output_prefix+".tmp.colvars.state\"\n\
+\""+output_prefix()+".tmp.colvars.state\"\n\
to:\n\
\""+ proxy->input_prefix()+".colvars.state\"\n");
- output_prefix = output_prefix+".tmp";
+ output_prefix() = output_prefix()+".tmp";
write_output_files();
cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
}
@@ -1104,8 +1237,8 @@ int colvarmodule::write_output_files()
// if this is a simulation run (i.e. not a postprocessing), output data
// must be written to be able to restart the simulation
std::string const out_name =
- (output_prefix.size() ?
- std::string(output_prefix+".colvars.state") :
+ (output_prefix().size() ?
+ std::string(output_prefix()+".colvars.state") :
std::string("colvars.state"));
cvm::log("Saving collective variables state to \""+out_name+"\".\n");
@@ -1340,7 +1473,8 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
void cvm::log(std::string const &message)
{
- size_t const d = depth();
+ // allow logging when the module is not fully initialized
+ size_t const d = (cvm::main() != NULL) ? depth() : 0;
if (d > 0)
proxy->log((std::string(2*d, ' '))+message);
else
@@ -1365,19 +1499,20 @@ void cvm::decrease_depth()
size_t & cvm::depth()
{
// NOTE: do not call log() or error() here, to avoid recursion
- size_t const nt = proxy->smp_num_threads();
+ colvarmodule *cv = cvm::main();
if (proxy->smp_enabled() == COLVARS_OK) {
- if (depth_v.size() != nt) {
- // update array of depths
+ int const nt = proxy->smp_num_threads();
+ if (int(cv->depth_v.size()) != nt) {
proxy->smp_lock();
- if (depth_v.size() > 0) { depth_s = depth_v[0]; }
- depth_v.clear();
- depth_v.assign(nt, depth_s);
+ // update array of depths
+ if (cv->depth_v.size() > 0) { cv->depth_s = cv->depth_v[0]; }
+ cv->depth_v.clear();
+ cv->depth_v.assign(nt, cv->depth_s);
proxy->smp_unlock();
}
- return depth_v[proxy->smp_thread_id()];
+ return cv->depth_v[proxy->smp_thread_id()];
}
- return depth_s;
+ return cv->depth_s;
}
@@ -1451,8 +1586,8 @@ int cvm::read_index_file(char const *filename)
FILE_ERROR);
}
}
- cvm::index_group_names.push_back(group_name);
- cvm::index_groups.push_back(std::vector<int> ());
+ index_group_names.push_back(group_name);
+ index_groups.push_back(std::vector<int>());
} else {
cvm::error("Error: in parsing index file \""+
std::string(filename)+"\".\n",
@@ -1462,7 +1597,7 @@ int cvm::read_index_file(char const *filename)
int atom_number = 1;
size_t pos = is.tellg();
while ( (is >> atom_number) && (atom_number > 0) ) {
- (cvm::index_groups.back()).push_back(atom_number);
+ (index_groups.back()).push_back(atom_number);
pos = is.tellg();
}
is.clear();
@@ -1532,8 +1667,8 @@ int cvm::load_coords_xyz(char const *filename,
+ std::string(filename) + ".\n", INPUT_ERROR);
}
// skip comment line
- std::getline(xyz_is, line);
- std::getline(xyz_is, line);
+ cvm::getline(xyz_is, line);
+ cvm::getline(xyz_is, line);
xyz_is.width(255);
std::vector<atom_pos>::iterator pos_i = pos.begin();
@@ -1543,7 +1678,7 @@ int cvm::load_coords_xyz(char const *filename,
for ( ; pos_i != pos.end() ; pos_i++, index++) {
while ( next < *index ) {
- std::getline(xyz_is, line);
+ cvm::getline(xyz_is, line);
next++;
}
xyz_is >> symbol;
@@ -1558,15 +1693,9 @@ int cvm::load_coords_xyz(char const *filename,
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
-// static pointers
-std::vector<colvar *> colvarmodule::colvars;
-std::vector<colvarbias *> colvarmodule::biases;
-size_t colvarmodule::n_abf_biases = 0;
-size_t colvarmodule::n_rest_biases = 0;
-size_t colvarmodule::n_histo_biases = 0;
-size_t colvarmodule::n_meta_biases = 0;
-colvarproxy *colvarmodule::proxy = NULL;
+// shared pointer to the proxy object
+colvarproxy *colvarmodule::proxy = NULL;
// static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
@@ -1575,14 +1704,9 @@ 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_s = 0;
-std::vector<size_t> colvarmodule::depth_v(0);
bool colvarmodule::b_analysis = false;
-std::list<std::string> colvarmodule::index_group_names;
-std::list<std::vector<int> > colvarmodule::index_groups;
-bool colvarmodule::use_scripted_forces = false;
-bool colvarmodule::scripting_after_biases = true;
-std::string colvarmodule::output_prefix = "";
+bool colvarmodule::use_scripted_forces = false;
+bool colvarmodule::scripting_after_biases = true;
// i/o constants
size_t const colvarmodule::it_width = 12;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 0e98709adb79c7e746d76ed3fa331a65c2809e65..b4f80e0b5e47bcd0dd4475ebe9a22ddc90981ed5 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -11,7 +11,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-12-27"
+#define COLVARS_VERSION "2017-03-09"
#endif
#ifndef COLVARS_DEBUG
@@ -161,12 +161,37 @@ public:
/// dt)
static real debug_gradients_step_size;
+private:
+
/// Prefix for all output files for this run
- static std::string output_prefix;
+ std::string cvm_output_prefix;
+
+public:
+ /// Accessor for the above
+ static inline std::string &output_prefix()
+ {
+ colvarmodule *cv = colvarmodule::main();
+ return cv->cvm_output_prefix;
+ }
+private:
+
+ /// Array of collective variables
+ std::vector<colvar *> colvars;
/// Array of collective variables
- static std::vector<colvar *> colvars;
+ std::vector<colvar *> colvars_active;
+
+ /// Collective variables to be calculated on different threads;
+ /// colvars with multple items (e.g. multiple active CVCs) are duplicated
+ std::vector<colvar *> colvars_smp;
+ /// Indexes of the items to calculate for each colvar
+ std::vector<int> colvars_smp_items;
+
+public:
+
+ /// Array of collective variables
+ std::vector<colvar *> *variables();
/* TODO: implement named CVCs
/// Array of named (reusable) collective variable components
@@ -177,26 +202,31 @@ public:
}
*/
+ /// Collective variables with the active flag on
+ std::vector<colvar *> *variables_active();
+
/// Collective variables to be calculated on different threads;
/// colvars with multple items (e.g. multiple active CVCs) are duplicated
- std::vector<colvar *> colvars_smp;
+ std::vector<colvar *> *variables_active_smp();
+
/// Indexes of the items to calculate for each colvar
- std::vector<int> colvars_smp_items;
+ std::vector<int> *variables_active_smp_items();
/// Array of collective variable biases
- static std::vector<colvarbias *> biases;
- /// \brief Number of ABF biases initialized (in normal conditions
- /// should be 1)
- static size_t n_abf_biases;
- /// \brief Number of metadynamics biases initialized (in normal
- /// conditions should be 1)
- static size_t n_meta_biases;
- /// \brief Number of restraint biases initialized (no limit on the
- /// number)
- static size_t n_rest_biases;
- /// \brief Number of histograms initialized (no limit on the
- /// number)
- static size_t n_histo_biases;
+ std::vector<colvarbias *> biases;
+
+ /// Energy of built-in and scripted biases, summed per time-step
+ real total_bias_energy;
+
+private:
+
+ /// Array of active collective variable biases
+ std::vector<colvarbias *> biases_active_;
+
+public:
+
+ /// Array of active collective variable biases
+ std::vector<colvarbias *> *biases_active();
/// \brief Whether debug output should be enabled (compile-time option)
static inline bool debug()
@@ -205,10 +235,7 @@ public:
}
/// \brief How many objects are configured yet?
- inline size_t size() const
- {
- return colvars.size() + biases.size();
- }
+ size_t size() const;
/// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name
@@ -230,9 +257,12 @@ public:
/// \brief Parse a "clean" config string (no comments)
int parse_config(std::string &conf);
-
// Parse functions (setup internal data based on a string)
+ /// Allow reading from Windows text files using using std::getline
+ /// (which can still be used when the text is produced by Colvars itself)
+ static std::istream & getline(std::istream &is, std::string &line);
+
/// Parse the few module's global parameters
int parse_global_params(std::string const &conf);
@@ -242,14 +272,33 @@ public:
/// Parse and initialize collective variable biases
int parse_biases(std::string const &conf);
+ /// \brief Add new configuration during parsing (e.g. to implement
+ /// back-compatibility); cannot be nested, i.e. conf should not contain
+ /// anything that triggers another call
+ int append_new_config(std::string const &conf);
+
+private:
+
+ /// Auto-generated configuration during parsing (e.g. to implement
+ /// back-compatibility)
+ std::string extra_conf;
+
/// Parse and initialize collective variable biases of a specific type
template <class bias_type>
- int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count);
+ int parse_biases_type(std::string const &conf, char const *keyword);
/// Test error condition and keyword parsing
/// on error, delete new bias
bool check_new_bias(std::string &conf, char const *key);
+public:
+
+ /// Return how many biases have this feature enabled
+ static int num_biases_feature(int feature_id);
+
+ /// Return how many biases are defined with this type
+ static int num_biases_type(std::string const &type);
+
private:
/// Useful wrapper to interrupt parsing if any error occurs
int catch_input_errors(int result);
@@ -449,13 +498,13 @@ public:
/// \brief Names of groups from a Gromacs .ndx file to be read at startup
- static std::list<std::string> index_group_names;
+ std::list<std::string> index_group_names;
/// \brief Groups from a Gromacs .ndx file read at startup
- static std::list<std::vector<int> > index_groups;
+ std::list<std::vector<int> > index_groups;
/// \brief Read a Gromacs .ndx file
- static int read_index_file(char const *filename);
+ int read_index_file(char const *filename);
/// \brief Create atoms from a file \param filename name of the file
@@ -515,13 +564,13 @@ protected:
/// Output restart file
colvarmodule::ofstream restart_out_os;
-protected:
+private:
/// Counter for the current depth in the object hierarchy (useg e.g. in output)
- static size_t depth_s;
+ size_t depth_s;
/// Thread-specific depth
- static std::vector<size_t> depth_v;
+ std::vector<size_t> depth_v;
public:
@@ -552,6 +601,10 @@ public:
/// from the hosting program; it is static in order to be accessible
/// from static functions in the colvarmodule class
static colvarproxy *proxy;
+
+ /// \brief Accessor for the above
+ static colvarmodule *main();
+
};
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index 4d2e17f098dbb55457e5b1d8bc44f2cab79fdf70..8055d925dbbd3a403649c091563eb6e4e5e1cd0e 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -503,7 +503,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
std::string line;
std::istringstream is(conf);
- while (std::getline(is, line)) {
+ while (cvm::getline(is, line)) {
if (line.size() == 0)
continue;
if (line.find_first_not_of(white_space) ==
@@ -539,10 +539,9 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
std::istream & colvarparse::getline_nocomments(std::istream &is,
- std::string &line,
- char const delim)
+ std::string &line)
{
- std::getline(is, line, delim);
+ cvm::getline(is, line);
size_t const comment = line.find('#');
if (comment != std::string::npos) {
line.erase(comment);
diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h
index acdc46ddc9ee42d39f5889f2ebf4ab01c1781156..9f116caafdf8a5b96ec19ccc21469d2abc60741c 100644
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@@ -304,8 +304,7 @@ public:
/// \brief Works as std::getline() but also removes everything
/// between a comment character and the following newline
static std::istream & getline_nocomments(std::istream &is,
- std::string &s,
- char const delim = '\n');
+ std::string &s);
/// Check if the content of the file has matching braces
bool brace_check(std::string const &conf,
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index bdf65c0b1f0a617a09427769e151a6db688dd154..5b216c9d41528f390f3f99dfd86184cd847e1186 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -336,8 +336,6 @@ public:
/// Are total forces being used?
virtual bool total_forces_enabled() const
{
- cvm::error("Error: total forces are currently not implemented.\n",
- COLVARS_NOT_IMPLEMENTED);
return false;
}
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index ffb83e97801df0d4bb448b75747925e729fdb4ff..f192dcb7c05e377ba2dc44e8cd7ebf6e98481f73 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -12,6 +12,7 @@
#include <string.h>
#include "colvarscript.h"
+#include "colvardeps.h"
colvarscript::colvarscript(colvarproxy *p)
@@ -221,6 +222,16 @@ int colvarscript::run(int argc, char const *argv[]) {
}
}
+ if (cmd == "addenergy") {
+ if (argc == 3) {
+ colvars->total_bias_energy += strtod(argv[2], NULL);
+ return COLVARS_OK;
+ } else {
+ result = "Wrong arguments to command \"addenergy\"\n" + help_string();
+ return COLVARSCRIPT_ERROR;
+ }
+ }
+
result = "Syntax error\n" + help_string();
return COLVARSCRIPT_ERROR;
}
@@ -263,10 +274,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
}
if (subcmd == "delete") {
- if (cv->biases.size() > 0) {
- result = "Cannot delete a colvar currently used by biases, delete those biases first";
- return COLVARSCRIPT_ERROR;
+ size_t i;
+ for (i = 0; i < cv->biases.size(); i++) {
+ delete cv->biases[i];
}
+ cv->biases.resize(0);
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
@@ -338,9 +350,8 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
return COLVARS_OK;
}
- if (subcmd == "state") {
- cv->print_state();
- return COLVARS_OK;
+ if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
+ return proc_features(cv, argc, argv);
}
result = "Syntax error\n" + help_string();
@@ -379,11 +390,6 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
return COLVARS_OK;
}
- if (subcmd == "state") {
- b->print_state();
- return COLVARS_OK;
- }
-
// Subcommands for MW ABF
if (subcmd == "bin") {
int r = b->current_bin();
@@ -420,6 +426,10 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
return COLVARS_OK;
}
+ if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
+ return proc_features(b, argc, argv);
+ }
+
if (argc >= 4) {
std::string param = argv[3];
if (subcmd == "count") {
@@ -441,6 +451,83 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
}
+int colvarscript::proc_features(colvardeps *obj,
+ int argc, char const *argv[]) {
+ // size was already checked before calling
+ std::string subcmd = argv[2];
+
+ if (argc == 3) {
+ if (subcmd == "state") {
+ // TODO make this returned as result?
+ obj->print_state();
+ return COLVARS_OK;
+ }
+
+ // get and set commands require more arguments
+ result = "Syntax error\n" + help_string();
+ return COLVARSCRIPT_ERROR;
+ }
+
+ if ((subcmd == "get") || (subcmd == "set")) {
+ std::vector<colvardeps::feature *> &features = obj->features();
+ std::string const req_feature(argv[3]);
+ colvardeps::feature *f = NULL;
+ int fid = 0;
+ for (fid = 0; fid < int(features.size()); fid++) {
+ if (features[fid]->description ==
+ colvarparse::to_lower_cppstr(req_feature)) {
+ f = features[fid];
+ break;
+ }
+ }
+
+ if (f == NULL) {
+
+ result = "Error: feature \""+req_feature+"\" does not exist.\n";
+ return COLVARSCRIPT_ERROR;
+
+ } else {
+
+ if (! obj->is_available(fid)) {
+ result = "Error: feature \""+req_feature+"\" is unavailable.\n";
+ return COLVARSCRIPT_ERROR;
+ }
+
+ if (subcmd == "get") {
+ result = cvm::to_str(obj->is_enabled(fid) ? 1 : 0);
+ return COLVARS_OK;
+ }
+
+ if (subcmd == "set") {
+ if (argc == 5) {
+ std::string const yesno =
+ colvarparse::to_lower_cppstr(std::string(argv[4]));
+ if ((yesno == std::string("yes")) ||
+ (yesno == std::string("on")) ||
+ (yesno == std::string("1"))) {
+ obj->enable(fid);
+ return COLVARS_OK;
+ } else if ((yesno == std::string("no")) ||
+ (yesno == std::string("off")) ||
+ (yesno == std::string("0"))) {
+ // TODO disable() function does not exist yet,
+ // dependencies will not be resolved
+ // obj->disable(fid);
+ obj->set_enabled(fid, false);
+ return COLVARS_OK;
+ }
+ }
+ result = "Syntax error\n" + help_string();
+ return COLVARSCRIPT_ERROR;
+ }
+ }
+ }
+
+ result = "Syntax error\n" + help_string();
+ return COLVARSCRIPT_ERROR;
+}
+
+
std::string colvarscript::help_string()
{
std::string buf;
@@ -459,6 +546,7 @@ Input and output:\n\
load <file name> -- load a state file (requires configuration)\n\
save <file name> -- save a state file (requires configuration)\n\
update -- recalculate colvars and biases\n\
+ addenergy <E> -- add <E> to the total bias energy\n\
printframe -- return a summary of the current frame\n\
printframelabels -- return labels to annotate printframe's output\n";
@@ -478,12 +566,17 @@ Accessing collective variables:\n\
colvar <name> addforce <F> -- apply given force on colvar <name>\n\
colvar <name> getconfig -- return config string of colvar <name>\n\
colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\
+ colvar <name> get <f> -- get the value of the colvar feature <f>\n\
+ colvar <name> set <f> <val> -- set the value of the colvar feature <f>\n\
\n\
Accessing biases:\n\
bias <name> energy -- return the current energy of bias <name>\n\
bias <name> update -- recalculate bias <name>\n\
bias <name> delete -- delete bias <name>\n\
- bias <name> getconfig -- return config string of bias <name>\n";
+ bias <name> getconfig -- return config string of bias <name>\n\
+ bias <name> get <f> -- get the value of the bias feature <f>\n\
+ bias <name> set <f> <val> -- set the value of the bias feature <f>\n\
+";
return buf;
}
diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h
index bfb8381e3b6c72bec02d3b6fe41e5037ab969bed..46b1ddd203afdcf72beb5da2e1ff2ffc2f01c65b 100644
--- a/lib/colvars/colvarscript.h
+++ b/lib/colvars/colvarscript.h
@@ -51,6 +51,10 @@ private:
/// Run subcommands on bias
int proc_bias(int argc, char const *argv[]);
+ /// Run subcommands on base colvardeps object (colvar, bias, ...)
+ int proc_features(colvardeps *obj,
+ int argc, char const *argv[]);
+
/// Builds and return a short help
std::string help_string(void);
};
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index 7320263ba7b27f3d6a79caf40002934405df1e91..89453ded9be77b795ea1b005e2262927442a1088 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -352,12 +352,17 @@ int colvarproxy_lammps::smp_enabled()
int colvarproxy_lammps::smp_colvars_loop()
{
colvarmodule *cv = this->colvars;
+ colvarproxy_lammps *proxy = (colvarproxy_lammps *) cv->proxy;
#pragma omp parallel for
- for (size_t i = 0; i < cv->colvars_smp.size(); i++) {
+ for (size_t i = 0; i < cv->variables_active_smp()->size(); i++) {
+ colvar *x = (*(cv->variables_active_smp()))[i];
+ int x_item = (*(cv->variables_active_smp_items()))[i];
if (cvm::debug()) {
- cvm::log("Calculating colvar \""+cv->colvars_smp[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n");
+ cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
+ "]: calc_colvars_items_smp(), i = "+cvm::to_str(i)+", cv = "+
+ x->name+", cvc = "+cvm::to_str(x_item)+"\n");
}
- cv->colvars_smp[i]->calc_cvcs(cv->colvars_smp_items[i], 1);
+ x->calc_cvcs(x_item, 1);
}
return cvm::get_error();
}
@@ -367,11 +372,13 @@ int colvarproxy_lammps::smp_biases_loop()
{
colvarmodule *cv = this->colvars;
#pragma omp parallel for
- for (size_t i = 0; i < cv->biases.size(); i++) {
+ for (size_t i = 0; i < cv->biases_active()->size(); i++) {
+ colvarbias *b = (*(cv->biases_active()))[i];
if (cvm::debug()) {
- cvm::log("Calculating bias \""+cv->biases[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n");
+ cvm::log("Calculating bias \""+b->name+"\" on thread "+
+ cvm::to_str(smp_thread_id())+"\n");
}
- cv->biases[i]->update();
+ b->update();
}
return cvm::get_error();
}
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index c0b4ea20fecb37dfad428b0e7058b4fa885652bc..ad63eb2f9eb6b740131d38aaac641b6f25df23e5 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -29,7 +29,7 @@
#endif
#ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2016-12-27"
+#define COLVARPROXY_VERSION "2017-01-09"
#endif
/* struct for packed data communication of coordinates and forces. */