jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update Colvars module to version 2017-03-10

This commit is contained in:
Giacomo Fiorin 2017-03-10 09:16:58 -05:00
parent 79b005dc3d
commit 5fe6206638
29 changed files with 1371 additions and 868 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

4
doc/src/fix_colvars.txt
View File

@ -59,7 +59,7 @@ always apply to the entire system and there can only be one instance
of the colvars fix at a time. The colvars fix will only communicate
the minimum information necessary and the colvars library supports
multiple, completely independent collective variables, so there is
no restriction to functionality by limiting the number of colvars fixes.
no restriction to functionaliry by limiting the number of colvars fixes.
The {input} keyword allows to specify a state file that would contain
the restart information required in order to continue a calculation from
@ -100,7 +100,7 @@ output"_thermo_style.html.
This fix computes a global scalar which can be accessed by various
"output commands"_Section_howto.html#howto_15. The scalar is the
cumulative energy change due to this fix. The scalar value
cummulative energy change due to this fix. The scalar value
calculated by this fix is "extensive".
[Restrictions:]

362
lib/colvars/colvar.cpp
View File

@ -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;
get_keyval(conf, "name", this->name,
(std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));
colvarmodule *cv = cvm::main();
if (cvm::colvar_by_name(this->name) != NULL) {
get_keyval(conf, "name", this->name,
(std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));
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);
cvm::log("Enabling the extended Lagrangian term for colvar \""+
this->name+"\".\n");
// The extended mass will apply forces
enable(f_cv_gradient);
enable(f_cv_extended_Lagrangian);
xr.type(value());
vr.type(value());
fr.type(value());
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);
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, "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");
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);
}
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());
}
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;

14
lib/colvars/colvar.h
View File

@ -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,

14
lib/colvars/colvaratoms.cpp
View File

@ -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);

51
lib/colvars/colvarbias.cpp
View File

@ -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]);
}
}

39
lib/colvars/colvarbias.h
View File

@ -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)

14
lib/colvars/colvarbias_abf.cpp
View File

@ -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) {

5
lib/colvars/colvarbias_alb.cpp
View File

@ -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;
}

10
lib/colvars/colvarbias_histogram.cpp
View File

@ -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 + "\"");
}
}

473
lib/colvars/colvarbias_meta.cpp
View File

@ -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();
std::vector<int> const curr_bin = values.size() ?
hills_energy->get_colvars_index(values) :
hills_energy->get_colvars_index();
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");
}
}
if (use_grids) {
// 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();
}
}
}
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 (hills_energy->index_ok(curr_bin)) {
// within the grid: add the energy and the forces from there
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 {
// off the grid: compute analytically only the hills at the grid's edges
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 (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);
}
}
} 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);
}
}
// 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())
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);
}
if (cvm::debug())
cvm::log("Hills energy = "+cvm::to_str(bias_energy)+
", hills forces = "+cvm::to_str(colvar_forces)+".\n");
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");
}
}
return COLVARS_OK;
}
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();
}
}
std::vector<int> const curr_bin = values.size() ?
hills_energy->get_colvars_index(values) :
hills_energy->get_colvars_index();
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);
}
}
}
// now include the hills that have not been binned yet (starting
// from new_hills_begin)
if (cvm::debug()) {
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding the forces from the other replicas.\n");
}
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");
}
}
}
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 &center = 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");

26
lib/colvars/colvarbias_meta.h
View File

@ -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

214
lib/colvars/colvarbias_restraint.cpp
View File

@ -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++) {

6
lib/colvars/colvarbias_restraint.h
View File

@ -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;

20
lib/colvars/colvarcomp.cpp
View File

@ -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;

9
lib/colvars/colvarcomp.h
View File

@ -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;
};

8
lib/colvars/colvarcomp_angles.cpp
View File

@ -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;

41
lib/colvars/colvarcomp_distances.cpp
View File

@ -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)

272
lib/colvars/colvardeps.cpp
View File

@ -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);
}
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);
init_feature(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);
init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
f_req_children(f_cvb_scalar_variables, f_cv_scalar);
init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
}
f_description(f_cvb_active, "active");
f_req_children(f_cvb_active, f_cv_active);
f_description(f_cvb_apply_force, "apply force");
f_req_children(f_cvb_apply_force, f_cv_gradient);
f_description(f_cvb_get_total_force, "obtain total force");
f_req_children(f_cvb_get_total_force, f_cv_total_force);
f_description(f_cvb_history_dependent, "history-dependent");
f_description(f_cvb_scalar_variables, "require scalar variables");
f_req_children(f_cvb_scalar_variables, f_cv_scalar);
// 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);
init_feature(f_cv_runave, "running average", f_type_user);
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_corrfunc, "correlation function", f_type_user);
f_description(f_cv_runave, "running average");
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++) {

86
lib/colvars/colvardeps.h
View File

@ -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

302
lib/colvars/colvarmodule.cpp
View File

@ -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;

115
lib/colvars/colvarmodule.h
View File

@ -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;
/// Prefix for all output files for this run
static std::string output_prefix;
private:
/// Prefix for all output files for this run
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
static std::vector<colvar *> colvars;
std::vector<colvar *> colvars;
/// Array of collective variables
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();
};

7
lib/colvars/colvarparse.cpp
View File

@ -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);

3
lib/colvars/colvarparse.h
View File

@ -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,

2
lib/colvars/colvarproxy.h
View File

@ -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;
}

117
lib/colvars/colvarscript.cpp
View File

@ -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;
}

4
lib/colvars/colvarscript.h
View File

@ -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);
};

19
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -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();
}

2
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -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. */