diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 5acaf8260c66e8d3269b6742fac7760b563eec63..5b776defca36759dadaaabebcc6df409000d52cd 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/doc/src/fix_colvars.txt b/doc/src/fix_colvars.txt
index 1991496d19c457323ab798185a251fa8e90d603d..d16853861f3c389a294b709ecbcc5cc7073a4a00 100644
--- a/doc/src/fix_colvars.txt
+++ b/doc/src/fix_colvars.txt
@@ -10,7 +10,7 @@ fix colvars command :h3
[Syntax:]
-fix ID group-ID colvars configfile keyword values ... :pre
+fix ID group-ID Colvars configfile keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
colvars = style name of this fix command :l
@@ -31,21 +31,19 @@ fix abf all colvars colvars.inp tstat 1 :pre
[Description:]
-This fix interfaces LAMMPS to a "collective variables" or "colvars"
-module library which allows to calculate potentials of mean force
+This fix interfaces LAMMPS to the collective variables "Colvars"
+library, which allows to calculate potentials of mean force
(PMFs) for any set of colvars, using different sampling methods:
currently implemented are the Adaptive Biasing Force (ABF) method,
metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
-(US) via a flexible harmonic restraint bias. The colvars library is
-hosted at "http://colvars.github.io/"_http://colvars.github.io/
+(US) via a flexible harmonic restraint bias.
This documentation describes only the fix colvars command itself and
LAMMPS specific parts of the code. The full documentation of the
colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf
-A detailed discussion of the implementation of the portable collective
-variable library is in "(Fiorin)"_#Fiorin. Additional information can
-be found in "(Henin)"_#Henin.
+The Colvars library is developed at "https://github.com/colvars/colvars"_https://github.com/colvars/colvars
+A detailed discussion of its implementation is in "(Fiorin)"_#Fiorin.
There are some example scripts for using this package with LAMMPS in the
examples/USER/colvars directory.
@@ -129,8 +127,3 @@ and tstat = NULL.
:link(Fiorin)
[(Fiorin)] Fiorin , Klein, Henin, Mol. Phys., DOI:10.1080/00268976.2013.813594
-
-:link(Henin)
-[(Henin)] Henin, Fiorin, Chipot, Klein, J. Chem. Theory Comput., 6,
-35-47 (2010)
-
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index b9c42d924752ff2a43ee370430931a5c88aafc30..95a40084f5854297cf2ac977d358a1a239536e5c 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// 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"
@@ -163,38 +170,38 @@ colvar::colvar(std::string const &conf)
}
feature_states[f_cv_homogeneous]->enabled = homogeneous;
}
- // Colvar is deemed periodic iff:
+
+ // Colvar is deemed periodic if:
// - it is homogeneous
// - all cvcs are periodic
// - all cvcs have the same period
-
- b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous);
- period = cvcs[0]->period;
- for (i = 1; i < cvcs.size(); i++) {
- if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
- b_periodic = false;
- period = 0.0;
+ if (cvcs[0]->b_periodic) { // TODO make this a CVC feature
+ bool b_periodic = true;
+ period = cvcs[0]->period;
+ for (i = 1; i < cvcs.size(); i++) {
+ if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
+ b_periodic = false;
+ period = 0.0;
+ cvm::log("Warning: although one component is periodic, this colvar will "
+ "not be treated as periodic, either because the exponent is not "
+ "1, or because components of different periodicity are defined. "
+ "Make sure that you know what you are doing!");
+ }
}
+ feature_states[f_cv_periodic]->enabled = b_periodic;
}
- feature_states[f_cv_periodic]->enabled = b_periodic;
// check that cvcs are compatible
for (i = 0; i < cvcs.size(); i++) {
- if ((cvcs[i])->b_periodic && !b_periodic) {
- cvm::log("Warning: although this component is periodic, the colvar will "
- "not be treated as periodic, either because the exponent is not "
- "1, or because multiple components are present. Make sure that "
- "you know what you are doing!");
- }
// components may have different types only for scripted functions
if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),
- cvcs[0]->value())) ) {
+ cvcs[0]->value())) ) {
cvm::error("ERROR: you are definining this collective variable "
"by using components of different types. "
"You must use the same type in order to "
- " sum them together.\n", INPUT_ERROR);
+ "sum them together.\n", INPUT_ERROR);
return;
}
}
@@ -207,16 +214,15 @@ colvar::colvar(std::string const &conf)
// at this point, the colvar's type is defined
f.type(value());
f_accumulated.type(value());
- fb.type(value());
+
+ reset_bias_force();
get_keyval(conf, "width", width, 1.0);
if (width <= 0.0) {
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
+ return;
}
- // NOTE: not porting wall stuff to new deps, as this will change to a separate bias
- // the grid functions will wait a little as well
-
lower_boundary.type(value());
lower_wall.type(value());
@@ -400,6 +406,9 @@ colvar::colvar(std::string const &conf)
f_old.type(value());
f_old.reset();
+ x_restart.type(value());
+ after_restart = false;
+
if (cvm::b_analysis)
parse_analysis(conf);
@@ -761,9 +770,13 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
return error_code;
}
+ if (cvm::step_relative() > 0) {
+ // Total force depends on Jacobian derivative from previous timestep
+ error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
+ }
+ // atom coordinates are updated by the next line
error_code |= calc_cvc_values(first_cvc, num_cvcs);
error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
- error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
if (cvm::debug())
@@ -780,9 +793,12 @@ int colvar::collect_cvc_data()
int error_code = COLVARS_OK;
+ if (cvm::step_relative() > 0) {
+ // Total force depends on Jacobian derivative from previous timestep
+ error_code |= collect_cvc_total_forces();
+ }
error_code |= collect_cvc_values();
error_code |= collect_cvc_gradients();
- error_code |= collect_cvc_total_forces();
error_code |= collect_cvc_Jacobians();
error_code |= calc_colvar_properties();
@@ -872,6 +888,22 @@ int colvar::collect_cvc_values()
cvm::log("Colvar \""+this->name+"\" has value "+
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) {
+ cvm::error("Error: the calculated value of colvar \""+name+
+ "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
+ "last read from the state file:\n"+cvm::to_str(x_restart)+
+ "\nPossible causes are changes in configuration, "
+ "wrong state file, or how PBC wrapping is handled.\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ }
+ }
+
return COLVARS_OK;
}
@@ -979,22 +1011,17 @@ int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
if (cvm::debug())
cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
- // if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
- // Disabled check to allow for explicit total force calculation
- // even with extended Lagrangian
+ cvm::increase_depth();
- if (cvm::step_relative() > 0) {
- cvm::increase_depth();
- // get from the cvcs the total forces from the PREVIOUS step
- for (i = first_cvc, cvc_count = 0;
- (i < cvcs.size()) && (cvc_count < cvc_max_count);
- i++) {
- if (!cvcs[i]->is_enabled()) continue;
- cvc_count++;
- (cvcs[i])->calc_force_invgrads();
- }
- cvm::decrease_depth();
+ for (i = first_cvc, cvc_count = 0;
+ (i < cvcs.size()) && (cvc_count < cvc_max_count);
+ i++) {
+ if (!cvcs[i]->is_enabled()) continue;
+ cvc_count++;
+ (cvcs[i])->calc_force_invgrads();
}
+ cvm::decrease_depth();
+
if (cvm::debug())
cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
@@ -1013,6 +1040,11 @@ int colvar::collect_cvc_total_forces()
// get from the cvcs the total forces from the PREVIOUS step
for (size_t i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
+ if (cvm::debug())
+ cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+ " within colvar \""+this->name+"\" has total force "+
+ cvm::to_str((cvcs[i])->total_force(),
+ cvm::cv_width, cvm::cv_prec)+".\n");
// linear combination is assumed
ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
}
@@ -1056,6 +1088,11 @@ int colvar::collect_cvc_Jacobians()
fj.reset();
for (size_t i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
+ if (cvm::debug())
+ cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+ " within colvar \""+this->name+"\" has Jacobian derivative"+
+ cvm::to_str((cvcs[i])->Jacobian_derivative(),
+ cvm::cv_width, cvm::cv_prec)+".\n");
// linear combination is assumed
fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
}
@@ -1128,15 +1165,16 @@ cvm::real colvar::update_forces_energy()
if (is_enabled(f_cv_Jacobian)) {
// the instantaneous Jacobian force was not included in the reported total force;
// instead, it is subtracted from the applied force (silent Jacobian correction)
+ // This requires the Jacobian term for the *current* timestep
if (is_enabled(f_cv_hide_Jacobian))
f -= fj;
}
- if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
+ // Wall force
+ colvarvalue fw(x);
+ fw.reset();
- // 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");
@@ -1144,12 +1182,11 @@ cvm::real colvar::update_forces_energy()
// 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_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {
+ (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
- cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
+ cvm::real const grad = this->dist2_lgrad(x, lower_wall);
if (grad < 0.0) {
fw = -0.5 * lower_wall_k * grad;
- f += fw;
if (cvm::debug())
cvm::log("Applying a lower wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
@@ -1157,10 +1194,9 @@ cvm::real colvar::update_forces_energy()
} else {
- cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
+ cvm::real const grad = this->dist2_lgrad(x, upper_wall);
if (grad > 0.0) {
fw = -0.5 * upper_wall_k * grad;
- f += fw;
if (cvm::debug())
cvm::log("Applying an upper wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
@@ -1168,6 +1204,9 @@ cvm::real colvar::update_forces_energy()
}
}
+ // 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)) {
cvm::real dt = cvm::dt();
@@ -1175,11 +1214,12 @@ cvm::real colvar::update_forces_energy()
f_ext.reset();
// the total force is applied to the fictitious mass, while the
- // atoms only feel the harmonic force
+ // atoms only feel the harmonic force + wall force
// fr: bias force on extended variable (without harmonic spring), for output in trajectory
// f_ext: total force on extended variable (including harmonic spring)
- // f: - initially, external biasing force (including wall forces)
- // - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
+ // f: - initially, external biasing force
+ // - after this code block, colvar force to be applied to atomic coordinates
+ // ie. spring force + wall force
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);
@@ -1209,9 +1249,16 @@ cvm::real colvar::update_forces_energy()
vr += (0.5 * dt) * f_ext / ext_mass;
xr += dt * vr;
xr.apply_constraints();
- if (this->b_periodic) this->wrap(xr);
+ 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;
+
+ // Store force to be applied, possibly summed over several timesteps
f_accumulated += f;
if (is_enabled(f_cv_fdiff_velocity)) {
@@ -1428,14 +1475,15 @@ std::istream & colvar::read_restart(std::istream &is)
}
}
- if ( !(get_keyval(conf, "x", x,
- colvarvalue(x.type()), colvarparse::parse_silent)) ) {
+ if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"the value of the colvar \""+
name+"\" .\n");
} else {
cvm::log("Restarting collective variable \""+name+"\" from value: "+
cvm::to_str(x)+"\n");
+ x_restart = x;
+ after_restart = true;
}
if (is_enabled(f_cv_extended_Lagrangian)) {
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index d7b1e6076f981600303aeba1238f2b3731212e44..2cf3d2dac512b7492452940393067e8d6f7088ef 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVAR_H
#define COLVAR_H
@@ -170,6 +177,9 @@ public:
/// the biases are updated
colvarvalue fb;
+ /// \brief Bias force to the actual value (only useful with extended Lagrangian)
+ colvarvalue fb_actual;
+
/// \brief Total \em applied force; fr (if extended_lagrangian
/// is defined), fb (if biases are applied) and the walls' forces
/// (if defined) contribute to it
@@ -183,13 +193,9 @@ public:
colvarvalue ft;
- /// Period, if it is a constant
+ /// Period, if this variable is periodic
cvm::real period;
- /// \brief Same as above, but also takes into account components
- /// with a variable period, such as distanceZ
- bool b_periodic;
-
/// \brief Expand the boundaries of multiples of width, to keep the
/// value always within range
@@ -290,6 +296,9 @@ public:
/// Add to the total force from biases
void add_bias_force(colvarvalue const &force);
+ /// Apply a force to the actual value (only meaningful with extended Lagrangian)
+ void add_bias_force_actual_value(colvarvalue const &force);
+
/// \brief Collect all forces on this colvar, integrate internal
/// equations of motion of internal degrees of freedom; see also
/// colvar::communicate_forces()
@@ -386,6 +395,12 @@ protected:
/// Previous value (to calculate velocities during analysis)
colvarvalue x_old;
+ /// Value read from the most recent state file (if any)
+ colvarvalue x_restart;
+
+ /// True if a state file was just read
+ bool after_restart;
+
/// Time series of values and velocities used in correlation
/// functions
std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
@@ -577,9 +592,20 @@ inline void colvar::add_bias_force(colvarvalue const &force)
}
+inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
+{
+ if (cvm::debug()) {
+ cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
+ }
+ fb_actual += force;
+}
+
+
inline void colvar::reset_bias_force() {
fb.type(value());
fb.reset();
+ fb_actual.type(value());
+ fb_actual.reset();
}
#endif
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index e8046a97d5d2db92b3048d6106b01da8d2ff0452..48c16e887af95c9e623aa3b8b525b54594e10ca2 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvaratoms.h"
@@ -171,7 +178,10 @@ int cvm::atom_group::remove_atom(cvm::atom_iter ai)
int cvm::atom_group::init()
{
- if (!key.size()) key = "atoms";
+ if (!key.size()) key = "unnamed";
+ description = "atom group " + key;
+ // These will be overwritten by parse(), if initializing from a config string
+
atoms.clear();
// TODO: check with proxy whether atom forces etc are available
@@ -179,6 +189,7 @@ int cvm::atom_group::init()
index = -1;
+ b_dummy = false;
b_center = false;
b_rotate = false;
b_user_defined_fit = false;
@@ -440,6 +451,7 @@ int cvm::atom_group::parse(std::string const &conf)
if (b_print_atom_ids) {
cvm::log("Internal definition of the atom group:\n");
+ cvm::log(print_atom_ids());
}
cvm::decrease_depth();
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index a2a771a8a3f3f8a181c2f05035a826cc68b6a2a9..85f6212951bec0160e504cde644add9787c202a8 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARATOMS_H
#define COLVARATOMS_H
@@ -64,7 +71,7 @@ public:
/// \brief Initialize an atom for collective variable calculation
/// and get its internal identifier \param atom_number Atom index in
- /// the system topology (starting from 1)
+ /// the system topology (1-based)
atom(int atom_number);
/// \brief Initialize an atom for collective variable calculation
@@ -453,6 +460,8 @@ public:
/// are not used, either because they were not defined (e.g because
/// the colvar has not a scalar value) or the biases require to
/// micromanage the force.
+ /// This function will be phased out eventually, in favor of
+ /// apply_colvar_force() once that is implemented for non-scalar values
void apply_force(cvm::rvector const &force);
};
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 76d8489de995f8ca50552ae388d65c563349b0c5..fdd2b6254ca12dd8b450244f72c75f8b570f8a2f 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias.h"
@@ -29,6 +36,7 @@ colvarbias::colvarbias(char const *key)
has_data = false;
b_output_energy = false;
reset();
+ state_file_step = 0;
// Start in active state by default
enable(f_cvb_active);
@@ -141,20 +149,25 @@ int colvarbias::clear()
int colvarbias::add_colvar(std::string const &cv_name)
{
if (colvar *cv = cvm::colvar_by_name(cv_name)) {
- // Removed this as nor all biases apply forces eg histogram
- // cv->enable(colvar::task_gradients);
+
if (cvm::debug()) {
cvm::log("Applying this bias to collective variable \""+
cv->name+"\".\n");
}
+
colvars.push_back(cv);
colvar_forces.push_back(colvarvalue());
colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
+ colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
colvar_forces.back().reset();
cv->biases.push_back(this); // add back-reference to this bias to colvar
+ if (is_enabled(f_cvb_apply_force)) {
+ cv->enable(f_cv_gradient);
+ }
+
// Add dependency link.
// All biases need at least the value of each colvar
// although possibly not at all timesteps
@@ -165,6 +178,7 @@ int colvarbias::add_colvar(std::string const &cv_name)
cv_name+"\".\n", INPUT_ERROR);
return INPUT_ERROR;
}
+
return COLVARS_OK;
}
@@ -190,16 +204,19 @@ void colvarbias::communicate_forces()
}
-void colvarbias::change_configuration(std::string const &conf)
+int colvarbias::change_configuration(std::string const &conf)
{
- cvm::error("Error: change_configuration() not implemented.\n");
+ cvm::error("Error: change_configuration() not implemented.\n",
+ COLVARS_NOT_IMPLEMENTED);
+ return COLVARS_NOT_IMPLEMENTED;
}
cvm::real colvarbias::energy_difference(std::string const &conf)
{
- cvm::error("Error: energy_difference() not implemented.\n");
- return 0.;
+ cvm::error("Error: energy_difference() not implemented.\n",
+ COLVARS_NOT_IMPLEMENTED);
+ return 0.0;
}
@@ -225,6 +242,118 @@ int colvarbias::replica_share()
return COLVARS_NOT_IMPLEMENTED;
}
+
+std::string const colvarbias::get_state_params() const
+{
+ std::ostringstream os;
+ os << "step " << cvm::step_absolute() << "\n"
+ << "name " << this->name << "\n";
+ return os.str();
+}
+
+
+int colvarbias::set_state_params(std::string const &conf)
+{
+ std::string new_name = "";
+ if (colvarparse::get_keyval(conf, "name", new_name,
+ std::string(""), colvarparse::parse_silent) &&
+ (new_name != this->name)) {
+ cvm::error("Error: in the state file, the "
+ "\""+bias_type+"\" block has a different name, \""+new_name+
+ "\": different system?\n", INPUT_ERROR);
+ }
+
+ if (name.size() == 0) {
+ cvm::error("Error: \""+bias_type+"\" block within the restart file "
+ "has no identifiers.\n", INPUT_ERROR);
+ }
+
+ colvarparse::get_keyval(conf, "step", state_file_step,
+ cvm::step_absolute(), colvarparse::parse_silent);
+
+ return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias::write_state(std::ostream &os)
+{
+ if (cvm::debug()) {
+ cvm::log("Writing state file for bias \""+name+"\"\n");
+ }
+ os.setf(std::ios::scientific, std::ios::floatfield);
+ os.precision(cvm::cv_prec);
+ os << bias_type << " {\n"
+ << " configuration {\n";
+ std::istringstream is(get_state_params());
+ std::string line;
+ while (std::getline(is, line)) {
+ os << " " << line << "\n";
+ }
+ os << " }\n";
+ write_state_data(os);
+ os << "}\n\n";
+ return os;
+}
+
+
+std::istream & colvarbias::read_state(std::istream &is)
+{
+ size_t const start_pos = is.tellg();
+
+ std::string key, brace, conf;
+ if ( !(is >> key) || !(key == bias_type) ||
+ !(is >> brace) || !(brace == "{") ||
+ !(is >> colvarparse::read_block("configuration", conf)) ||
+ (set_state_params(conf) != COLVARS_OK) ) {
+ cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
+ this->name+"\" at position "+
+ cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+ is.clear();
+ is.seekg(start_pos, std::ios::beg);
+ is.setstate(std::ios::failbit);
+ return is;
+ }
+
+ if (!read_state_data(is)) {
+ cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
+ this->name+"\" at position "+
+ cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+ is.clear();
+ is.seekg(start_pos, std::ios::beg);
+ is.setstate(std::ios::failbit);
+ }
+
+ is >> brace;
+ if (brace != "}") {
+ cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
+ this->name+"\": no matching brace at position "+
+ cvm::to_str(is.tellg())+" in stream.\n");
+ is.setstate(std::ios::failbit);
+ }
+
+ return is;
+}
+
+
+std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
+{
+ size_t const start_pos = is.tellg();
+ std::string key_in;
+ if ( !(is >> key_in) ||
+ !(key_in == to_lower_cppstr(std::string(key))) ) {
+ cvm::error("Error: in reading restart configuration for "+
+ bias_type+" bias \""+this->name+"\" at position "+
+ cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+ is.clear();
+ is.seekg(start_pos, std::ios::beg);
+ is.setstate(std::ios::failbit);
+ return is;
+ }
+ return is;
+}
+
+
+
std::ostream & colvarbias::write_traj_label(std::ostream &os)
{
os << " ";
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index 8238f58f402bf1aeedb65a9b7adebe3ae7a4cfee..12397dcb8fa0ba8aa2356fd776d1a8ccb16b4546 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_H
#define COLVARBIAS_H
@@ -9,7 +16,8 @@
/// \brief Collective variable bias, base class
-class colvarbias : public colvarparse, public colvardeps {
+class colvarbias
+ : public virtual colvarparse, public virtual colvardeps {
public:
/// Name of this bias
@@ -24,6 +32,12 @@ 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
+ {
+ return colvars.size();
+ }
+
/// Retrieve colvar values and calculate their biasing forces
/// Return bias energy
virtual int update();
@@ -31,7 +45,7 @@ public:
// TODO: move update_bias here (share with metadynamics)
/// Load new configuration - force constant and/or centers only
- virtual void change_configuration(std::string const &conf);
+ virtual int change_configuration(std::string const &conf);
/// Calculate change in energy from using alternate configuration
virtual cvm::real energy_difference(std::string const &conf);
@@ -49,7 +63,7 @@ public:
virtual void analyze() {}
/// Send forces to the collective variables
- void communicate_forces();
+ virtual void communicate_forces();
/// \brief Constructor
colvarbias(char const *key);
@@ -60,13 +74,11 @@ public:
/// \brief Set to zero all mutable data
virtual int reset();
-protected:
+private:
/// Default constructor
colvarbias();
-private:
-
/// Copy constructor
colvarbias(colvarbias &);
@@ -78,28 +90,59 @@ public:
/// Destructor
virtual ~colvarbias();
- /// Read the bias configuration from a restart file
- virtual std::istream & read_restart(std::istream &is) = 0;
+ /// Write the values of specific mutable properties to a string
+ virtual std::string const get_state_params() const;
+
+ /// Read the values of specific mutable properties from a string
+ virtual int set_state_params(std::string const &state_conf);
+
+ /// Write all mutable data not already written by get_state_params()
+ virtual std::ostream & write_state_data(std::ostream &os)
+ {
+ return os;
+ }
+
+ /// Read all mutable data not already set by set_state_params()
+ virtual std::istream & read_state_data(std::istream &is)
+ {
+ return is;
+ }
+
+ /// Read a keyword from the state data (typically a header)
+ std::istream & read_state_data_key(std::istream &is, char const *key);
- /// Write the bias configuration to a restart file
- virtual std::ostream & write_restart(std::ostream &os) = 0;
+ /// Write the bias configuration to a restart file or other stream
+ virtual std::ostream & write_state(std::ostream &os);
+
+ /// Read the bias configuration from a restart file or other stream
+ virtual std::istream & read_state(std::istream &is);
/// Write a label to the trajectory file (comment line)
virtual std::ostream & write_traj_label(std::ostream &os);
- /// (Re)initialize the output files (does not write them yet)
- virtual int setup_output() { return COLVARS_OK; }
-
/// Output quantities such as the bias energy to the trajectory file
virtual std::ostream & write_traj(std::ostream &os);
- /// Write output files (if defined, e.g. in analysis mode)
+ /// (Re)initialize the output files (does not write them yet)
+ virtual int setup_output()
+ {
+ return COLVARS_OK;
+ }
+
+ /// Write any output files that this bias may have (e.g. PMF files)
virtual int write_output_files()
{
return COLVARS_OK;
}
- inline cvm::real get_energy() {
+ /// If this bias is communicating with other replicas through files, send it to them
+ virtual int write_state_to_replicas()
+ {
+ return COLVARS_OK;
+ }
+
+ inline cvm::real get_energy()
+ {
return bias_energy;
}
@@ -107,9 +150,11 @@ public:
static std::vector<feature *> cvb_features;
/// \brief Implementation of the feature list accessor for colvarbias
- virtual std::vector<feature *> &features() {
+ virtual std::vector<feature *> &features()
+ {
return cvb_features;
}
+
protected:
/// \brief Pointers to collective variables to which the bias is
@@ -130,6 +175,9 @@ protected:
/// (for history-dependent biases)
bool has_data;
+ /// \brief Step number read from the last state file
+ size_t state_file_step;
+
};
#endif
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index 1665ea61a9c577afea8cbe3a3d8c680729581aa2..ccfa401697693bb97a2d403c6cc06da3f2f32d9e 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvar.h"
#include "colvarbias_abf.h"
@@ -23,6 +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);
// TODO relax this in case of VMD plugin
@@ -590,16 +599,10 @@ void colvarbias_abf::read_gradients_samples()
}
-std::ostream & colvarbias_abf::write_restart(std::ostream& os)
+std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
{
-
std::ios::fmtflags flags(os.flags());
- os << "abf {\n"
- << " configuration {\n"
- << " name " << this->name << "\n";
- os << " }\n";
-
os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
os << "\nsamples\n";
samples->write_raw(os, 8);
@@ -617,117 +620,47 @@ std::ostream & colvarbias_abf::write_restart(std::ostream& os)
z_gradients->write_raw(os, 8);
}
- os << "}\n\n";
-
os.flags(flags);
return os;
}
-std::istream & colvarbias_abf::read_restart(std::istream& is)
+std::istream & colvarbias_abf::read_state_data(std::istream& is)
{
if ( input_prefix.size() > 0 ) {
cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR);
}
- size_t const start_pos = is.tellg();
-
- cvm::log("Restarting ABF bias \""+
- this->name+"\".\n");
- std::string key, brace, conf;
-
- if ( !(is >> key) || !(key == "abf") ||
- !(is >> brace) || !(brace == "{") ||
- !(is >> colvarparse::read_block("configuration", conf)) ) {
- cvm::log("Error: in reading restart configuration for ABF bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
- return is;
- }
-
- std::string name = "";
- if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
- (name != this->name) )
- cvm::error("Error: in the restart file, the "
- "\"abf\" block has wrong name(" + name + ")\n");
- if ( name == "" ) {
- cvm::error("Error: \"abf\" block in the restart file has no name.\n");
- }
-
- if ( !(is >> key) || !(key == "samples")) {
- cvm::log("Error: in reading restart configuration for ABF bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
+ if (! read_state_data_key(is, "samples")) {
return is;
}
if (! samples->read_raw(is)) {
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
return is;
}
- if ( !(is >> key) || !(key == "gradient")) {
- cvm::log("Error: in reading restart configuration for ABF bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
+ if (! read_state_data_key(is, "gradient")) {
return is;
}
if (! gradients->read_raw(is)) {
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
return is;
}
if (z_gradients) {
- if ( !(is >> key) || !(key == "z_samples")) {
- cvm::log("Error: in reading restart configuration for ABF bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
+
+ if (! read_state_data_key(is, "z_samples")) {
return is;
}
if (! z_samples->read_raw(is)) {
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
return is;
}
- if ( !(is >> key) || !(key == "z_gradient")) {
- cvm::log("Error: in reading restart configuration for ABF bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
+ if (! read_state_data_key(is, "z_gradient")) {
return is;
}
if (! z_gradients->read_raw(is)) {
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
return is;
}
}
- is >> brace;
- if (brace != "}") {
- cvm::error("Error: corrupt restart information for ABF bias \""+
- this->name+"\": no matching brace at position "+
- cvm::to_str(is.tellg())+" in the restart file.\n");
- is.setstate(std::ios::failbit);
- }
+
return is;
}
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index ee7624096835e52bf27dc02cbb1b96d10c6e9c81..41a5475fa70087ca709e4ca2be9f0f118874a311 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_ABF_H
#define COLVARBIAS_ABF_H
@@ -107,8 +114,8 @@ private:
/// Read human-readable FE gradients and sample count (if not using restart)
void read_gradients_samples();
- std::istream& read_restart(std::istream&);
- std::ostream& write_restart(std::ostream&);
+ std::istream& read_state_data(std::istream&);
+ std::ostream& write_state_data(std::ostream&);
};
#endif
diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp
index 8c8956082e61ebb1f938e9cb6e57c567e125d441..9e414d4db49637087dc954ff8239a5465287664e 100644
--- a/lib/colvars/colvarbias_alb.cpp
+++ b/lib/colvars/colvarbias_alb.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
@@ -33,6 +40,9 @@ 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;
@@ -239,37 +249,8 @@ int colvarbias_alb::update()
}
-std::istream & colvarbias_alb::read_restart(std::istream &is)
+int colvarbias_alb::set_state_params(std::string const &conf)
{
- size_t const start_pos = is.tellg();
-
- cvm::log("Restarting adaptive linear bias \""+
- this->name+"\".\n");
-
- std::string key, brace, conf;
- if ( !(is >> key) || !(key == "ALB") ||
- !(is >> brace) || !(brace == "{") ||
- !(is >> colvarparse::read_block("configuration", conf)) ) {
-
- cvm::log("Error: in reading restart configuration for restraint bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
- return is;
- }
-
- std::string name = "";
- if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
- (name != this->name) )
- cvm::fatal_error("Error: in the restart file, the "
- "\"ALB\" block has a wrong name\n");
- if (name.size() == 0) {
- cvm::fatal_error("Error: \"ALB\" block in the restart file "
- "has no identifiers.\n");
- }
-
if (!get_keyval(conf, "setCoupling", set_coupling))
cvm::fatal_error("Error: current setCoupling is missing from the restart.\n");
@@ -299,23 +280,13 @@ std::istream & colvarbias_alb::read_restart(std::istream &is)
if (!get_keyval(conf, "b_equilibration", b_equilibration))
cvm::fatal_error("Error: current updateCalls is missing from the restart.\n");
- is >> brace;
- if (brace != "}") {
- cvm::fatal_error("Error: corrupt restart information for adaptive linear bias \""+
- this->name+"\": no matching brace at position "+
- cvm::to_str(is.tellg())+" in the restart file.\n");
- is.setstate(std::ios::failbit);
- }
-
- return is;
+ return COLVARS_OK;
}
-std::ostream & colvarbias_alb::write_restart(std::ostream &os)
+std::string const colvarbias_alb::get_state_params() const
{
- os << "ALB {\n"
- << " configuration {\n"
- << " name " << this->name << "\n";
+ std::ostringstream os;
os << " setCoupling ";
size_t i;
for (i = 0; i < colvars.size(); i++) {
@@ -358,10 +329,7 @@ std::ostream & colvarbias_alb::write_restart(std::ostream &os)
else
os << " b_equilibration no\n";
- os << " }\n"
- << "}\n\n";
-
- return os;
+ return os.str();
}
diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h
index a3c6133db2503110c450c97523ad44b87f81199b..d003c244e96855486b37bf6bb2ec94a4b881a594 100644
--- a/lib/colvars/colvarbias_alb.h
+++ b/lib/colvars/colvarbias_alb.h
@@ -1,10 +1,18 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_ALB_H
#define COLVARBIAS_ALB_H
#include "colvar.h"
-#include "colvarbias_restraint.h"
+#include "colvarbias.h"
+
class colvarbias_alb : public colvarbias {
@@ -15,8 +23,8 @@ public:
virtual int init(std::string const &conf);
virtual int update();
- virtual std::istream & read_restart(std::istream &is);
- virtual std::ostream & write_restart(std::ostream &os);
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
virtual std::ostream & write_traj_label(std::ostream &os);
virtual std::ostream & write_traj(std::ostream &os);
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index eef4b506e30e87a00076ed38846d840ad11918d0..900ad213d51d6d0301063d6edba45561e285a0e7 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvar.h"
#include "colvarbias_histogram.h"
@@ -17,6 +24,9 @@ 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);
@@ -196,78 +206,25 @@ int colvarbias_histogram::write_output_files()
}
-std::istream & colvarbias_histogram::read_restart(std::istream& is)
+std::istream & colvarbias_histogram::read_state_data(std::istream& is)
{
- size_t const start_pos = is.tellg();
-
- cvm::log("Restarting collective variable histogram \""+
- this->name+"\".\n");
- std::string key, brace, conf;
-
- if ( !(is >> key) || !(key == "histogram") ||
- !(is >> brace) || !(brace == "{") ||
- !(is >> colvarparse::read_block("configuration", conf)) ) {
- cvm::log("Error: in reading restart configuration for histogram \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
- return is;
- }
-
- int id = -1;
- std::string name = "";
- if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
- (name != this->name) )
- cvm::error("Error: in the restart file, the "
- "\"histogram\" block has a wrong name: different system?\n");
- if ( (id == -1) && (name == "") ) {
- cvm::error("Error: \"histogram\" block in the restart file "
- "has no name.\n");
- }
-
- if ( !(is >> key) || !(key == "grid")) {
- cvm::error("Error: in reading restart configuration for histogram \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
+ if (! read_state_data_key(is, "grid")) {
return is;
}
if (! grid->read_raw(is)) {
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
return is;
}
- is >> brace;
- if (brace != "}") {
- cvm::error("Error: corrupt restart information for ABF bias \""+
- this->name+"\": no matching brace at position "+
- cvm::to_str(is.tellg())+" in the restart file.\n");
- is.setstate(std::ios::failbit);
- }
return is;
}
-std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
+
+std::ostream & colvarbias_histogram::write_state_data(std::ostream& os)
{
std::ios::fmtflags flags(os.flags());
os.setf(std::ios::fmtflags(0), std::ios::floatfield);
-
- os << "histogram {\n"
- << " configuration {\n"
- << " name " << this->name << "\n";
- os << " }\n";
-
os << "grid\n";
grid->write_raw(os, 8);
-
- os << "}\n\n";
-
os.flags(flags);
return os;
}
diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h
index cca49b9901834faf7e52a3cabca21ce8006a9611..b9c1b499502d9de0995396e3a8d5a83de2cf8ca7 100644
--- a/lib/colvars/colvarbias_histogram.h
+++ b/lib/colvars/colvarbias_histogram.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_HISTOGRAM_H
#define COLVARBIAS_HISTOGRAM_H
@@ -35,8 +42,8 @@ protected:
/// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram
std::vector<cvm::real> weights;
- virtual std::istream& read_restart(std::istream&);
- virtual std::ostream& write_restart(std::ostream&);
+ virtual std::istream & read_state_data(std::istream &is);
+ virtual std::ostream & write_state_data(std::ostream &os);
};
#endif
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index 285420932b385e8509f8ce00a5e35a2d8367bacc..01a04d1a016dbb116f707296cfd54bb99ad86f83 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <iostream>
#include <sstream>
#include <fstream>
@@ -26,10 +33,9 @@
colvarbias_meta::colvarbias_meta(char const *key)
- : colvarbias(key),
- new_hills_begin(hills.end()),
- state_file_step(0)
+ : colvarbias(key)
{
+ new_hills_begin = hills.end();
}
@@ -164,11 +170,32 @@ int colvarbias_meta::init(std::string const &conf)
target_dist = NULL;
get_keyval(conf, "ebMeta", ebmeta, false);
if(ebmeta){
+ if (use_grids && expand_grids) {
+ cvm::fatal_error("Error: expandBoundaries is not supported with "
+ "ebMeta please allocate wide enough boundaries for "
+ "each colvar ahead of time and set targetdistfile "
+ "accordingly. \n");
+ }
target_dist = new colvar_grid_scalar();
target_dist->init_from_colvars(colvars);
get_keyval(conf, "targetdistfile", target_dist_file);
std::ifstream targetdiststream(target_dist_file.c_str());
target_dist->read_multicol(targetdiststream);
+ cvm::real min_val = target_dist->minimum_value();
+ if(min_val<0){
+ cvm::error("Error: Target distribution of ebMeta "
+ "has negative values!.\n", INPUT_ERROR);
+ }
+ cvm::real min_pos_val = target_dist->minimum_pos_value();
+ if(min_pos_val<=0){
+ cvm::error("Error: Target distribution of ebMeta has negative "
+ "or zero minimum positive value!.\n", INPUT_ERROR);
+ }
+ if(min_val==0){
+ cvm::log("WARNING: Target distribution has zero values.\n");
+ cvm::log("Zeros will be converted to the minimum positive value.\n");
+ target_dist->remove_zeros(min_pos_val);
+ }
// normalize target distribution and multiply by effective volume = exp(differential entropy)
target_dist->multiply_constant(1.0/target_dist->integral());
cvm::real volume = std::exp(target_dist->entropy());
@@ -405,7 +432,9 @@ int colvarbias_meta::update()
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(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(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;
}
}
@@ -972,7 +1001,7 @@ void colvarbias_meta::read_replica_files()
(replicas[ir])->replica_state_file+"\".\n");
std::ifstream is((replicas[ir])->replica_state_file.c_str());
- if (! (replicas[ir])->read_restart(is)) {
+ if (! (replicas[ir])->read_state(is)) {
cvm::log("Reading from file \""+(replicas[ir])->replica_state_file+
"\" failed or incomplete: will try again in "+
cvm::to_str(replica_update_freq)+" steps.\n");
@@ -1068,63 +1097,24 @@ void colvarbias_meta::read_replica_files()
}
-// **********************************************************************
-// input functions
-// **********************************************************************
-
-
-std::istream & colvarbias_meta::read_restart(std::istream& is)
+int colvarbias_meta::set_state_params(std::string const &state_conf)
{
- size_t const start_pos = is.tellg();
-
- if (comm == single_replica) {
- // if using a multiple replicas scheme, output messages
- // are printed before and after calling this function
- cvm::log("Restarting metadynamics bias \""+this->name+"\""+
- ".\n");
- }
- std::string key, brace, conf;
- if ( !(is >> key) || !(key == "metadynamics") ||
- !(is >> brace) || !(brace == "{") ||
- !(is >> colvarparse::read_block("configuration", conf)) ) {
-
- if (comm == single_replica)
- cvm::log("Error: in reading restart configuration for metadynamics bias \""+
- this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- (replica_state_file_in_sync ? ("at position "+
- cvm::to_str(start_pos)+
- " in the state file") : "")+".\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
- return is;
+ std::string new_replica = "";
+ if (colvarparse::get_keyval(state_conf, "replicaID", new_replica,
+ std::string(""), colvarparse::parse_silent) &&
+ (new_replica != this->replica_id)) {
+ cvm::error("Error: in the state file, the "
+ "\"metadynamics\" block has a different replicaID: different system?\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
}
- std::string name = "";
- if ( colvarparse::get_keyval(conf, "name", name,
- std::string(""), colvarparse::parse_silent) &&
- (name != this->name) )
- cvm::fatal_error("Error: in the restart file, the "
- "\"metadynamics\" block has a different name: different system?\n");
-
- if (name.size() == 0) {
- cvm::fatal_error("Error: \"metadynamics\" block within the restart file "
- "has no identifiers.\n");
- }
-
- if (comm != single_replica) {
- std::string replica = "";
- if (colvarparse::get_keyval(conf, "replicaID", replica,
- std::string(""), colvarparse::parse_silent) &&
- (replica != this->replica_id))
- cvm::fatal_error("Error: in the restart file, the "
- "\"metadynamics\" block has a different replicaID: different system?\n");
+ return COLVARS_OK;
+}
- colvarparse::get_keyval(conf, "step", state_file_step,
- cvm::step_absolute(), colvarparse::parse_silent);
- }
+std::istream & colvarbias_meta::read_state_data(std::istream& is)
+{
bool grids_from_restart_file = use_grids;
if (use_grids) {
@@ -1155,6 +1145,7 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
}
size_t const hills_energy_pos = is.tellg();
+ std::string key;
if (!(is >> key)) {
if (hills_energy_backup != NULL) {
delete hills_energy;
@@ -1316,17 +1307,6 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
}
}
- is >> brace;
- if (brace != "}") {
- cvm::log("Incomplete restart information for metadynamics bias \""+
- this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- ": no closing brace at position "+
- cvm::to_str(is.tellg())+" in the file.\n");
- is.setstate(std::ios::failbit);
- return is;
- }
-
if (cvm::debug())
cvm::log("colvarbias_meta::read_restart() done\n");
@@ -1424,13 +1404,6 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
}
-
-
-// **********************************************************************
-// output functions
-// **********************************************************************
-
-
int colvarbias_meta::setup_output()
{
@@ -1537,16 +1510,17 @@ int colvarbias_meta::setup_output()
}
-std::ostream & colvarbias_meta::write_restart(std::ostream& os)
+std::string const colvarbias_meta::get_state_params() const
{
- os << "metadynamics {\n"
- << " configuration {\n"
- << " step " << cvm::step_absolute() << "\n"
- << " name " << this->name << "\n";
+ std::ostringstream os;
if (this->comm != single_replica)
- os << " replicaID " << this->replica_id << "\n";
- os << " }\n\n";
+ os << "replicaID " << this->replica_id << "\n";
+ return (colvarbias::get_state_params() + os.str());
+}
+
+std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
+{
if (use_grids) {
// this is a very good time to project hills, if you haven't done
@@ -1578,8 +1552,12 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
}
}
- os << "}\n\n";
+ return os;
+}
+
+int colvarbias_meta::write_state_to_replicas()
+{
if (comm != single_replica) {
write_replica_state_file();
// schedule to reread the state files of the other replicas (they
@@ -1588,12 +1566,16 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
(replicas[ir])->replica_state_file_in_sync = false;
}
}
+ return COLVARS_OK;
+}
+
+int colvarbias_meta::write_output_files()
+{
if (dump_fes) {
write_pmf();
}
-
- return os;
+ return COLVARS_OK;
}
@@ -1665,53 +1647,67 @@ void colvarbias_meta::write_pmf()
-void colvarbias_meta::write_replica_state_file()
+int colvarbias_meta::write_replica_state_file()
{
- // write down also the restart for the other replicas: TODO: this
- // is duplicated code, that could be cleaned up later
- cvm::backup_file(replica_state_file.c_str());
- cvm::ofstream rep_state_os(replica_state_file.c_str());
- if (!rep_state_os.is_open())
- cvm::fatal_error("Error: in opening file \""+
- replica_state_file+"\" for writing.\n");
-
- rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
- rep_state_os << "\n"
- << "metadynamics {\n"
- << " configuration {\n"
- << " name " << this->name << "\n"
- << " step " << cvm::step_absolute() << "\n";
- if (this->comm != single_replica) {
- rep_state_os << " replicaID " << this->replica_id << "\n";
- }
- rep_state_os << " }\n\n";
- rep_state_os << " hills_energy\n";
- rep_state_os << std::setprecision(cvm::cv_prec)
- << std::setw(cvm::cv_width);
- hills_energy->write_restart(rep_state_os);
- rep_state_os << " hills_energy_gradients\n";
- rep_state_os << std::setprecision(cvm::cv_prec)
- << std::setw(cvm::cv_width);
- hills_energy_gradients->write_restart(rep_state_os);
-
- if ( (!use_grids) || keep_hills ) {
- // write all hills currently in memory
- for (std::list<hill>::const_iterator h = this->hills.begin();
- h != this->hills.end();
- h++) {
- rep_state_os << *h;
- }
- } else {
- // write just those that are near the grid boundaries
- for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
- h != this->hills_off_grid.end();
- h++) {
- rep_state_os << *h;
- }
+ if (cvm::debug()) {
+ cvm::log("Writing replica state file for bias \""+name+"\"\n");
}
-
- rep_state_os << "}\n\n";
- rep_state_os.close();
+ // write down also the restart for the other replicas
+ cvm::backup_file(replica_state_file.c_str());
+ std::ostream *rep_state_os = cvm::proxy->output_stream(replica_state_file);
+ if (rep_state_os == NULL) {
+ cvm::error("Error: in opening file \""+
+ replica_state_file+"\" for writing.\n", FILE_ERROR);
+ return FILE_ERROR;
+ }
+
+ rep_state_os->setf(std::ios::scientific, std::ios::floatfield);
+
+ if (!write_state(*rep_state_os)) {
+ cvm::error("Error: in writing to file \""+
+ replica_state_file+"\".\n", FILE_ERROR);
+ cvm::proxy->close_output_stream(replica_state_file);
+ return FILE_ERROR;
+ }
+
+ cvm::proxy->close_output_stream(replica_state_file);
+
+ // rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
+ // rep_state_os << "\n"
+ // << "metadynamics {\n"
+ // << " configuration {\n"
+ // << " name " << this->name << "\n"
+ // << " step " << cvm::step_absolute() << "\n";
+ // if (this->comm != single_replica) {
+ // rep_state_os << " replicaID " << this->replica_id << "\n";
+ // }
+ // rep_state_os << " }\n\n";
+ // rep_state_os << " hills_energy\n";
+ // rep_state_os << std::setprecision(cvm::cv_prec)
+ // << std::setw(cvm::cv_width);
+ // hills_energy->write_restart(rep_state_os);
+ // rep_state_os << " hills_energy_gradients\n";
+ // rep_state_os << std::setprecision(cvm::cv_prec)
+ // << std::setw(cvm::cv_width);
+ // hills_energy_gradients->write_restart(rep_state_os);
+
+ // if ( (!use_grids) || keep_hills ) {
+ // // write all hills currently in memory
+ // for (std::list<hill>::const_iterator h = this->hills.begin();
+ // h != this->hills.end();
+ // h++) {
+ // rep_state_os << *h;
+ // }
+ // } else {
+ // // write just those that are near the grid boundaries
+ // for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
+ // h != this->hills_off_grid.end();
+ // h++) {
+ // rep_state_os << *h;
+ // }
+ // }
+ // rep_state_os << "}\n\n";
+ // rep_state_os.close();
// reopen the hills file
replica_hills_os.close();
@@ -1721,8 +1717,11 @@ void colvarbias_meta::write_replica_state_file()
cvm::fatal_error("Error: in opening file \""+
replica_hills_file+"\" for writing.\n");
replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
+
+ return COLVARS_OK;
}
+
std::string colvarbias_meta::hill::output_traj()
{
std::ostringstream os;
diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h
index 6936e9795454d702a8d825f8f2f0a278f28c1b42..a88a34ba00a7749895390cd81e5179687ebb5599 100644
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_META_H
#define COLVARBIAS_META_H
@@ -31,10 +38,16 @@ public:
virtual int init(std::string const &conf);
virtual ~colvarbias_meta();
virtual int update();
- virtual std::istream & read_restart(std::istream &is);
- virtual std::ostream & write_restart(std::ostream &os);
+
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &state_conf);
+ virtual std::ostream & write_state_data(std::ostream &os);
+ virtual std::istream & read_state_data(std::istream &os);
+
virtual int setup_output();
+ virtual int write_output_files();
virtual void write_pmf();
+ virtual int write_state_to_replicas();
class hill;
typedef std::list<hill>::iterator hill_iter;
@@ -77,13 +90,6 @@ protected:
/// Read a hill from a file
std::istream & read_hill(std::istream &is);
- /// \brief step present in a state file
- ///
- /// When using grids and reading state files containing them
- /// (multiple replicas), this is used to check whether a hill is
- /// newer or older than the grids
- size_t state_file_step;
-
/// \brief Add a new hill; if a .hills trajectory is written,
/// write it there; if there is more than one replica, communicate
/// it to the others
@@ -187,7 +193,7 @@ protected:
virtual void read_replica_files();
/// \brief Write data to other replicas
- virtual void write_replica_state_file();
+ virtual int write_replica_state_file();
/// \brief Additional, "mirror" metadynamics biases, to collect info
/// from the other replicas
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index bd37a430a1b7a088559ec3eb04890691cfcbda1f..84630984e5085fb9b07e96ca9ad8da9480608b62 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -1,10 +1,18 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias_restraint.h"
+
colvarbias_restraint::colvarbias_restraint(char const *key)
: colvarbias(key)
{
@@ -14,600 +22,1146 @@ colvarbias_restraint::colvarbias_restraint(char const *key)
int colvarbias_restraint::init(std::string const &conf)
{
colvarbias::init(conf);
+ enable(f_cvb_apply_force);
if (cvm::debug())
cvm::log("Initializing a new restraint bias.\n");
- // TODO move these initializations to constructor and let get_keyval
- // only override existing values
- target_nstages = 0;
- target_nsteps = 0;
- force_k = 0.0;
+ return COLVARS_OK;
+}
- get_keyval(conf, "forceConstant", force_k, 1.0);
- {
- // get the initial restraint centers
- colvar_centers.resize(colvars.size());
- colvar_centers_raw.resize(colvars.size());
- size_t i;
+int colvarbias_restraint::update()
+{
+ bias_energy = 0.0;
+
+ if (cvm::debug())
+ cvm::log("Updating the restraint bias \""+this->name+"\".\n");
+
+ // Force and energy calculation
+ for (size_t i = 0; i < colvars.size(); i++) {
+ colvar_forces[i].type(colvars[i]->value());
+ colvar_forces[i].is_derivative();
+ colvar_forces[i] = restraint_force(i);
+ bias_energy += restraint_potential(i);
+ }
+
+ if (cvm::debug())
+ cvm::log("Done updating the restraint bias \""+this->name+"\".\n");
+
+ if (cvm::debug())
+ cvm::log("Current forces for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(colvar_forces)+".\n");
+
+ return COLVARS_OK;
+}
+
+
+colvarbias_restraint::~colvarbias_restraint()
+{
+ if (cvm::n_rest_biases > 0)
+ cvm::n_rest_biases -= 1;
+}
+
+
+std::string const colvarbias_restraint::get_state_params() const
+{
+ return colvarbias::get_state_params();
+}
+
+
+int colvarbias_restraint::set_state_params(std::string const &conf)
+{
+ return colvarbias::set_state_params(conf);
+}
+
+
+std::ostream & colvarbias_restraint::write_traj_label(std::ostream &os)
+{
+ return colvarbias::write_traj_label(os);
+}
+
+
+std::ostream & colvarbias_restraint::write_traj(std::ostream &os)
+{
+ return colvarbias::write_traj(os);
+}
+
+
+
+colvarbias_restraint_centers::colvarbias_restraint_centers(char const *key)
+ : colvarbias(key), colvarbias_restraint(key)
+{
+}
+
- enable(f_cvb_apply_force);
+int colvarbias_restraint_centers::init(std::string const &conf)
+{
+ size_t i;
+ 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[i].reset();
colvar_centers_raw[i].type(colvars[i]->value());
+ colvar_centers_raw[i].reset();
+ }
+ }
+
+ if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
+ for (i = 0; i < colvars.size(); i++) {
if (cvm::debug()) {
- cvm::log("colvarbias_restraint: center size = "+
- cvm::to_str(colvar_centers[i].vector1d_value.size())+"\n");
+ cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n");
}
+ colvar_centers_raw[i] = colvar_centers[i];
+ colvar_centers[i].apply_constraints();
}
- if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
- for (i = 0; i < colvars.size(); i++) {
- if (cvm::debug()) {
- cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n");
- }
+ null_centers = false;
+ }
- colvar_centers[i].apply_constraints();
- colvar_centers_raw[i] = colvar_centers[i];
- }
- } else {
- colvar_centers.clear();
- cvm::error("Error: must define the initial centers of the restraints.\n");
+ if (null_centers) {
+ colvar_centers.clear();
+ cvm::error("Error: must define the initial centers of the restraints.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+
+ if (colvar_centers.size() != colvars.size()) {
+ cvm::error("Error: number of centers does not match "
+ "that of collective variables.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+
+ return COLVARS_OK;
+}
+
+
+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());
+ colvar_centers[i].apply_constraints();
+ colvar_centers_raw[i].type(colvars[i]->value());
+ colvar_centers_raw[i] = colvar_centers[i];
}
+ }
+ return COLVARS_OK;
+}
- if (colvar_centers.size() != colvars.size()) {
- cvm::error("Error: number of centers does not match "
- "that of collective variables.\n");
+
+
+colvarbias_restraint_k::colvarbias_restraint_k(char const *key)
+ : colvarbias(key), colvarbias_restraint(key)
+{
+ force_k = -1.0;
+}
+
+
+int colvarbias_restraint_k::init(std::string const &conf)
+{
+ get_keyval(conf, "forceConstant", force_k, (force_k > 0.0 ? force_k : 1.0));
+ if (force_k < 0.0) {
+ cvm::error("Error: undefined or invalid force constant.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+ return COLVARS_OK;
+}
+
+
+int colvarbias_restraint_k::change_configuration(std::string const &conf)
+{
+ get_keyval(conf, "forceConstant", force_k, force_k);
+ return COLVARS_OK;
+}
+
+
+
+colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key)
+{
+ target_nstages = 0;
+ target_nsteps = 0;
+ stage = 0;
+ b_chg_centers = false;
+ b_chg_force_k = false;
+}
+
+
+int colvarbias_restraint_moving::init(std::string const &conf)
+{
+ if (b_chg_centers && b_chg_force_k) {
+ cvm::error("Error: cannot specify both targetCenters and targetForceConstant.\n",
+ INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+
+ if (b_chg_centers || b_chg_force_k) {
+
+ get_keyval(conf, "targetNumSteps", target_nsteps, target_nsteps);
+ if (!target_nsteps) {
+ cvm::error("Error: targetNumSteps must be non-zero.\n", INPUT_ERROR);
+ return cvm::get_error();
+ }
+
+ if (get_keyval(conf, "targetNumStages", target_nstages, target_nstages) &&
+ lambda_schedule.size()) {
+ cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", INPUT_ERROR);
+ return cvm::get_error();
}
}
- {
- if (cvm::debug()) {
- cvm::log("colvarbias_restraint: parsing target centers.\n");
+ return COLVARS_OK;
+}
+
+
+std::string const colvarbias_restraint_moving::get_state_params() const
+{
+ std::ostringstream os;
+ os.setf(std::ios::scientific, std::ios::floatfield);
+ if (b_chg_centers || b_chg_force_k) {
+ // TODO move this
+ if (target_nstages) {
+ os << "stage " << std::setw(cvm::it_width)
+ << stage << "\n";
}
+ }
+ return os.str();
+}
- size_t i;
- if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
- if (colvar_centers.size() != colvars.size()) {
- cvm::error("Error: number of target centers does not match "
- "that of collective variables.\n");
- }
- b_chg_centers = true;
- for (i = 0; i < target_centers.size(); i++) {
- target_centers[i].apply_constraints();
- }
- } else {
- b_chg_centers = false;
- target_centers.clear();
+
+int colvarbias_restraint_moving::set_state_params(std::string const &conf)
+{
+ if (b_chg_centers || b_chg_force_k) {
+ if (target_nstages) {
+ // cvm::log ("Reading current stage from the restart.\n");
+ if (!get_keyval(conf, "stage", stage))
+ cvm::error("Error: current stage is missing from the restart.\n");
}
}
+ return COLVARS_OK;
+}
- if (get_keyval(conf, "targetForceConstant", target_force_k, 0.0)) {
- if (b_chg_centers)
- cvm::error("Error: cannot specify both targetCenters and targetForceConstant.\n");
- starting_force_k = force_k;
- b_chg_force_k = true;
- get_keyval(conf, "targetEquilSteps", target_equil_steps, 0);
+colvarbias_restraint_centers_moving::colvarbias_restraint_centers_moving(char const *key)
+ : colvarbias(key),
+ colvarbias_restraint(key),
+ colvarbias_restraint_centers(key),
+ colvarbias_restraint_moving(key)
+{
+ b_chg_centers = false;
+ b_output_centers = false;
+ b_output_acc_work = false;
+ acc_work = 0.0;
+}
+
+
+int colvarbias_restraint_centers_moving::init(std::string const &conf)
+{
+ colvarbias_restraint_centers::init(conf);
- get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule);
- if (lambda_schedule.size()) {
- // There is one more lambda-point than stages
- target_nstages = lambda_schedule.size() - 1;
+ if (cvm::debug()) {
+ cvm::log("colvarbias_restraint: parsing target centers.\n");
+ }
+
+ size_t i;
+ if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
+ if (colvar_centers.size() != colvars.size()) {
+ cvm::error("Error: number of target centers does not match "
+ "that of collective variables.\n");
}
+ b_chg_centers = true;
+ for (i = 0; i < target_centers.size(); i++) {
+ target_centers[i].apply_constraints();
+ }
+ }
+
+ if (b_chg_centers) {
+ // parse moving restraint options
+ colvarbias_restraint_moving::init(conf);
} else {
- b_chg_force_k = false;
+ target_centers.clear();
+ return COLVARS_OK;
}
- if (b_chg_centers || b_chg_force_k) {
- get_keyval(conf, "targetNumSteps", target_nsteps, 0);
- if (!target_nsteps)
- cvm::error("Error: targetNumSteps must be non-zero.\n");
+ get_keyval(conf, "outputCenters", b_output_centers, b_output_centers);
+ get_keyval(conf, "outputAccumulatedWork", b_output_acc_work, b_output_acc_work);
- if (get_keyval(conf, "targetNumStages", target_nstages, target_nstages) &&
- lambda_schedule.size()) {
- cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n");
+ return COLVARS_OK;
+}
+
+
+int colvarbias_restraint_centers_moving::update()
+{
+ if (b_chg_centers) {
+
+ if (cvm::debug()) {
+ cvm::log("Updating centers for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
+ }
+
+ if (!centers_incr.size()) {
+ // if this is the first calculation, calculate the advancement
+ // 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[i] = (target_centers[i] - colvar_centers_raw[i]) /
+ cvm::real( target_nstages ? (target_nstages - stage) :
+ (target_nsteps - cvm::step_absolute()));
+ }
+ if (cvm::debug()) {
+ cvm::log("Center increment for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(centers_incr)+" at stage "+cvm::to_str(stage)+ ".\n");
+ }
}
if (target_nstages) {
- // This means that either numStages of lambdaSchedule has been provided
- stage = 0;
- restraint_FE = 0.0;
+ if ((cvm::step_relative() > 0)
+ && (cvm::step_absolute() % target_nsteps) == 0
+ && stage < target_nstages) {
+
+ for (size_t i = 0; i < colvars.size(); i++) {
+ colvar_centers_raw[i] += centers_incr[i];
+ colvar_centers[i] = colvar_centers_raw[i];
+ colvars[i]->wrap(colvar_centers[i]);
+ colvar_centers[i].apply_constraints();
+ }
+ stage++;
+ cvm::log("Moving restraint \"" + this->name +
+ "\" stage " + cvm::to_str(stage) +
+ " : setting centers to " + cvm::to_str(colvar_centers) +
+ " at step " + cvm::to_str(cvm::step_absolute()));
+ }
+ } 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++) {
+ colvar_centers_raw[i] += centers_incr[i];
+ colvar_centers[i] = colvar_centers_raw[i];
+ colvars[i]->wrap(colvar_centers[i]);
+ colvar_centers[i].apply_constraints();
+ }
+ }
+
+ if (cvm::debug()) {
+ cvm::log("New centers for the restraint bias \""+
+ this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
}
+ }
+
+ return COLVARS_OK;
+}
- if (get_keyval(conf, "targetForceExponent", force_k_exp, 1.0)) {
- if (! b_chg_force_k)
- cvm::log("Warning: not changing force constant: targetForceExponent will be ignored\n");
- if (force_k_exp < 1.0)
- cvm::log("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n");
+
+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++) {
+ // project forces on the calculated increments at this step
+ acc_work += colvar_forces[i] * centers_incr[i];
+ }
}
}
+ return COLVARS_OK;
+}
+
+
+std::string const colvarbias_restraint_centers_moving::get_state_params() const
+{
+ std::ostringstream os;
+ os.setf(std::ios::scientific, std::ios::floatfield);
- get_keyval(conf, "outputCenters", b_output_centers, false);
if (b_chg_centers) {
- get_keyval(conf, "outputAccumulatedWork", b_output_acc_work, false);
- } else {
- b_output_acc_work = false;
+ size_t i;
+ os << "centers ";
+ for (i = 0; i < colvars.size(); 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++) {
+ os << " "
+ << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+ << colvar_centers_raw[i];
+ }
+ os << "\n";
+
+ if (b_output_acc_work) {
+ os << "accumulatedWork "
+ << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+ << acc_work << "\n";
+ }
}
- acc_work = 0.0;
- if (cvm::debug())
- cvm::log("Done initializing a new restraint bias.\n");
+ return colvarbias_restraint_moving::get_state_params() + os.str();
+}
+
+
+int colvarbias_restraint_centers_moving::set_state_params(std::string const &conf)
+{
+ colvarbias_restraint::set_state_params(conf);
+
+ if (b_chg_centers) {
+ // cvm::log ("Reading the updated restraint centers from the restart.\n");
+ if (!get_keyval(conf, "centers", colvar_centers))
+ cvm::error("Error: restraint centers are missing from the restart.\n");
+ if (!get_keyval(conf, "centers_raw", colvar_centers_raw))
+ cvm::error("Error: \"raw\" restraint centers are missing from the restart.\n");
+ if (b_output_acc_work) {
+ if (!get_keyval(conf, "accumulatedWork", acc_work))
+ cvm::error("Error: accumulatedWork is missing from the restart.\n");
+ }
+ }
return COLVARS_OK;
}
-colvarbias_restraint::~colvarbias_restraint()
+std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os)
{
- if (cvm::n_rest_biases > 0)
- cvm::n_rest_biases -= 1;
+ 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);
+ os << " x0_"
+ << cvm::wrap_string(colvars[i]->name, this_cv_width-3);
+ }
+ }
+
+ if (b_output_acc_work) {
+ os << " W_"
+ << cvm::wrap_string(this->name, cvm::en_width-2);
+ }
+
+ return os;
}
-void colvarbias_restraint::change_configuration(std::string const &conf)
+std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
{
- get_keyval(conf, "forceConstant", force_k, force_k);
- if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
+ if (b_output_centers) {
for (size_t i = 0; i < colvars.size(); i++) {
- colvar_centers[i].type(colvars[i]->value());
- colvar_centers[i].apply_constraints();
- colvar_centers_raw[i].type(colvars[i]->value());
- colvar_centers_raw[i] = colvar_centers[i];
+ os << " "
+ << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+ << colvar_centers[i];
}
}
+
+ if (b_output_acc_work) {
+ os << " "
+ << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+ << acc_work;
+ }
+
+ return os;
}
-cvm::real colvarbias_restraint::energy_difference(std::string const &conf)
+
+colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key)
+ : colvarbias(key),
+ colvarbias_restraint(key),
+ colvarbias_restraint_k(key),
+ colvarbias_restraint_moving(key)
{
- std::vector<colvarvalue> alt_colvar_centers;
- cvm::real alt_force_k;
- cvm::real alt_bias_energy = 0.0;
+ b_chg_force_k = false;
+ target_equil_steps = 0;
+ target_force_k = 0.0;
+ starting_force_k = 0.0;
+ force_k_exp = 1.0;
+ restraint_FE = 0.0;
+}
- get_keyval(conf, "forceConstant", alt_force_k, force_k);
- alt_colvar_centers.resize(colvars.size());
- size_t i;
- for (i = 0; i < colvars.size(); i++) {
- alt_colvar_centers[i].type(colvars[i]->value());
+int colvarbias_restraint_k_moving::init(std::string const &conf)
+{
+ colvarbias_restraint_k::init(conf);
+
+ if (get_keyval(conf, "targetForceConstant", target_force_k, target_force_k)) {
+ starting_force_k = force_k;
+ b_chg_force_k = true;
}
- if (get_keyval(conf, "centers", alt_colvar_centers, colvar_centers)) {
- for (i = 0; i < colvars.size(); i++) {
- alt_colvar_centers[i].apply_constraints();
+
+ if (b_chg_force_k) {
+ // parse moving restraint options
+ colvarbias_restraint_moving::init(conf);
+ } else {
+ return COLVARS_OK;
+ }
+
+ get_keyval(conf, "targetEquilSteps", target_equil_steps, target_equil_steps);
+
+ if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) &&
+ target_nstages > 0) {
+ cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", INPUT_ERROR);
+ return cvm::get_error();
+ }
+
+ if (lambda_schedule.size()) {
+ // There is one more lambda-point than stages
+ target_nstages = lambda_schedule.size() - 1;
+ }
+
+ if (get_keyval(conf, "targetForceExponent", force_k_exp, force_k_exp)) {
+ if (! b_chg_force_k)
+ cvm::log("Warning: not changing force constant: targetForceExponent will be ignored\n");
+ }
+ if (force_k_exp < 1.0) {
+ cvm::log("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n");
+ }
+
+ return COLVARS_OK;
+}
+
+
+int colvarbias_restraint_k_moving::update()
+{
+ if (b_chg_force_k) {
+
+ cvm::real lambda;
+
+ if (target_nstages) {
+
+ if (cvm::step_absolute() == 0) {
+ // Setup first stage of staged variable force constant calculation
+ if (lambda_schedule.size()) {
+ lambda = lambda_schedule[0];
+ } else {
+ lambda = 0.0;
+ }
+ 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));
+ }
+
+ // TI calculation: estimate free energy derivative
+ // need current lambda
+ if (lambda_schedule.size()) {
+ lambda = lambda_schedule[stage];
+ } else {
+ lambda = cvm::real(stage) / cvm::real(target_nstages);
+ }
+
+ if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) {
+ // Start averaging after equilibration period, if requested
+
+ // Square distance normalized by square colvar width
+ cvm::real dist_sq = 0.0;
+ for (size_t i = 0; i < colvars.size(); i++) {
+ dist_sq += d_restraint_potential_dk(i);
+ }
+
+ restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0)
+ * (target_force_k - starting_force_k) * dist_sq;
+ }
+
+ // Finish current stage...
+ if (cvm::step_absolute() % target_nsteps == 0 &&
+ cvm::step_absolute() > 0) {
+
+ cvm::log("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
+ if (stage < target_nstages) {
+
+ restraint_FE = 0.0;
+ stage++;
+ if (lambda_schedule.size()) {
+ lambda = lambda_schedule[stage];
+ } else {
+ lambda = cvm::real(stage) / cvm::real(target_nstages);
+ }
+ 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));
+ }
+ }
+
+ } else if (cvm::step_absolute() <= target_nsteps) {
+
+ // update force constant (slow growth)
+ lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
+ force_k = starting_force_k + (target_force_k - starting_force_k)
+ * std::pow(lambda, force_k_exp);
}
}
- for (i = 0; i < colvars.size(); i++) {
- alt_bias_energy += restraint_potential(restraint_convert_k(alt_force_k, colvars[i]->width),
- colvars[i],
- alt_colvar_centers[i]);
- }
+ return COLVARS_OK;
+}
+
+
+std::string const colvarbias_restraint_k_moving::get_state_params() const
+{
+ std::ostringstream os;
+ os.setf(std::ios::scientific, std::ios::floatfield);
+ if (b_chg_force_k) {
+ os << "forceConstant "
+ << std::setprecision(cvm::en_prec)
+ << std::setw(cvm::en_width) << force_k << "\n";
+ }
+ return colvarbias_restraint_moving::get_state_params() + os.str();
+}
+
+
+int colvarbias_restraint_k_moving::set_state_params(std::string const &conf)
+{
+ colvarbias_restraint::set_state_params(conf);
+
+ if (b_chg_force_k) {
+ // cvm::log ("Reading the updated force constant from the restart.\n");
+ if (!get_keyval(conf, "forceConstant", force_k, force_k))
+ cvm::error("Error: force constant is missing from the restart.\n");
+ }
+
+ return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias_restraint_k_moving::write_traj_label(std::ostream &os)
+{
+ return os;
+}
+
+
+std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os)
+{
+ return os;
+}
+
+
+
+// redefined due to legacy state file keyword "harmonic"
+std::istream & colvarbias_restraint::read_state(std::istream &is)
+{
+ size_t const start_pos = is.tellg();
+
+ std::string key, brace, conf;
+ if ( !(is >> key) || !(key == "restraint" || key == "harmonic") ||
+ !(is >> brace) || !(brace == "{") ||
+ !(is >> colvarparse::read_block("configuration", conf)) ||
+ (set_state_params(conf) != COLVARS_OK) ) {
+ cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
+ this->name+"\" at position "+
+ cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+ is.clear();
+ is.seekg(start_pos, std::ios::beg);
+ is.setstate(std::ios::failbit);
+ return is;
+ }
+
+ if (!read_state_data(is)) {
+ cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
+ this->name+"\" at position "+
+ cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+ is.clear();
+ is.seekg(start_pos, std::ios::beg);
+ is.setstate(std::ios::failbit);
+ }
+
+ is >> brace;
+ if (brace != "}") {
+ cvm::log("brace = "+brace+"\n");
+ cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
+ this->name+"\": no matching brace at position "+
+ cvm::to_str(is.tellg())+" in stream.\n");
+ is.setstate(std::ios::failbit);
+ }
+
+ return is;
+}
+
+
+std::ostream & colvarbias_restraint::write_state(std::ostream &os)
+{
+ os.setf(std::ios::scientific, std::ios::floatfield);
+ os << "restraint {\n"
+ << " configuration {\n";
+ std::istringstream is(get_state_params());
+ std::string line;
+ while (std::getline(is, line)) {
+ os << " " << line << "\n";
+ }
+ os << " }\n";
+ write_state_data(os);
+ os << "}\n\n";
+ return os;
+}
+
+
+
+colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
+ : colvarbias(key),
+ colvarbias_restraint(key),
+ colvarbias_restraint_centers(key),
+ colvarbias_restraint_moving(key),
+ colvarbias_restraint_k(key),
+ colvarbias_restraint_centers_moving(key),
+ colvarbias_restraint_k_moving(key)
+{
+}
+
+
+int colvarbias_restraint_harmonic::init(std::string const &conf)
+{
+ colvarbias_restraint::init(conf);
+ colvarbias_restraint_moving::init(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+
+ "\" will be rescaled to "+
+ cvm::to_str(force_k / (colvars[i]->width * colvars[i]->width))+
+ " according to the specified width.\n");
+ }
+
+ return COLVARS_OK;
+}
+
+
+int colvarbias_restraint_harmonic::update()
+{
+ // update parameters (centers or force constant)
+ colvarbias_restraint_centers_moving::update();
+ colvarbias_restraint_k_moving::update();
+
+ // update restraint energy and forces
+ colvarbias_restraint::update();
+
+ // update accumulated work using the current forces
+ colvarbias_restraint_centers_moving::update_acc_work();
+
+ return COLVARS_OK;
+}
+
+
+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]);
+}
+
+
+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]);
+}
+
+
+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]);
+}
+
+
+std::string const colvarbias_restraint_harmonic::get_state_params() const
+{
+ return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_centers_moving::get_state_params() +
+ colvarbias_restraint_k_moving::get_state_params();
+}
+
+
+int colvarbias_restraint_harmonic::set_state_params(std::string const &conf)
+{
+ int error_code = COLVARS_OK;
+ error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
+ error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
+ return error_code;
+}
- return alt_bias_energy - bias_energy;
+
+std::ostream & colvarbias_restraint_harmonic::write_traj_label(std::ostream &os)
+{
+ colvarbias_restraint::write_traj_label(os);
+ colvarbias_restraint_centers_moving::write_traj_label(os);
+ colvarbias_restraint_k_moving::write_traj_label(os);
+ return os;
}
-int colvarbias_restraint::update()
+std::ostream & colvarbias_restraint_harmonic::write_traj(std::ostream &os)
{
- bias_energy = 0.0;
+ colvarbias_restraint::write_traj(os);
+ colvarbias_restraint_centers_moving::write_traj(os);
+ colvarbias_restraint_k_moving::write_traj(os);
+ return os;
+}
- if (cvm::debug())
- cvm::log("Updating the restraint bias \""+this->name+"\".\n");
- // Setup first stage of staged variable force constant calculation
- if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) {
- cvm::real lambda;
- if (lambda_schedule.size()) {
- lambda = lambda_schedule[0];
- } else {
- lambda = 0.0;
- }
- 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));
- }
+int colvarbias_restraint_harmonic::change_configuration(std::string const &conf)
+{
+ return colvarbias_restraint_centers::change_configuration(conf) |
+ colvarbias_restraint_k::change_configuration(conf);
+}
- if (b_chg_centers) {
- if (!centers_incr.size()) {
- // if this is the first calculation, calculate the advancement
- // 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[i] = (target_centers[i] - colvar_centers_raw[i]) /
- cvm::real( target_nstages ? (target_nstages - stage) :
- (target_nsteps - cvm::step_absolute()));
- }
- if (cvm::debug()) {
- cvm::log("Center increment for the restraint bias \""+
- this->name+"\": "+cvm::to_str(centers_incr)+" at stage "+cvm::to_str(stage)+ ".\n");
- }
- }
+cvm::real colvarbias_restraint_harmonic::energy_difference(std::string const &conf)
+{
+ cvm::real const old_bias_energy = bias_energy;
+ cvm::real const old_force_k = force_k;
+ std::vector<colvarvalue> const old_centers = colvar_centers;
- if (target_nstages) {
- if ((cvm::step_relative() > 0)
- && (cvm::step_absolute() % target_nsteps) == 0
- && stage < target_nstages) {
-
- for (size_t i = 0; i < colvars.size(); i++) {
- colvar_centers_raw[i] += centers_incr[i];
- colvar_centers[i] = colvar_centers_raw[i];
- colvars[i]->wrap(colvar_centers[i]);
- colvar_centers[i].apply_constraints();
- }
- stage++;
- cvm::log("Moving restraint \"" + this->name +
- "\" stage " + cvm::to_str(stage) +
- " : setting centers to " + cvm::to_str(colvar_centers) +
- " at step " + cvm::to_str(cvm::step_absolute()));
- }
- } 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++) {
- colvar_centers_raw[i] += centers_incr[i];
- colvar_centers[i] = colvar_centers_raw[i];
- colvars[i]->wrap(colvar_centers[i]);
- colvar_centers[i].apply_constraints();
- }
- }
+ change_configuration(conf);
+ update();
- if (cvm::debug())
- cvm::log("Current centers for the restraint bias \""+
- this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
- }
+ cvm::real const result = (bias_energy - old_bias_energy);
- if (b_chg_force_k) {
- // Coupling parameter, between 0 and 1
- cvm::real lambda;
+ bias_energy = old_bias_energy;
+ force_k = old_force_k;
+ colvar_centers = old_centers;
- if (target_nstages) {
- // TI calculation: estimate free energy derivative
- // need current lambda
- if (lambda_schedule.size()) {
- lambda = lambda_schedule[stage];
- } else {
- lambda = cvm::real(stage) / cvm::real(target_nstages);
- }
+ return result;
+}
- if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) {
- // Start averaging after equilibration period, if requested
- // Square distance normalized by square colvar width
- cvm::real dist_sq = 0.0;
- for (size_t i = 0; i < colvars.size(); i++) {
- dist_sq += colvars[i]->dist2(colvars[i]->value(), colvar_centers[i])
- / (colvars[i]->width * colvars[i]->width);
- }
- restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0)
- * (target_force_k - starting_force_k) * dist_sq;
- }
+colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char const *key)
+ : colvarbias(key),
+ colvarbias_restraint(key),
+ colvarbias_restraint_k(key),
+ colvarbias_restraint_moving(key),
+ colvarbias_restraint_k_moving(key)
+{
+}
- // Finish current stage...
- if (cvm::step_absolute() % target_nsteps == 0 &&
- cvm::step_absolute() > 0) {
- cvm::log("Lambda= " + cvm::to_str(lambda) + " dA/dLambda= "
- + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps)));
+int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
+{
+ colvarbias_restraint::init(conf);
+ colvarbias_restraint_moving::init(conf);
+ colvarbias_restraint_k_moving::init(conf);
- // ...and move on to the next one
- if (stage < target_nstages) {
+ provide(f_cvb_scalar_variables);
+ enable(f_cvb_scalar_variables);
- restraint_FE = 0.0;
- stage++;
- if (lambda_schedule.size()) {
- lambda = lambda_schedule[stage];
- } else {
- lambda = cvm::real(stage) / cvm::real(target_nstages);
- }
- 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));
- }
- }
- } else if (cvm::step_absolute() <= target_nsteps) {
- // update force constant (slow growth)
- lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
- force_k = starting_force_k + (target_force_k - starting_force_k)
- * std::pow(lambda, force_k_exp);
+ size_t i;
+
+ 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[i].reset();
}
}
+ if (!get_keyval(conf, "lowerWalls", lower_walls, lower_walls) &&
+ b_null_lower_walls) {
+ cvm::log("Lower walls were not provided.\n");
+ lower_walls.resize(0);
+ }
- if (cvm::debug())
- cvm::log("Done updating the restraint bias \""+this->name+"\".\n");
+ 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[i].reset();
+ }
+ }
+ if (!get_keyval(conf, "upperWalls", upper_walls, upper_walls) &&
+ b_null_upper_walls) {
+ cvm::log("Upper walls were not provided.\n");
+ upper_walls.resize(0);
+ }
- // Force and energy calculation
- for (size_t i = 0; i < colvars.size(); i++) {
- colvar_forces[i].type(colvars[i]->value());
- colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(force_k, colvars[i]->width),
- colvars[i],
- colvar_centers[i]);
- bias_energy += restraint_potential(restraint_convert_k(force_k, colvars[i]->width),
- colvars[i],
- colvar_centers[i]);
- if (cvm::debug()) {
- cvm::log("dist_grad["+cvm::to_str(i)+
- "] = "+cvm::to_str(colvars[i]->dist2_lgrad(colvars[i]->value(),
- colvar_centers[i]))+"\n");
+ if ((lower_walls.size() == 0) && (upper_walls.size() == 0)) {
+ cvm::error("Error: no walls provided.\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
+
+ if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) {
+ for (i = 0; i < colvars.size(); i++) {
+ if (colvars[i]->is_enabled(f_cv_periodic)) {
+ cvm::error("Error: at least one variable is periodic, "
+ "both walls must be provided .\n", INPUT_ERROR);
+ return INPUT_ERROR;
+ }
}
}
- if (b_output_acc_work) {
- if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) {
- for (size_t i = 0; i < colvars.size(); i++) {
- // project forces on the calculated increments at this step
- acc_work += colvar_forces[i] * centers_incr[i];
+ if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) {
+ for (i = 0; i < colvars.size(); i++) {
+ if (lower_walls[i] >= upper_walls[i]) {
+ cvm::error("Error: one upper wall, "+
+ cvm::to_str(upper_walls[i])+
+ ", is not higher than the lower wall, "+
+ cvm::to_str(lower_walls[i])+".\n",
+ INPUT_ERROR);
}
}
}
- if (cvm::debug())
- cvm::log("Current forces for the restraint bias \""+
- this->name+"\": "+cvm::to_str(colvar_forces)+".\n");
+ for (i = 0; i < colvars.size(); i++) {
+ if (colvars[i]->width != 1.0)
+ cvm::log("The force constant for colvar \""+colvars[i]->name+
+ "\" will be rescaled to "+
+ cvm::to_str(force_k / (colvars[i]->width * colvars[i]->width))+
+ " according to the specified width.\n");
+ }
return COLVARS_OK;
}
-std::istream & colvarbias_restraint::read_restart(std::istream &is)
+int colvarbias_restraint_harmonic_walls::update()
{
- size_t const start_pos = is.tellg();
-
- cvm::log("Restarting restraint bias \""+
- this->name+"\".\n");
-
- std::string key, brace, conf;
- if ( !(is >> key) || !(key == "restraint" || key == "harmonic") ||
- !(is >> brace) || !(brace == "{") ||
- !(is >> colvarparse::read_block("configuration", conf)) ) {
-
- cvm::log("Error: in reading restart configuration for restraint bias \""+
- this->name+"\" at position "+
- cvm::to_str(is.tellg())+" in stream.\n");
- is.clear();
- is.seekg(start_pos, std::ios::beg);
- is.setstate(std::ios::failbit);
- return is;
- }
-
-// int id = -1;
- std::string name = "";
-// if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) &&
-// (id != this->id) ) ||
- if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) &&
- (name != this->name) )
- cvm::error("Error: in the restart file, the "
- "\"restraint\" block has a wrong name\n");
-// if ( (id == -1) && (name == "") ) {
- if (name.size() == 0) {
- cvm::error("Error: \"restraint\" block in the restart file "
- "has no identifiers.\n");
- }
-
- if (b_chg_centers) {
-// cvm::log ("Reading the updated restraint centers from the restart.\n");
- if (!get_keyval(conf, "centers", colvar_centers))
- cvm::error("Error: restraint centers are missing from the restart.\n");
- if (!get_keyval(conf, "centers_raw", colvar_centers_raw))
- cvm::error("Error: \"raw\" restraint centers are missing from the restart.\n");
- }
+ colvarbias_restraint_k_moving::update();
- if (b_chg_force_k) {
-// cvm::log ("Reading the updated force constant from the restart.\n");
- if (!get_keyval(conf, "forceConstant", force_k))
- cvm::error("Error: force constant is missing from the restart.\n");
- }
+ colvarbias_restraint::update();
- if (target_nstages) {
-// cvm::log ("Reading current stage from the restart.\n");
- if (!get_keyval(conf, "stage", stage))
- cvm::error("Error: current stage is missing from the restart.\n");
- }
+ return COLVARS_OK;
+}
- if (b_output_acc_work) {
- if (!get_keyval(conf, "accumulatedWork", acc_work))
- cvm::error("Error: accumulatedWork is missing from the restart.\n");
- }
- is >> brace;
- if (brace != "}") {
- cvm::error("Error: corrupt restart information for restraint bias \""+
- this->name+"\": no matching brace at position "+
- cvm::to_str(is.tellg())+" in the restart file.\n");
- is.setstate(std::ios::failbit);
+void colvarbias_restraint_harmonic_walls::communicate_forces()
+{
+ for (size_t i = 0; i < colvars.size(); i++) {
+ if (cvm::debug()) {
+ cvm::log("Communicating a force to colvar \""+
+ colvars[i]->name+"\".\n");
+ }
+ colvars[i]->add_bias_force_actual_value(colvar_forces[i]);
}
- return is;
}
-std::ostream & colvarbias_restraint::write_restart(std::ostream &os)
+cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
{
- os << "restraint {\n"
- << " configuration {\n"
- // << " id " << this->id << "\n"
- << " name " << this->name << "\n";
-
- if (b_chg_centers) {
- size_t i;
- os << " centers ";
- for (i = 0; i < colvars.size(); i++) {
- os << " " << colvar_centers[i];
- }
- os << "\n";
- os << " centers_raw ";
- for (i = 0; i < colvars.size(); i++) {
- os << " " << colvar_centers_raw[i];
+ colvar *cv = colvars[i];
+ colvarvalue const &cvv = colvars[i]->actual_value();
+
+ // For a periodic colvar, both walls may be applicable at the same time
+ // in which case we pick the closer one
+
+ if (cv->is_enabled(f_cv_periodic)) {
+ cvm::real const lower_wall_dist2 = cv->dist2(cvv, lower_walls[i]);
+ 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; }
+ } else {
+ cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
+ if (grad > 0.0) { return grad; }
}
- os << "\n";
- }
-
- if (b_chg_force_k) {
- os << " forceConstant "
- << std::setprecision(cvm::en_prec)
- << std::setw(cvm::en_width) << force_k << "\n";
+ return 0.0;
}
- if (target_nstages) {
- os << " stage " << std::setw(cvm::it_width)
- << stage << "\n";
+ if (lower_walls.size() > 0) {
+ cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]);
+ if (grad < 0.0) { return grad; }
}
-
- if (b_output_acc_work) {
- os << " accumulatedWork " << acc_work << "\n";
+ if (upper_walls.size() > 0) {
+ cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
+ if (grad > 0.0) { return grad; }
}
+ return 0.0;
+}
- os << " }\n"
- << "}\n\n";
- return os;
+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) *
+ dist * dist;
}
-std::ostream & colvarbias_restraint::write_traj_label(std::ostream &os)
+colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const
{
- os << " ";
+ cvm::real const dist = colvar_distance(i);
+ return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
+ dist;
+}
- if (b_output_energy)
- os << " E_"
- << cvm::wrap_string(this->name, cvm::en_width-2);
- 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);
- os << " x0_"
- << cvm::wrap_string(colvars[i]->name, this_cv_width-3);
- }
+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) *
+ dist * dist;
+}
- if (b_output_acc_work)
- os << " W_"
- << cvm::wrap_string(this->name, cvm::en_width-2);
- return os;
+std::string const colvarbias_restraint_harmonic_walls::get_state_params() const
+{
+ return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_k_moving::get_state_params();
}
-std::ostream & colvarbias_restraint::write_traj(std::ostream &os)
+int colvarbias_restraint_harmonic_walls::set_state_params(std::string const &conf)
{
- os << " ";
+ int error_code = COLVARS_OK;
+ error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
+ return error_code;
+}
- if (b_output_energy)
- os << " "
- << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
- << bias_energy;
- if (b_output_centers)
- for (size_t i = 0; i < colvars.size(); i++) {
- os << " "
- << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
- << colvar_centers[i];
- }
+std::ostream & colvarbias_restraint_harmonic_walls::write_traj_label(std::ostream &os)
+{
+ colvarbias_restraint::write_traj_label(os);
+ colvarbias_restraint_k_moving::write_traj_label(os);
+ return os;
+}
- if (b_output_acc_work)
- os << " "
- << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
- << acc_work;
+std::ostream & colvarbias_restraint_harmonic_walls::write_traj(std::ostream &os)
+{
+ colvarbias_restraint::write_traj(os);
+ colvarbias_restraint_k_moving::write_traj(os);
return os;
}
-colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
- : colvarbias_restraint(key)
+
+colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
+ : colvarbias(key),
+ colvarbias_restraint(key),
+ colvarbias_restraint_centers(key),
+ colvarbias_restraint_moving(key),
+ colvarbias_restraint_k(key),
+ colvarbias_restraint_centers_moving(key),
+ colvarbias_restraint_k_moving(key)
{
}
-int colvarbias_restraint_harmonic::init(std::string const &conf)
+int colvarbias_restraint_linear::init(std::string const &conf)
{
colvarbias_restraint::init(conf);
+ colvarbias_restraint_moving::init(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)) {
+ 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+
"\" will be rescaled to "+
- cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
+ cvm::to_str(force_k / colvars[i]->width)+
" according to the specified width.\n");
}
+
return COLVARS_OK;
}
-cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k,
- colvar const *x,
- colvarvalue const &xcenter) const
+int colvarbias_restraint_linear::update()
{
- return 0.5 * k * x->dist2(x->value(), xcenter);
+ // update parameters (centers or force constant)
+ colvarbias_restraint_centers_moving::update();
+ colvarbias_restraint_k_moving::update();
+
+ // update restraint energy and forces
+ colvarbias_restraint::update();
+
+ // update accumulated work using the current forces
+ colvarbias_restraint_centers_moving::update_acc_work();
+
+ return COLVARS_OK;
}
-colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k,
- colvar const *x,
- colvarvalue const &xcenter) const
+int colvarbias_restraint_linear::change_configuration(std::string const &conf)
{
- return 0.5 * k * x->dist2_lgrad(x->value(), xcenter);
+ // Only makes sense to change the force constant
+ return colvarbias_restraint_k::change_configuration(conf);
}
-cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k,
- cvm::real dist_measure) const
+cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf)
{
- return k / (dist_measure * dist_measure);
+ cvm::real const old_bias_energy = bias_energy;
+ cvm::real const old_force_k = force_k;
+
+ change_configuration(conf);
+ update();
+
+ cvm::real const result = (bias_energy - old_bias_energy);
+
+ bias_energy = old_bias_energy;
+ force_k = old_force_k;
+
+ return result;
}
+cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
+{
+ return force_k / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]);
+}
+
-colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
- : colvarbias_restraint(key)
+colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
{
+ return -1.0 * force_k / colvars[i]->width;
}
-int colvarbias_restraint_linear::init(std::string const &conf)
+cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const
{
- colvarbias_restraint::init(conf);
+ return 1.0 / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]);
+}
- 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+
- "\" will be rescaled to "+
- cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
- " according to the specified width.\n");
- }
- return COLVARS_OK;
+
+std::string const colvarbias_restraint_linear::get_state_params() const
+{
+ return colvarbias_restraint::get_state_params() +
+ colvarbias_restraint_centers_moving::get_state_params() +
+ colvarbias_restraint_k_moving::get_state_params();
}
-cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k,
- colvar const *x,
- colvarvalue const &xcenter) const
+int colvarbias_restraint_linear::set_state_params(std::string const &conf)
{
- return k * (x->value() - xcenter);
+ int error_code = COLVARS_OK;
+ error_code |= colvarbias_restraint::set_state_params(conf);
+ error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
+ error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
+ return error_code;
}
-colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k,
- colvar const *x,
- colvarvalue const &xcenter) const
+std::ostream & colvarbias_restraint_linear::write_traj_label(std::ostream &os)
{
- return k;
+ colvarbias_restraint::write_traj_label(os);
+ colvarbias_restraint_centers_moving::write_traj_label(os);
+ colvarbias_restraint_k_moving::write_traj_label(os);
+ return os;
}
-cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,
- cvm::real dist_measure) const
+std::ostream & colvarbias_restraint_linear::write_traj(std::ostream &os)
{
- return k / dist_measure;
+ colvarbias_restraint::write_traj(os);
+ colvarbias_restraint_centers_moving::write_traj(os);
+ colvarbias_restraint_k_moving::write_traj(os);
+ return os;
}
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index 54097bd5240907d07565b1c57a191f96eba7b858..cbdd9c44d163a3c4684eafafaeb07cfa0a9b0495 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -1,36 +1,43 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARBIAS_RESTRAINT_H
#define COLVARBIAS_RESTRAINT_H
#include "colvarbias.h"
-/// \brief Bias restraint, optionally moving towards a target
+/// \brief Most general definition of a colvar restraint:
+/// see derived classes for specific types
/// (implementation of \link colvarbias \endlink)
-class colvarbias_restraint : public colvarbias {
+class colvarbias_restraint
+ : public virtual colvarbias
+{
public:
/// Retrieve colvar values and calculate their biasing forces
virtual int update();
- // TODO the following can be supplanted by a new call to init()
/// Load new configuration - force constant and/or centers only
- virtual void change_configuration(std::string const &conf);
+ virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
/// Calculate change in energy from using alternate configuration
- virtual cvm::real energy_difference(std::string const &conf);
-
- /// Read the bias configuration from a restart file
- virtual std::istream & read_restart(std::istream &is);
+ virtual cvm::real energy_difference(std::string const &conf) { return 0.0; }
- /// Write the bias configuration to a restart file
- virtual std::ostream & write_restart(std::ostream &os);
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ // virtual std::ostream & write_state_data(std::ostream &os);
+ // virtual std::istream & read_state_data(std::istream &os);
+ virtual std::ostream & write_state(std::ostream &os);
+ virtual std::istream & read_state(std::istream &is);
- /// Write a label to the trajectory file (comment line)
virtual std::ostream & write_traj_label(std::ostream &os);
-
- /// Output quantities such as the bias energy to the trajectory file
virtual std::ostream & write_traj(std::ostream &os);
/// \brief Constructor
@@ -42,26 +49,110 @@ public:
protected:
- /// \brief Potential function
- virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const = 0;
+ /// \brief Potential function for the i-th colvar
+ virtual cvm::real restraint_potential(size_t i) const = 0;
- /// \brief Force function
- virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const = 0;
+ /// \brief Force function for the i-th colvar
+ virtual colvarvalue const restraint_force(size_t i) const = 0;
+
+ /// \brief Derivative of the potential function with respect to the force constant
+ virtual cvm::real d_restraint_potential_dk(size_t i) const = 0;
+};
- ///\brief Unit scaling
- virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0;
+
+/// Definition and parsing of the restraint centers
+class colvarbias_restraint_centers
+ : public virtual colvarbias_restraint
+{
+public:
+
+ colvarbias_restraint_centers(char const *key);
+ virtual int init(std::string const &conf);
+ virtual int change_configuration(std::string const &conf);
+
+protected:
/// \brief Restraint centers
std::vector<colvarvalue> colvar_centers;
- /// \brief Restraint centers without wrapping or constraints applied
+ /// \brief Restraint centers outside the domain of the colvars (no wrapping or constraints applied)
std::vector<colvarvalue> colvar_centers_raw;
+};
+
+
+/// Definition and parsing of the force constant
+class colvarbias_restraint_k
+ : public virtual colvarbias_restraint
+{
+public:
+
+ colvarbias_restraint_k(char const *key);
+ virtual int init(std::string const &conf);
+ virtual int change_configuration(std::string const &conf);
+
+protected:
+ /// \brief Restraint force constant
+ cvm::real force_k;
+};
+
+
+/// Options to change the restraint configuration over time (shared between centers and k moving)
+class colvarbias_restraint_moving
+ : public virtual colvarparse {
+public:
+
+ colvarbias_restraint_moving(char const *key);
+ // Note: despite the diamond inheritance, most of this function gets only executed once
+ virtual int init(std::string const &conf);
+ virtual int update() { return COLVARS_OK; }
+ virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+
+protected:
/// \brief Moving target?
bool b_chg_centers;
+ /// \brief Changing force constant?
+ bool b_chg_force_k;
+
+ /// \brief Number of stages over which to perform the change
+ /// If zero, perform a continuous change
+ int target_nstages;
+
+ /// \brief Number of current stage of the perturbation
+ int stage;
+
+ /// \brief Lambda-schedule for custom varying force constant
+ std::vector<cvm::real> lambda_schedule;
+
+ /// \brief Number of steps required to reach the target force constant
+ /// or restraint centers
+ long target_nsteps;
+};
+
+
+/// Options to change the restraint centers over time
+class colvarbias_restraint_centers_moving
+ : public virtual colvarbias_restraint_centers,
+ public virtual colvarbias_restraint_moving
+{
+public:
+
+ colvarbias_restraint_centers_moving(char const *key);
+ virtual int init(std::string const &conf);
+ virtual int update();
+ virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_traj_label(std::ostream &os);
+ virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
/// \brief New restraint centers
std::vector<colvarvalue> target_centers;
@@ -78,11 +169,29 @@ protected:
/// \brief Accumulated work
cvm::real acc_work;
- /// \brief Restraint force constant
- cvm::real force_k;
+ /// Update the accumulated work
+ int update_acc_work();
+};
- /// \brief Changing force constant?
- bool b_chg_force_k;
+
+/// Options to change the restraint force constant over time
+class colvarbias_restraint_k_moving
+ : public virtual colvarbias_restraint_k,
+ public virtual colvarbias_restraint_moving
+{
+public:
+
+ colvarbias_restraint_k_moving(char const *key);
+ virtual int init(std::string const &conf);
+ virtual int update();
+ virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_traj_label(std::ostream &os);
+ virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
/// \brief Restraint force constant (target value)
cvm::real target_force_k;
@@ -90,9 +199,6 @@ protected:
/// \brief Restraint force constant (starting value)
cvm::real starting_force_k;
- /// \brief Lambda-schedule for custom varying force constant
- std::vector<cvm::real> lambda_schedule;
-
/// \brief Exponent for varying the force constant
cvm::real force_k_exp;
@@ -100,71 +206,91 @@ protected:
/// (in TI, would be the accumulating FE derivative)
cvm::real restraint_FE;
-
/// \brief Equilibration steps for restraint FE calculation through TI
cvm::real target_equil_steps;
+};
- /// \brief Number of stages over which to perform the change
- /// If zero, perform a continuous change
- int target_nstages;
- /// \brief Number of current stage of the perturbation
- int stage;
+/// \brief Harmonic bias restraint
+/// (implementation of \link colvarbias_restraint \endlink)
+class colvarbias_restraint_harmonic
+ : public colvarbias_restraint_centers_moving,
+ public colvarbias_restraint_k_moving
+{
+public:
+ colvarbias_restraint_harmonic(char const *key);
+ virtual int init(std::string const &conf);
+ virtual int update();
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_traj_label(std::ostream &os);
+ virtual std::ostream & write_traj(std::ostream &os);
+ virtual int change_configuration(std::string const &conf);
+ virtual cvm::real energy_difference(std::string const &conf);
- /// \brief Number of steps required to reach the target force constant
- /// or restraint centers
- long target_nsteps;
+protected:
+
+ virtual cvm::real restraint_potential(size_t i) const;
+ virtual colvarvalue const restraint_force(size_t i) const;
+ virtual cvm::real d_restraint_potential_dk(size_t i) const;
};
-/// \brief Harmonic bias restraint
+/// \brief Wall restraint
/// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_harmonic : public colvarbias_restraint {
-
+class colvarbias_restraint_harmonic_walls
+ : public colvarbias_restraint_k_moving
+{
public:
- colvarbias_restraint_harmonic(char const *key);
+ colvarbias_restraint_harmonic_walls(char const *key);
virtual int init(std::string const &conf);
- // no additional members, destructor not needed
+ virtual int update();
+ virtual void communicate_forces();
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_traj_label(std::ostream &os);
+ virtual std::ostream & write_traj(std::ostream &os);
protected:
- /// \brief Potential function
- virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const;
-
- /// \brief Force function
- virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const;
+ /// \brief Location of the lower walls
+ std::vector<colvarvalue> lower_walls;
- ///\brief Unit scaling
- virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
+ /// \brief Location of the upper walls
+ std::vector<colvarvalue> upper_walls;
+ 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;
+ virtual cvm::real d_restraint_potential_dk(size_t i) const;
};
/// \brief Linear bias restraint
/// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_linear : public colvarbias_restraint {
+class colvarbias_restraint_linear
+ : public colvarbias_restraint_centers_moving,
+ public colvarbias_restraint_k_moving
+{
public:
colvarbias_restraint_linear(char const *key);
virtual int init(std::string const &conf);
- // no additional members, destructor not needed
-
-protected:
-
- /// \brief Potential function
- virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const;
+ virtual int update();
+ virtual int change_configuration(std::string const &conf);
+ virtual cvm::real energy_difference(std::string const &conf);
- /// \brief Force function
- virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
- colvarvalue const &xcenter) const;
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_traj_label(std::ostream &os);
+ virtual std::ostream & write_traj(std::ostream &os);
- ///\brief Unit scaling
- virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
+protected:
+ virtual cvm::real restraint_potential(size_t i) const;
+ virtual colvarvalue const restraint_force(size_t i) const;
+ virtual cvm::real d_restraint_potential_dk(size_t i) const;
};
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 050e03f78ed520d7d4068a16734c28eb24ebffc4..d99da56386f05970d4d90f30d1cee27811390aa0 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvar.h"
@@ -154,15 +161,17 @@ void colvar::cvc::read_data()
void colvar::cvc::calc_force_invgrads()
{
- cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
- "for colvar components of type \""+function_type+"\".\n");
+ cvm::error("Error: calculation of inverse gradients is not implemented "
+ "for colvar components of type \""+function_type+"\".\n",
+ COLVARS_NOT_IMPLEMENTED);
}
void colvar::cvc::calc_Jacobian_derivative()
{
- cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
- "for colvar components of type \""+function_type+"\".\n");
+ cvm::error("Error: calculation of inverse gradients is not implemented "
+ "for colvar components of type \""+function_type+"\".\n",
+ COLVARS_NOT_IMPLEMENTED);
}
@@ -281,6 +290,33 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
}
+cvm::real colvar::cvc::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x1.dist2(x2);
+}
+
+
+colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x1.dist2_grad(x2);
+}
+
+
+colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x2.dist2_grad(x1);
+}
+
+
+void colvar::cvc::wrap(colvarvalue &x) const
+{
+ return;
+}
+
+
// Static members
std::vector<colvardeps::feature *> colvar::cvc::cvc_features;
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 0649895e452ca44f19a2cec99bfe2da9871aaf2f..a230cae8a9a04bf9b42e4b769f8e693062b8926d 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARCOMP_H
#define COLVARCOMP_H
@@ -105,8 +112,8 @@ public:
/// options from the provided configuration string
/// Returns reference to new group
cvm::atom_group *parse_group(std::string const &conf,
- char const *group_key,
- bool optional = false);
+ char const *group_key,
+ bool optional = false);
/// \brief Parse options pertaining to total force calculation
virtual int init_total_force_params(std::string const &conf);
@@ -130,7 +137,7 @@ public:
}
/// \brief Obtain data needed for the calculation for the backend
- void read_data();
+ virtual void read_data();
/// \brief Calculate the variable
virtual void calc_value() = 0;
@@ -151,17 +158,14 @@ public:
/// \brief Return the previously calculated value
- virtual colvarvalue const & value() const;
-
- // /// \brief Return const pointer to the previously calculated value
- // virtual const colvarvalue *p_value() const;
+ colvarvalue const & value() const;
/// \brief Return the previously calculated total force
- virtual colvarvalue const & total_force() const;
+ colvarvalue const & total_force() const;
/// \brief Return the previously calculated divergence of the
/// inverse atomic gradients
- virtual colvarvalue const & Jacobian_derivative() const;
+ colvarvalue const & Jacobian_derivative() const;
/// \brief Apply the collective variable force, by communicating the
/// atomic forces to the simulation program (\b Note: the \link ft
@@ -247,52 +251,24 @@ protected:
};
-
-
inline colvarvalue const & colvar::cvc::value() const
{
return x;
}
-// inline const colvarvalue * colvar::cvc::p_value() const
-// {
-// return &x;
-// }
inline colvarvalue const & colvar::cvc::total_force() const
{
return ft;
}
+
inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
{
return jd;
}
-inline cvm::real colvar::cvc::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x1.dist2(x2);
-}
-
-inline colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x1.dist2_grad(x2);
-}
-
-inline colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x2.dist2_grad(x1);
-}
-
-inline void colvar::cvc::wrap(colvarvalue &x) const
-{
- return;
-}
-
/// \brief Colvar component: distance between the centers of mass of
/// two groups (colvarvalue::type_scalar type, range [0:*))
@@ -312,7 +288,7 @@ protected:
public:
distance(std::string const &conf);
distance();
- virtual inline ~distance() {}
+ virtual ~distance() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -327,6 +303,7 @@ public:
};
+
// \brief Colvar component: distance vector between centers of mass
// of two groups (\link colvarvalue::type_3vector \endlink type,
// range (-*:*)x(-*:*)x(-*:*))
@@ -336,7 +313,7 @@ class colvar::distance_vec
public:
distance_vec(std::string const &conf);
distance_vec();
- virtual inline ~distance_vec() {}
+ virtual ~distance_vec() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -352,6 +329,7 @@ public:
};
+
/// \brief Colvar component: distance unit vector (direction) between
/// centers of mass of two groups (colvarvalue::type_unit3vector type,
/// range [-1:1]x[-1:1]x[-1:1])
@@ -361,19 +339,14 @@ class colvar::distance_dir
public:
distance_dir(std::string const &conf);
distance_dir();
- virtual inline ~distance_dir() {}
+ virtual ~distance_dir() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
- virtual cvm::real dist2(colvarvalue const &x1,
- colvarvalue const &x2) const;
- virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const;
- virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const;
};
+
/// \brief Colvar component: projection of the distance vector along
/// an axis(colvarvalue::type_scalar type, range (-*:*))
class colvar::distance_z
@@ -399,7 +372,7 @@ protected:
public:
distance_z(std::string const &conf);
distance_z();
- virtual inline ~distance_z() {}
+ virtual ~distance_z() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -416,6 +389,7 @@ public:
};
+
/// \brief Colvar component: projection of the distance vector on a
/// plane (colvarvalue::type_scalar type, range [0:*))
class colvar::distance_xy
@@ -429,7 +403,7 @@ protected:
public:
distance_xy(std::string const &conf);
distance_xy();
- virtual inline ~distance_xy() {}
+ virtual ~distance_xy() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -444,6 +418,7 @@ public:
};
+
/// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
/// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
class colvar::distance_inv
@@ -455,7 +430,7 @@ protected:
public:
distance_inv(std::string const &conf);
distance_inv();
- virtual inline ~distance_inv() {}
+ virtual ~distance_inv() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -468,6 +443,7 @@ public:
};
+
/// \brief Colvar component: N1xN2 vector of pairwise distances
/// (colvarvalue::type_vector type, range (0:*) for each component)
class colvar::distance_pairs
@@ -483,16 +459,10 @@ protected:
public:
distance_pairs(std::string const &conf);
distance_pairs();
- virtual inline ~distance_pairs() {}
+ virtual ~distance_pairs() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
- virtual cvm::real dist2(colvarvalue const &x1,
- colvarvalue const &x2) const;
- virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const;
- virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -509,7 +479,7 @@ public:
/// Constructor
gyration(std::string const &conf);
gyration();
- virtual inline ~gyration() {}
+ virtual ~gyration() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -524,6 +494,7 @@ public:
};
+
/// \brief Colvar component: moment of inertia of an atom group
/// (colvarvalue::type_scalar type, range [0:*))
class colvar::inertia
@@ -533,7 +504,7 @@ public:
/// Constructor
inertia(std::string const &conf);
inertia();
- virtual inline ~inertia() {}
+ virtual ~inertia() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -546,6 +517,7 @@ public:
};
+
/// \brief Colvar component: moment of inertia of an atom group
/// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
class colvar::inertia_z
@@ -558,7 +530,7 @@ public:
/// Constructor
inertia_z(std::string const &conf);
inertia_z();
- virtual inline ~inertia_z() {}
+ virtual ~inertia_z() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -571,6 +543,7 @@ public:
};
+
/// \brief Colvar component: projection of 3N coordinates onto an
/// eigenvector(colvarvalue::type_scalar type, range (-*:*))
class colvar::eigenvector
@@ -597,7 +570,7 @@ public:
/// Constructor
eigenvector(std::string const &conf);
- virtual inline ~eigenvector() {}
+ virtual ~eigenvector() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -645,7 +618,7 @@ public:
/// \brief Initialize the three groups after three atoms
angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
angle();
- virtual inline ~angle() {}
+ virtual ~angle() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -659,6 +632,8 @@ public:
colvarvalue const &x2) const;
};
+
+
/// \brief Colvar component: angle between the dipole of a molecule and an axis
/// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
class colvar::dipole_angle
@@ -691,7 +666,7 @@ public:
/// \brief Initialize the three groups after three atoms
dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
dipole_angle();
- virtual inline ~dipole_angle() {}
+ virtual ~dipole_angle() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force (colvarvalue const &force);
@@ -703,6 +678,8 @@ public:
colvarvalue const &x2) const;
};
+
+
/// \brief Colvar component: dihedral between the centers of mass of
/// four groups (colvarvalue::type_scalar type, range [-PI:PI])
class colvar::dihedral
@@ -732,7 +709,7 @@ public:
/// \brief Initialize the four groups after four atoms
dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
dihedral();
- virtual inline ~dihedral() {}
+ virtual ~dihedral() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -753,6 +730,7 @@ public:
};
+
/// \brief Colvar component: coordination number between two groups
/// (colvarvalue::type_scalar type, range [0:N1*N2])
class colvar::coordnum
@@ -781,7 +759,7 @@ public:
/// Constructor
coordnum(std::string const &conf);
coordnum();
- virtual inline ~coordnum() {}
+ virtual ~coordnum() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -812,6 +790,8 @@ public:
colvarvalue const &x2) const;
};
+
+
/// \brief Colvar component: self-coordination number within a group
/// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
class colvar::selfcoordnum
@@ -830,7 +810,7 @@ public:
/// Constructor
selfcoordnum(std::string const &conf);
selfcoordnum();
- virtual inline ~selfcoordnum() {}
+ virtual ~selfcoordnum() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -852,6 +832,7 @@ public:
};
+
/// \brief Colvar component: coordination number between two groups
/// (colvarvalue::type_scalar type, range [0:N1*N2])
class colvar::groupcoordnum
@@ -873,7 +854,7 @@ public:
/// Constructor
groupcoordnum(std::string const &conf);
groupcoordnum();
- virtual inline ~groupcoordnum() {}
+ virtual ~groupcoordnum() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -896,6 +877,7 @@ public:
static cvm::real switching_function(cvm::rvector const &r0_vec,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
+ */
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
@@ -903,10 +885,10 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
- */
};
+
/// \brief Colvar component: hydrogen bond, defined as the product of
/// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
/// (colvarvalue::type_scalar type, range [0:1])
@@ -941,6 +923,7 @@ public:
};
+
/// \brief Colvar component: alpha helix content of a contiguous
/// segment of 5 or more residues, implemented as a sum of phi/psi
/// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
@@ -969,7 +952,7 @@ public:
// alpha_dihedrals (std::string const &conf);
// alpha_dihedrals();
-// virtual inline ~alpha_dihedrals() {}
+// virtual ~alpha_dihedrals() {}
// virtual void calc_value();
// virtual void calc_gradients();
// virtual void apply_force (colvarvalue const &force);
@@ -982,6 +965,7 @@ public:
// };
+
/// \brief Colvar component: alpha helix content of a contiguous
/// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
/// angles and hydrogen bonds (colvarvalue::type_scalar type, range
@@ -1022,6 +1006,8 @@ public:
colvarvalue const &x2) const;
};
+
+
/// \brief Colvar component: dihedPC
/// Projection of the config onto a dihedral principal component
/// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
@@ -1050,6 +1036,8 @@ public:
colvarvalue const &x2) const;
};
+
+
/// \brief Colvar component: orientation in space of an atom group,
/// with respect to a set of reference coordinates
/// (colvarvalue::type_quaternion type, range
@@ -1078,7 +1066,7 @@ public:
orientation(std::string const &conf);
orientation();
- virtual inline ~orientation() {}
+ virtual ~orientation() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1091,6 +1079,7 @@ public:
};
+
/// \brief Colvar component: angle of rotation with respect to a set
/// of reference coordinates (colvarvalue::type_scalar type, range
/// [0:PI))
@@ -1101,7 +1090,7 @@ public:
orientation_angle(std::string const &conf);
orientation_angle();
- virtual inline ~orientation_angle() {}
+ virtual ~orientation_angle() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1114,6 +1103,7 @@ public:
};
+
/// \brief Colvar component: cosine of the angle of rotation with respect to a set
/// of reference coordinates (colvarvalue::type_scalar type, range
/// [-1:1])
@@ -1124,7 +1114,7 @@ public:
orientation_proj(std::string const &conf);
orientation_proj();
- virtual inline ~orientation_proj() {}
+ virtual ~orientation_proj() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1137,6 +1127,7 @@ public:
};
+
/// \brief Colvar component: projection of the orientation vector onto
/// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
class colvar::tilt
@@ -1150,7 +1141,7 @@ public:
tilt(std::string const &conf);
tilt();
- virtual inline ~tilt() {}
+ virtual ~tilt() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1177,7 +1168,7 @@ public:
spin_angle(std::string const &conf);
spin_angle();
- virtual inline ~spin_angle() {}
+ virtual ~spin_angle() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1215,7 +1206,7 @@ public:
/// Constructor
rmsd(std::string const &conf);
- virtual inline ~rmsd() {}
+ virtual ~rmsd() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void calc_force_invgrads();
@@ -1230,6 +1221,7 @@ public:
};
+
// \brief Colvar component: flat vector of Cartesian coordinates
// Mostly useful to compute scripted colvar values
class colvar::cartesian
@@ -1243,7 +1235,7 @@ protected:
public:
cartesian(std::string const &conf);
cartesian();
- virtual inline ~cartesian() {}
+ virtual ~cartesian() {}
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
@@ -1260,255 +1252,26 @@ public:
#define simple_scalar_dist_functions(TYPE) \
\
- inline cvm::real colvar::TYPE::dist2(colvarvalue const &x1, \
- colvarvalue const &x2) const \
+ \
+ cvm::real colvar::TYPE::dist2(colvarvalue const &x1, \
+ colvarvalue const &x2) const \
{ \
return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
} \
\
- inline colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1, \
- colvarvalue const &x2) const \
+ \
+ colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1, \
+ colvarvalue const &x2) const \
{ \
return 2.0 * (x1.real_value - x2.real_value); \
} \
\
- inline colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1, \
- colvarvalue const &x2) const \
+ \
+ colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1, \
+ colvarvalue const &x2) const \
{ \
return this->dist2_lgrad(x2, x1); \
} \
\
-simple_scalar_dist_functions(distance)
-// NOTE: distance_z has explicit functions, see below
-simple_scalar_dist_functions(distance_xy)
-simple_scalar_dist_functions(distance_inv)
-simple_scalar_dist_functions(angle)
-simple_scalar_dist_functions(dipole_angle)
-simple_scalar_dist_functions(coordnum)
-simple_scalar_dist_functions(selfcoordnum)
-simple_scalar_dist_functions(h_bond)
-simple_scalar_dist_functions(gyration)
-simple_scalar_dist_functions(inertia)
-simple_scalar_dist_functions(inertia_z)
-simple_scalar_dist_functions(rmsd)
-simple_scalar_dist_functions(orientation_angle)
-simple_scalar_dist_functions(orientation_proj)
-simple_scalar_dist_functions(tilt)
-simple_scalar_dist_functions(eigenvector)
-// simple_scalar_dist_functions (alpha_dihedrals)
-simple_scalar_dist_functions(alpha_angles)
-simple_scalar_dist_functions(dihedPC)
-
-
-// metrics functions for cvc implementations with a periodicity
-
-inline cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return diff * diff;
-}
-
-inline colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return 2.0 * diff;
-}
-
-inline colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return (-2.0) * diff;
-}
-
-inline void colvar::dihedral::wrap(colvarvalue &x) const
-{
-if ((x.real_value - wrap_center) >= 180.0) {
-x.real_value -= 360.0;
-return;
-}
-
-if ((x.real_value - wrap_center) < -180.0) {
-x.real_value += 360.0;
-return;
-}
-
-return;
-}
-
-inline cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return diff * diff;
-}
-
-inline colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return 2.0 * diff;
-}
-
-inline colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return (-2.0) * diff;
-}
-
-inline void colvar::spin_angle::wrap(colvarvalue &x) const
-{
-if ((x.real_value - wrap_center) >= 180.0) {
-x.real_value -= 360.0;
-return;
-}
-
-if ((x.real_value - wrap_center) < -180.0) {
-x.real_value += 360.0;
-return;
-}
-
-return;
-}
-
-
-// Projected distance
-// Differences should always be wrapped around 0 (ignoring wrap_center)
-inline cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::real diff = x1.real_value - x2.real_value;
- if (period != 0.0) {
- cvm::real shift = std::floor(diff/period + 0.5);
- diff -= shift * period;
- }
- return diff * diff;
-}
-
-inline colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::real diff = x1.real_value - x2.real_value;
- if (period != 0.0) {
- cvm::real shift = std::floor(diff/period + 0.5);
- diff -= shift * period;
- }
- return 2.0 * diff;
-}
-
-inline colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::real diff = x1.real_value - x2.real_value;
- if (period != 0.0) {
- cvm::real shift = std::floor(diff/period + 0.5);
- diff -= shift * period;
- }
- return (-2.0) * diff;
-}
-
-inline void colvar::distance_z::wrap(colvarvalue &x) const
-{
- if (! this->b_periodic) {
- // don't wrap if the period has not been set
- return;
- }
-
- cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
- x.real_value -= shift * period;
- return;
-}
-
-
-// distance between three dimensional vectors
-//
-// TODO apply PBC to distance_vec
-// Note: differences should be centered around (0, 0, 0)!
-
-inline cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
-}
-
-inline colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
-}
-
-inline colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
-}
-
-inline cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return (x1.rvector_value - x2.rvector_value).norm2();
-}
-
-inline colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
-}
-
-inline colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
-}
-
-inline cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return (x1.vector1d_value - x2.vector1d_value).norm2();
-}
-
-inline colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return colvarvalue((x1.vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
-}
-
-inline colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return colvarvalue((x2.vector1d_value - x1.vector1d_value), colvarvalue::type_vector);
-}
-
-
-// distance between quaternions
-
-inline cvm::real colvar::orientation::dist2(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x1.quaternion_value.dist2(x2);
-}
-
-inline colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x1.quaternion_value.dist2_grad(x2);
-}
-
-inline colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return x2.quaternion_value.dist2_grad(x1);
-}
-
-
#endif
diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp
index 6b05fb16eb50522756f32b2cc531005c7af9ef45..4c3390a8bd1de9bf0fb37288d114d6cb8177daca 100644
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include "colvarmodule.h"
#include "colvar.h"
#include "colvarcomp.h"
@@ -7,6 +14,7 @@
#include <cmath>
+
colvar::angle::angle(std::string const &conf)
: cvc(conf)
{
@@ -85,6 +93,7 @@ void colvar::angle::calc_gradients()
group3->set_weighted_gradient(dxdr3);
}
+
void colvar::angle::calc_force_invgrads()
{
// This uses a force measurement on groups 1 and 3 only
@@ -107,6 +116,7 @@ void colvar::angle::calc_force_invgrads()
return;
}
+
void colvar::angle::calc_Jacobian_derivative()
{
// det(J) = (2 pi) r^2 * sin(theta)
@@ -129,6 +139,8 @@ void colvar::angle::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(angle)
+
colvar::dipole_angle::dipole_angle(std::string const &conf)
@@ -235,6 +247,8 @@ void colvar::dipole_angle::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(dipole_angle)
+
colvar::dihedral::dihedral(std::string const &conf)
@@ -453,3 +467,46 @@ void colvar::dihedral::apply_force(colvarvalue const &force)
}
+// metrics functions for cvc implementations with a periodicity
+
+cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return diff * diff;
+}
+
+
+colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return 2.0 * diff;
+}
+
+
+colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return (-2.0) * diff;
+}
+
+
+void colvar::dihedral::wrap(colvarvalue &x) const
+{
+ if ((x.real_value - wrap_center) >= 180.0) {
+ x.real_value -= 360.0;
+ return;
+ }
+
+ if ((x.real_value - wrap_center) < -180.0) {
+ x.real_value += 360.0;
+ return;
+ }
+
+ return;
+}
diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp
index d5e87e6b1e3f0388d0626d3370a342c5545e7929..987a16a816b0f3869025a1e612ab978afe954328 100644
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <cmath>
#include "colvarmodule.h"
@@ -13,10 +20,10 @@
template<bool calculate_gradients>
cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
- int const &en,
- int const &ed,
- cvm::atom &A1,
- cvm::atom &A2)
+ int const &en,
+ int const &ed,
+ cvm::atom &A1,
+ cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::real const l2 = diff.norm2()/(r0*r0);
@@ -42,10 +49,10 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
template<bool calculate_gradients>
cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
- int const &en,
- int const &ed,
- cvm::atom &A1,
- cvm::atom &A2)
+ int const &en,
+ int const &ed,
+ cvm::atom &A1,
+ cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
@@ -190,6 +197,7 @@ void colvar::coordnum::calc_gradients()
}
}
+
void colvar::coordnum::apply_force(colvarvalue const &force)
{
if (!group1->noforce)
@@ -200,6 +208,9 @@ void colvar::coordnum::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(coordnum)
+
+
// h_bond member functions
@@ -252,6 +263,7 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor,
atom_groups[0]->add_atom(donor);
}
+
colvar::h_bond::h_bond()
: cvc()
{
@@ -284,6 +296,8 @@ void colvar::h_bond::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(h_bond)
+
colvar::selfcoordnum::selfcoordnum(std::string const &conf)
@@ -339,6 +353,9 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(selfcoordnum)
+
+
// groupcoordnum member functions
colvar::groupcoordnum::groupcoordnum(std::string const &conf)
: distance(conf), b_anisotropic(false)
@@ -448,7 +465,6 @@ cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
#endif
-
void colvar::groupcoordnum::calc_value()
{
@@ -460,7 +476,6 @@ void colvar::groupcoordnum::calc_value()
x.real_value = coordnum::switching_function<false>(r0, en, ed,
group1_com_atom, group2_com_atom);
-
}
@@ -486,3 +501,6 @@ void colvar::groupcoordnum::apply_force(colvarvalue const &force)
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
}
+
+
+simple_scalar_dist_functions(groupcoordnum)
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 9ef96af2d3ab5bb4a2d3a7edbf55c3fa9e24c360..7481dd360ba170e327a31977d21c33d47f4507a1 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <cmath>
#include "colvarmodule.h"
@@ -91,6 +98,9 @@ void colvar::distance::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(distance)
+
+
colvar::distance_vec::distance_vec(std::string const &conf)
: distance(conf)
@@ -138,6 +148,27 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
}
+cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
+}
+
+
+colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
+}
+
+
+colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
+}
+
+
colvar::distance_z::distance_z(std::string const &conf)
: cvc(conf)
@@ -191,6 +222,7 @@ colvar::distance_z::distance_z(std::string const &conf)
}
+
colvar::distance_z::distance_z()
{
function_type = "distance_z";
@@ -200,6 +232,7 @@ colvar::distance_z::distance_z()
x.type(colvarvalue::type_scalar);
}
+
void colvar::distance_z::calc_value()
{
if (fixed_axis) {
@@ -227,6 +260,7 @@ void colvar::distance_z::calc_value()
this->wrap(x);
}
+
void colvar::distance_z::calc_gradients()
{
main->set_weighted_gradient( axis );
@@ -248,6 +282,7 @@ void colvar::distance_z::calc_gradients()
}
}
+
void colvar::distance_z::calc_force_invgrads()
{
main->read_total_forces();
@@ -260,11 +295,13 @@ void colvar::distance_z::calc_force_invgrads()
}
}
+
void colvar::distance_z::calc_Jacobian_derivative()
{
jd.real_value = 0.0;
}
+
void colvar::distance_z::apply_force(colvarvalue const &force)
{
if (!ref1->noforce)
@@ -278,6 +315,56 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
}
+// Differences should always be wrapped around 0 (ignoring wrap_center)
+cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ if (b_periodic) {
+ cvm::real shift = std::floor(diff/period + 0.5);
+ diff -= shift * period;
+ }
+ return diff * diff;
+}
+
+
+colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ if (b_periodic) {
+ cvm::real shift = std::floor(diff/period + 0.5);
+ diff -= shift * period;
+ }
+ return 2.0 * diff;
+}
+
+
+colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ if (b_periodic) {
+ cvm::real shift = std::floor(diff/period + 0.5);
+ diff -= shift * period;
+ }
+ return (-2.0) * diff;
+}
+
+
+void colvar::distance_z::wrap(colvarvalue &x) const
+{
+ if (!b_periodic) {
+ // don't wrap if the period has not been set
+ return;
+ }
+
+ cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
+ x.real_value -= shift * period;
+ return;
+}
+
+
colvar::distance_xy::distance_xy(std::string const &conf)
: distance_z(conf)
@@ -289,6 +376,7 @@ colvar::distance_xy::distance_xy(std::string const &conf)
x.type(colvarvalue::type_scalar);
}
+
colvar::distance_xy::distance_xy()
: distance_z()
{
@@ -299,6 +387,7 @@ colvar::distance_xy::distance_xy()
x.type(colvarvalue::type_scalar);
}
+
void colvar::distance_xy::calc_value()
{
if (b_no_PBC) {
@@ -321,6 +410,7 @@ void colvar::distance_xy::calc_value()
x.real_value = dist_v_ortho.norm();
}
+
void colvar::distance_xy::calc_gradients()
{
// Intermediate quantity (r_P3 / r_12 where P is the projection
@@ -348,6 +438,7 @@ void colvar::distance_xy::calc_gradients()
}
}
+
void colvar::distance_xy::calc_force_invgrads()
{
main->read_total_forces();
@@ -360,11 +451,13 @@ void colvar::distance_xy::calc_force_invgrads()
}
}
+
void colvar::distance_xy::calc_Jacobian_derivative()
{
jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
}
+
void colvar::distance_xy::apply_force(colvarvalue const &force)
{
if (!ref1->noforce)
@@ -378,6 +471,9 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(distance_xy)
+
+
colvar::distance_dir::distance_dir(std::string const &conf)
: distance(conf)
@@ -403,7 +499,7 @@ void colvar::distance_dir::calc_value()
dist_v = group2->center_of_mass() - group1->center_of_mass();
} else {
dist_v = cvm::position_distance(group1->center_of_mass(),
- group2->center_of_mass());
+ group2->center_of_mass());
}
x.rvector_value = dist_v.unit();
}
@@ -460,6 +556,7 @@ colvar::distance_inv::distance_inv(std::string const &conf)
x.type(colvarvalue::type_scalar);
}
+
colvar::distance_inv::distance_inv()
{
function_type = "distance_inv";
@@ -467,6 +564,7 @@ colvar::distance_inv::distance_inv()
x.type(colvarvalue::type_scalar);
}
+
void colvar::distance_inv::calc_value()
{
x.real_value = 0.0;
@@ -504,6 +602,7 @@ void colvar::distance_inv::calc_value()
x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));
}
+
void colvar::distance_inv::calc_gradients()
{
cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size());
@@ -515,6 +614,7 @@ void colvar::distance_inv::calc_gradients()
}
}
+
void colvar::distance_inv::apply_force(colvarvalue const &force)
{
if (!group1->noforce)
@@ -525,6 +625,9 @@ void colvar::distance_inv::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(distance_inv)
+
+
colvar::distance_pairs::distance_pairs(std::string const &conf)
: cvc(conf)
@@ -579,11 +682,13 @@ void colvar::distance_pairs::calc_value()
}
}
+
void colvar::distance_pairs::calc_gradients()
{
// will be calculated on the fly in apply_force()
}
+
void colvar::distance_pairs::apply_force(colvarvalue const &force)
{
if (b_no_PBC) {
@@ -608,6 +713,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
}
+
colvar::gyration::gyration(std::string const &conf)
: cvc(conf)
{
@@ -682,6 +788,9 @@ void colvar::gyration::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(gyration)
+
+
colvar::inertia::inertia(std::string const &conf)
: gyration(conf)
@@ -722,6 +831,10 @@ void colvar::inertia::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(inertia_z)
+
+
+
colvar::inertia_z::inertia_z(std::string const &conf)
: inertia(conf)
{
@@ -772,6 +885,10 @@ void colvar::inertia_z::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(inertia)
+
+
+
colvar::rmsd::rmsd(std::string const &conf)
: cvc(conf)
@@ -971,6 +1088,8 @@ void colvar::rmsd::calc_Jacobian_derivative()
}
+simple_scalar_dist_functions(rmsd)
+
colvar::eigenvector::eigenvector(std::string const &conf)
@@ -1255,6 +1374,10 @@ void colvar::eigenvector::calc_Jacobian_derivative()
}
+simple_scalar_dist_functions(eigenvector)
+
+
+
colvar::cartesian::cartesian(std::string const &conf)
: cvc(conf)
{
diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp
index e5db3a4b508deb626c56fb205f7937aacd1301da..393c7dcf9aa2a9533ad553251e38082961a2d544 100644
--- a/lib/colvars/colvarcomp_protein.cpp
+++ b/lib/colvars/colvarcomp_protein.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <cmath>
#include "colvarmodule.h"
@@ -119,6 +126,7 @@ colvar::alpha_angles::alpha_angles()
x.type(colvarvalue::type_scalar);
}
+
colvar::alpha_angles::~alpha_angles()
{
while (theta.size() != 0) {
@@ -131,6 +139,7 @@ colvar::alpha_angles::~alpha_angles()
}
}
+
void colvar::alpha_angles::calc_value()
{
x.real_value = 0.0;
@@ -222,6 +231,10 @@ void colvar::alpha_angles::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(alpha_angles)
+
+
+
//////////////////////////////////////////////////////////////////////
// dihedral principal component
//////////////////////////////////////////////////////////////////////
@@ -337,14 +350,22 @@ colvar::dihedPC::dihedPC(std::string const &conf)
for (size_t i = 0; i < residues.size()-1; i++) {
// Psi
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "N", sid),
- cvm::atom(r[i ], "CA", sid),
- cvm::atom(r[i ], "C", sid),
- cvm::atom(r[i+1], "N", sid)));
+ cvm::atom(r[i ], "CA", sid),
+ cvm::atom(r[i ], "C", sid),
+ cvm::atom(r[i+1], "N", sid)));
+ atom_groups.push_back(theta.back()->atom_groups[0]);
+ atom_groups.push_back(theta.back()->atom_groups[1]);
+ atom_groups.push_back(theta.back()->atom_groups[2]);
+ atom_groups.push_back(theta.back()->atom_groups[3]);
// Phi (next res)
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "C", sid),
- cvm::atom(r[i+1], "N", sid),
- cvm::atom(r[i+1], "CA", sid),
- cvm::atom(r[i+1], "C", sid)));
+ cvm::atom(r[i+1], "N", sid),
+ cvm::atom(r[i+1], "CA", sid),
+ cvm::atom(r[i+1], "C", sid)));
+ atom_groups.push_back(theta.back()->atom_groups[0]);
+ atom_groups.push_back(theta.back()->atom_groups[1]);
+ atom_groups.push_back(theta.back()->atom_groups[2]);
+ atom_groups.push_back(theta.back()->atom_groups[3]);
}
if (cvm::debug())
@@ -400,3 +421,6 @@ void colvar::dihedPC::apply_force(colvarvalue const &force)
coeffs[2*i+1] * dsindt) * force);
}
}
+
+
+simple_scalar_dist_functions(dihedPC)
diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp
index a3cbbb2908ce0dd143bba26eef90f96d3c9d7906..936e770169d060895d02e14538d189241dbab6d5 100644
--- a/lib/colvars/colvarcomp_rotations.cpp
+++ b/lib/colvars/colvarcomp_rotations.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <cmath>
#include "colvarmodule.h"
@@ -123,6 +130,27 @@ void colvar::orientation::apply_force(colvarvalue const &force)
}
+cvm::real colvar::orientation::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x1.quaternion_value.dist2(x2);
+}
+
+
+colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x1.quaternion_value.dist2_grad(x2);
+}
+
+
+colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ return x2.quaternion_value.dist2_grad(x1);
+}
+
+
colvar::orientation_angle::orientation_angle(std::string const &conf)
: orientation(conf)
@@ -176,6 +204,9 @@ void colvar::orientation_angle::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(orientation_angle)
+
+
colvar::orientation_proj::orientation_proj(std::string const &conf)
: orientation(conf)
@@ -220,6 +251,9 @@ void colvar::orientation_proj::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(orientation_proj)
+
+
colvar::tilt::tilt(std::string const &conf)
: orientation(conf)
@@ -278,6 +312,9 @@ void colvar::tilt::apply_force(colvarvalue const &force)
}
+simple_scalar_dist_functions(tilt)
+
+
colvar::spin_angle::spin_angle(std::string const &conf)
: orientation(conf)
@@ -339,3 +376,46 @@ void colvar::spin_angle::apply_force(colvarvalue const &force)
atoms->apply_colvar_force(fw);
}
}
+
+
+cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return diff * diff;
+}
+
+
+colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return 2.0 * diff;
+}
+
+
+colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
+ colvarvalue const &x2) const
+{
+ cvm::real diff = x1.real_value - x2.real_value;
+ diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+ return (-2.0) * diff;
+}
+
+
+void colvar::spin_angle::wrap(colvarvalue &x) const
+{
+ if ((x.real_value - wrap_center) >= 180.0) {
+ x.real_value -= 360.0;
+ return;
+ }
+
+ if ((x.real_value - wrap_center) < -180.0) {
+ x.real_value += 360.0;
+ return;
+ }
+
+ return;
+}
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index a44f86c2925661d410f4cc6a321ed1ccbdcf4f7b..e3ccd4ecf5815fcabc0bf7b337a06bf32340c1e7 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -1,3 +1,13 @@
+// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
#include "colvardeps.h"
@@ -219,6 +229,9 @@ void colvardeps::init_cvb_requires() {
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++) {
@@ -229,6 +242,9 @@ void colvardeps::init_cvb_requires() {
// 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;
}
@@ -321,6 +337,7 @@ void colvardeps::init_cv_requires() {
// The features below are set programmatically
f_description(f_cv_scripted, "scripted");
f_description(f_cv_periodic, "periodic");
+ 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");
@@ -407,6 +424,11 @@ void colvardeps::init_cvc_requires() {
// Each cvc specifies what other features are available
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;
+
+ // Features That are implemented only for certain simulation engine configurations
feature_states[f_cvc_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
feature_states[f_cvc_scalable]->available = feature_states[f_cvc_scalable_com]->available;
}
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index 31dfd0e3023e027310156a463730ad9b378f54c0..4ef27ded85a5b2a1bbcaaad54933a9294e11c51a 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARDEPS_H
#define COLVARDEPS_H
@@ -157,6 +164,7 @@ public:
f_cvb_apply_force, // will apply forces
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_ntot
};
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index ca2d935e1caf8c6578a6238b4f2c42b4722c1645..3b25acd2ef92141af45b0b016f305ed64534f7eb 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// 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"
@@ -73,6 +80,21 @@ cvm::real colvar_grid_scalar::minimum_value() const
return min;
}
+cvm::real colvar_grid_scalar::minimum_pos_value() const
+{
+ cvm::real minpos = data[0];
+ size_t i;
+ for (i = 0; i < nt; i++) {
+ if(data[i] > 0) {
+ minpos = data[i];
+ break;
+ }
+ }
+ for (i = 0; i < nt; i++) {
+ if (data[i] > 0 && data[i] < minpos) minpos = data[i];
+ }
+ return minpos;
+}
cvm::real colvar_grid_scalar::integral() const
{
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 135e405ee7b71e7b580e988ed9c65c53c4d3512a..d4b9295c6e1901c0be09beee3bce7fcffe385d93 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARGRID_H
#define COLVARGRID_H
@@ -378,6 +385,13 @@ public:
return value_to_bin_scalar(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
}
+ /// \brief Report the bin corresponding to the current value of variable i
+ /// and assign first or last bin if out of boundaries
+ inline int current_bin_scalar_bound(int const i) const
+ {
+ return value_to_bin_scalar_bound(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
+ }
+
/// \brief Report the bin corresponding to the current value of item iv in variable i
inline int current_bin_scalar(int const i, int const iv) const
{
@@ -393,6 +407,16 @@ public:
return (int) std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
}
+ /// \brief Use the lower boundary and the width to report which bin
+ /// the provided value is in and assign first or last bin if out of boundaries
+ inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
+ {
+ int bin_index = std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
+ if (bin_index < 0) bin_index=0;
+ if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
+ return (int) bin_index;
+ }
+
/// \brief Same as the standard version, but uses another grid definition
inline int value_to_bin_scalar(colvarvalue const &value,
colvarvalue const &new_offset,
@@ -514,6 +538,13 @@ public:
data[i] *= a;
}
+ /// \brief Assign all zero elements a scalar constant (fast loop)
+ inline void remove_zeros(cvm::real const &a)
+ {
+ for (size_t i = 0; i < nt; i++)
+ if(data[i]==0) data[i] = a;
+ }
+
/// \brief Get the bin indices corresponding to the provided values of
/// the colvars
@@ -537,6 +568,17 @@ public:
return index;
}
+ /// \brief Get the bin indices corresponding to the provided values of
+ /// the colvars and assign first or last bin if out of boundaries
+ inline std::vector<int> const get_colvars_index_bound() const
+ {
+ std::vector<int> index = new_index();
+ for (size_t i = 0; i < nd; i++) {
+ index[i] = current_bin_scalar_bound(i);
+ }
+ return index;
+ }
+
/// \brief Get the minimal distance (in number of bins) from the
/// boundaries; a negative number is returned if the given point is
/// off-grid
@@ -1169,42 +1211,46 @@ public:
inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
int n = 0)
{
- cvm::real A0, A1;
- std::vector<int> ix;
-
- // factor for mesh width, 2.0 for central finite difference
- // but only 1.0 on edges for non-PBC coordinates
- cvm::real factor;
+ int A0, A1, A2;
+ std::vector<int> ix = ix0;
if (periodic[n]) {
- factor = 2.;
- ix = ix0;
ix[n]--; wrap(ix);
A0 = data[address(ix)];
ix = ix0;
ix[n]++; wrap(ix);
A1 = data[address(ix)];
- } else {
- factor = 0.;
- ix = ix0;
- if (ix[n] > 0) { // not left edge
- ix[n]--;
- factor += 1.;
+ if (A0 * A1 == 0) {
+ return 0.; // can't handle empty bins
+ } else {
+ return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+ / (widths[n] * 2.);
}
+ } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
+ ix[n]--;
A0 = data[address(ix)];
ix = ix0;
- if (ix[n]+1 < nx[n]) { // not right edge
- ix[n]++;
- factor += 1.;
- }
+ ix[n]++;
A1 = data[address(ix)];
- }
- if (A0 == 0 || A1 == 0) {
- // can't handle empty bins
- return 0.;
+ if (A0 * A1 == 0) {
+ return 0.; // can't handle empty bins
+ } else {
+ return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+ / (widths[n] * 2.);
+ }
} else {
- return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
- / (widths[n] * factor);
+ // edge: use 2nd order derivative
+ int increment = (ix[n] == 0 ? 1 : -1);
+ // move right from left edge, or the other way around
+ A0 = data[address(ix)];
+ ix[n] += increment; A1 = data[address(ix)];
+ ix[n] += increment; A2 = data[address(ix)];
+ if (A0 * A1 * A2 == 0) {
+ return 0.; // can't handle empty bins
+ } else {
+ return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1)
+ - 0.5 * std::log((cvm::real)A2)) * increment / widths[n];
+ }
}
}
};
@@ -1322,6 +1368,9 @@ public:
/// \brief Return the lowest value
cvm::real minimum_value() const;
+ /// \brief Return the lowest positive value
+ cvm::real minimum_pos_value() const;
+
/// \brief Calculates the integral of the map (uses widths if they are defined)
cvm::real integral() const;
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 8e110275b21501208c725c270175ff46e1f11813..b9a435152b2cc56b50987b08adb7cf2bd75e4362 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <sstream>
#include <string.h>
@@ -296,6 +303,9 @@ int colvarmodule::parse_biases(std::string const &conf)
/// initialize harmonic restraints
parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);
+ /// initialize harmonic walls restraints
+ parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls", n_rest_biases);
+
/// initialize histograms
parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
@@ -562,7 +572,6 @@ int colvarmodule::calc_colvars()
colvars_smp_items.reserve(colvars.size());
// set up a vector containing all components
- size_t num_colvar_items = 0;
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
@@ -576,8 +585,6 @@ int colvarmodule::calc_colvars()
colvars_smp.push_back(*cvi);
colvars_smp_items.push_back(icvc);
}
-
- num_colvar_items += num_items;
}
cvm::decrease_depth();
@@ -641,7 +648,7 @@ int colvarmodule::calc_biases()
for (bi = biases.begin(); bi != biases.end(); bi++) {
error_code |= (*bi)->update();
if (cvm::get_error()) {
- return COLVARS_ERROR;
+ return error_code;
}
}
cvm::decrease_depth();
@@ -1007,7 +1014,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
- if (!((*bi)->read_restart(is))) {
+ if (!((*bi)->read_state(is))) {
cvm::error("Error: in reading restart configuration for bias \""+
(*bi)->name+"\".\n",
INPUT_ERROR);
@@ -1070,15 +1077,15 @@ continue the previous simulation.\n\n");
cvm::log(cvm::line_marker);
// update this ahead of time in this special case
- output_prefix = proxy->output_prefix();
+ output_prefix = proxy->input_prefix();
cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n");
- output_prefix = output_prefix+".tmp";
- write_output_files();
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\
to:\n\
-\""+output_prefix+".colvars.state\"\n");
+\""+ proxy->input_prefix()+".colvars.state\"\n");
+ output_prefix = output_prefix+".tmp";
+ write_output_files();
cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
}
@@ -1120,6 +1127,7 @@ int colvarmodule::write_output_files()
bi != biases.end();
bi++) {
(*bi)->write_output_files();
+ (*bi)->write_state_to_replicas();
}
cvm::decrease_depth();
@@ -1212,20 +1220,30 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
<< " version " << std::string(COLVARS_VERSION) << "\n"
<< "}\n\n";
+ int error_code = COLVARS_OK;
+
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_restart(os);
+ error_code |= (*cvi)->write_output_files();
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
- (*bi)->write_restart(os);
+ (*bi)->write_state(os);
+ error_code |= (*bi)->write_state_to_replicas();
+ error_code |= (*bi)->write_output_files();
}
cvm::decrease_depth();
+ if (error_code != COLVARS_OK) {
+ // TODO make this function return an int instead
+ os.setstate(std::ios::failbit);
+ }
+
return os;
}
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index d5907f57e9a13a34872046de64ce610c659883e3..c63d04bcc432b1aef691aad7a8f8935612caa149 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -1,10 +1,17 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-27"
+#define COLVARS_VERSION "2016-12-22"
#endif
#ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index 86cd00ad29e504138f05e3959ca03fb38e62f6da..4d2e17f098dbb55457e5b1d8bc44f2cab79fdf70 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -1,3 +1,11 @@
+// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
#include <sstream>
diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h
index 892f1632e1f74925549ed975c3713de9a1e47b85..acdc46ddc9ee42d39f5889f2ebf4ab01c1781156 100644
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARPARSE_H
#define COLVARPARSE_H
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 255fba688bd179b62d8f916af27beba870aad243..bdf65c0b1f0a617a09427769e151a6db688dd154 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARPROXY_H
#define COLVARPROXY_H
@@ -24,11 +31,18 @@ public:
/// Pointer to the main object
colvarmodule *colvars;
- /// Default constructor
- inline colvarproxy() : b_smp_active(true), script(NULL) {}
+ /// Constructor
+ colvarproxy()
+ {
+ colvars = NULL;
+ b_simulation_running = true;
+ b_smp_active = true;
+ script = NULL;
+ }
- /// Default destructor
- virtual ~colvarproxy() {}
+ /// Destructor
+ virtual ~colvarproxy()
+ {}
/// (Re)initialize required member data after construction
virtual int setup()
@@ -104,6 +118,19 @@ public:
return 0;
}
+protected:
+
+ /// Whether a simulation is running (and try to prevent irrecovarable errors)
+ bool b_simulation_running;
+
+public:
+
+ /// Whether a simulation is running (and try to prevent irrecovarable errors)
+ virtual bool simulation_running() const
+ {
+ return b_simulation_running;
+ }
+
protected:
/// \brief Currently opened output files: by default, these are ofstream objects.
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index a0269aa7a911d8a1f3321e8df9bd40a6283308b0..ffb83e97801df0d4bb448b75747925e729fdb4ff 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <cstdlib>
#include <stdlib.h>
#include <string.h>
diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h
index 7778c3757f2d0c816bfed779bde658f1aed998a3..bfb8381e3b6c72bec02d3b6fe41e5037ab969bed 100644
--- a/lib/colvars/colvarscript.h
+++ b/lib/colvars/colvarscript.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARSCRIPT_H
#define COLVARSCRIPT_H
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 6f773bd5f74146398d555d5ab66156a35190d269..5200d4d0415838c52c4dd34e101b619312e2c43f 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <stdlib.h>
#include <string.h>
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 2c3d49f686ce54ac5a7cd029283e617c00b62ab1..e0cebb83bc8f1d049abb2bcec066b1f483e01fff 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARTYPES_H
#define COLVARTYPES_H
diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp
index 28650c825049c3a7de380e8a449ab640ed0efa1c..deccc6b7e0ebbbecb37a1497fbdb83c09443cce9 100644
--- a/lib/colvars/colvarvalue.cpp
+++ b/lib/colvars/colvarvalue.cpp
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <vector>
#include <sstream>
#include <iostream>
@@ -72,6 +79,7 @@ void colvarvalue::set_elem(int const icv, colvarvalue const &x)
void colvarvalue::set_random()
{
+ size_t ic;
switch (this->type()) {
case colvarvalue::type_scalar:
this->real_value = cvm::rand_gaussian();
@@ -91,7 +99,7 @@ void colvarvalue::set_random()
this->quaternion_value.q3 = cvm::rand_gaussian();
break;
case colvarvalue::type_vector:
- for (size_t ic = 0; ic < this->vector1d_value.size(); ic++) {
+ for (ic = 0; ic < this->vector1d_value.size(); ic++) {
this->vector1d_value[ic] = cvm::rand_gaussian();
}
break;
@@ -103,56 +111,6 @@ void colvarvalue::set_random()
}
-colvarvalue colvarvalue::inverse() const
-{
- switch (value_type) {
- case colvarvalue::type_scalar:
- return colvarvalue(1.0/real_value);
- break;
- case colvarvalue::type_3vector:
- case colvarvalue::type_unit3vector:
- case colvarvalue::type_unit3vectorderiv:
- return colvarvalue(cvm::rvector(1.0/rvector_value.x,
- 1.0/rvector_value.y,
- 1.0/rvector_value.z));
- break;
- case colvarvalue::type_quaternion:
- case colvarvalue::type_quaternionderiv:
- return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0,
- 1.0/quaternion_value.q1,
- 1.0/quaternion_value.q2,
- 1.0/quaternion_value.q3));
- break;
- case colvarvalue::type_vector:
- {
- cvm::vector1d<cvm::real> result(vector1d_value);
- if (elem_types.size() > 0) {
- // if we have information about non-scalar types, use it
- size_t i;
- for (i = 0; i < elem_types.size(); i++) {
- result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i],
- cvm::vector1d<cvm::real>((this->get_elem(i)).inverse()));
- }
- } else {
- size_t i;
- for (i = 0; i < result.size(); i++) {
- if (result[i] != 0.0) {
- result = 1.0/result[i];
- }
- }
- }
- return colvarvalue(result, type_vector);
- }
- break;
- case colvarvalue::type_notset:
- default:
- undef_op();
- break;
- }
- return colvarvalue();
-}
-
-
// binary operations between two colvarvalues
colvarvalue operator + (colvarvalue const &x1,
@@ -321,7 +279,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
(-1.0) * sin_t * v2.z +
cos_t/sin_t * (v1.z - cos_t*v2.z)
),
- colvarvalue::type_unit3vector );
+ colvarvalue::type_unit3vectorderiv );
}
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h
index 3e5e2645a3f530def77570d667e1cd42f4722c63..e369feefcd256a5a29868d7acdbea25ef2f86a24 100644
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@@ -1,5 +1,12 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#ifndef COLVARVALUE_H
#define COLVARVALUE_H
@@ -128,30 +135,23 @@ public:
{}
/// \brief Copy constructor from rvector base type (Note: this sets
- /// automatically a type \link type_3vector \endlink , if you want a
+ /// by default a type \link type_3vector \endlink , if you want a
/// \link type_unit3vector \endlink you must set it explicitly)
- inline colvarvalue(cvm::rvector const &v)
- : value_type(type_3vector), rvector_value(v)
- {}
-
- /// \brief Copy constructor from rvector base type (additional
- /// argument to make possible to choose a \link type_unit3vector
- /// \endlink
- inline colvarvalue(cvm::rvector const &v, Type const &vti)
+ inline colvarvalue(cvm::rvector const &v, Type vti = type_3vector)
: value_type(vti), rvector_value(v)
{}
/// \brief Copy constructor from quaternion base type
- inline colvarvalue(cvm::quaternion const &q)
- : value_type(type_quaternion), quaternion_value(q)
+ inline colvarvalue(cvm::quaternion const &q, Type vti = type_quaternion)
+ : value_type(vti), quaternion_value(q)
{}
+ /// Copy constructor from vector1d base type
+ colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti = type_vector);
+
/// Copy constructor from another \link colvarvalue \endlink
colvarvalue(colvarvalue const &x);
- /// Copy constructor from vector1d base type
- colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti);
-
/// Set to the null value for the data type currently defined
void reset();
@@ -211,10 +211,6 @@ public:
return std::sqrt(this->norm2());
}
- /// \brief Return the value whose scalar product with this value is
- /// 1
- inline colvarvalue inverse() const;
-
/// Square distance between this \link colvarvalue \endlink and another
cvm::real dist2(colvarvalue const &x2) const;
@@ -536,7 +532,7 @@ inline colvarvalue::colvarvalue(colvarvalue const &x)
}
}
-inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti)
+inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
{
if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
@@ -620,11 +616,22 @@ inline int colvarvalue::check_types_assign(colvarvalue::Type const &vt1,
}
if (vt1 != type_notset) {
- if (vt1 != vt2) {
- cvm::error("Trying to assign a colvar value with type \""+
- type_desc(vt2)+"\" to one with type \""+
- type_desc(vt1)+"\".\n");
- return COLVARS_ERROR;
+ if (((vt1 == type_unit3vector) &&
+ (vt2 == type_unit3vectorderiv)) ||
+ ((vt2 == type_unit3vector) &&
+ (vt1 == type_unit3vectorderiv)) ||
+ ((vt1 == type_quaternion) &&
+ (vt2 == type_quaternionderiv)) ||
+ ((vt2 == type_quaternion) &&
+ (vt1 == type_quaternionderiv))) {
+ return COLVARS_OK;
+ } else {
+ if (vt1 != vt2) {
+ cvm::error("Trying to assign a colvar value with type \""+
+ type_desc(vt2)+"\" to one with type \""+
+ type_desc(vt1)+"\".\n");
+ return COLVARS_ERROR;
+ }
}
}
return COLVARS_OK;
diff --git a/src/USER-COLVARS/README b/src/USER-COLVARS/README
index 759556aeb5c42a84564abf3becf4f319af8d84aa..9090a07ead1a28539031787209e874184993a52c 100644
--- a/src/USER-COLVARS/README
+++ b/src/USER-COLVARS/README
@@ -1,49 +1,46 @@
-This package implements the "fix colvars" command which can be used
+This package implements the "fix colvars" command, which can be used
in a LAMMPS input script.
-This fix allows use of "collective variables" to implement Adaptive
-Biasing Force, Metadynamics, Steered MD, Umbrella Sampling and
-Restraints.
+The fix allows use of enhanced sampling methods based on a reduced set
+of "collective variables", including free-energy estimators based on
+thermodynamic forces, non-equilibrium work and probability
+distributions.
-This package uses an external library in lib/colvars which must be
-compiled before making LAMMPS. See the lib/colvars/README file and
-the LAMMPS manual for information on building LAMMPS with external
-libraries. The settings in the Makefile.lammps file in that directory
-must be correct for LAMMPS to build correctly with this package
-installed.
+The package uses the "Colvars" library, whose source code is included
+in the LAMMPS source code distribution and must be linked with LAMMPS.
+See the lib/colvars/README file and the LAMMPS manual for information
+on building LAMMPS with external libraries. The settings in the
+Makefile.lammps file in that directory must be correct for LAMMPS to
+build correctly with this package installed.
-The external library is a portable collective variable module library
-written and maintained by Giacomo Fiorin (ICMS, Temple University,
-Philadelphia, PA, USA) and Jerome Henin (IBPC, CNRS, Paris, France).
+The files in the USER-COLVARS package folder implement an interface
+between LAMMPS and Colvars, originally written by Axel Kohlmeyer
+(akohlmey@gmail.com) and maintained by Giacomo Fiorin
+(giacomo.fiorin@gmail.com).
-More info about this library can be found at: http://colvars.github.io
+More info about the Colvars library can be found at:
-and in these publications:
+https://github.com/colvars/colvars
-Using collective variables to drive molecular dynamics simulations,
-Giacomo Fiorin, Michael L. Klein & Jérôme Hénin: Molecular Physics,
-111, 3345-3362 (2013)
-
-Exploring Multidimensional Free Energy Landscapes Using Time-Dependent
-Biases on Collective Variables, J. Hénin, G. Fiorin, C. Chipot, and
-M. L. Klein, J. Chem. Theory Comput., 6, 35-47 (2010).
+and in the reference article:
-The colvars fix implementes a thin interface layer, which exchanges
-information between LAMMPS and the collective variable module. This
-interface was written and is maintained by Axel Kohlmeyer
-(akohlmey@gmail.com)
-
-See the doc page of fix colvars for more details.
+Using collective variables to drive molecular dynamics simulations,
+G. Fiorin, M. L. Klein, and J. Henin,
+Molecular Physics 111, 3345 (2013)
+http://dx.doi.org/10.1080/00268976.2013.813594
-There is a reference manual for the package included with the LAMMPS
-doc pages: doc/PDF/colvars-refman-lammps.pdf
+A reference manual for the package and library is included with the
+LAMMPS doc pages:
+doc/PDF/colvars-refman-lammps.pdf
+which also includes citations to the articles documenting the various
+methods that make use Colvars.
-There are example scripts for using this package in
-examples/USER/colvars.
+There are also example scripts for using this package in the folder
+examples/USER/colvars, as well as the GitHub page for Colvars.
-The person who created this package is Axel Kohlmeyer at Temple U
-(akohlmey at gmail.com). Contact him directly if you have questions.
+Please contact Giacomo Fiorin (giacomo.fiorin@gmail.com) for questions
+regarding this package.
---------------------------------
-Version: 2015-01-08
+Version: 2016-12-22
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index 350009f6bf9488b3f225caebe9064821343e1129..7320263ba7b27f3d6a79caf40002934405df1e91 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -1,3 +1,12 @@
+// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
#include <mpi.h>
#include "lammps.h"
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index f1636e112bdec53bcd568e8fa0a06cead563fe59..c0b4ea20fecb37dfad428b0e7058b4fa885652bc 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -1,5 +1,13 @@
// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
#ifndef COLVARPROXY_LAMMPS_H
#define COLVARPROXY_LAMMPS_H
@@ -21,7 +29,7 @@
#endif
#ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2016-10-05"
+#define COLVARPROXY_VERSION "2016-12-27"
#endif
/* struct for packed data communication of coordinates and forces. */
diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp
index 6a3613d76a8bf70110b5f5f11098581b14ef5f78..59e6c46b76113bd1594dcf2acb983ee75212b15b 100644
--- a/src/USER-COLVARS/fix_colvars.cpp
+++ b/src/USER-COLVARS/fix_colvars.cpp
@@ -1,3 +1,12 @@
+// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h
index a1dae0d757c0c8a91d1e1da2ba3530c66d871976..c00b18aa4668b56ea8a1a82a729a0d02a9b0433a 100644
--- a/src/USER-COLVARS/fix_colvars.h
+++ b/src/USER-COLVARS/fix_colvars.h
@@ -1,3 +1,12 @@
+// -*- 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
diff --git a/src/USER-COLVARS/group_ndx.cpp b/src/USER-COLVARS/group_ndx.cpp
index 06b404098cd8e20061c063a5112a714af0e55f0a..50769180390cec50f8f693444393506ca990522c 100644
--- a/src/USER-COLVARS/group_ndx.cpp
+++ b/src/USER-COLVARS/group_ndx.cpp
@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
diff --git a/src/USER-COLVARS/group_ndx.h b/src/USER-COLVARS/group_ndx.h
index e0104f176fe55e54d2e9735ac1b68109c57335b6..25c38dcf9f069f54495f51150fe12d96e40f3905 100644
--- a/src/USER-COLVARS/group_ndx.h
+++ b/src/USER-COLVARS/group_ndx.h
@@ -1,4 +1,6 @@
-/* -*- c++ -*- ----------------------------------------------------------
+// -*- c++ -*-
+
+/* ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
diff --git a/src/USER-COLVARS/ndx_group.cpp b/src/USER-COLVARS/ndx_group.cpp
index 31d8332c9c7dec4cc12b092ac6b5aa9f22797750..c7972e372ec940989ec58d2c727cf2ce4aba5a78 100644
--- a/src/USER-COLVARS/ndx_group.cpp
+++ b/src/USER-COLVARS/ndx_group.cpp
@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
@@ -10,9 +12,8 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
-
/* ----------------------------------------------------------------------
- Contributing author: Axel Kohlmeyer (Temple U)
+ Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "ndx_group.h"
@@ -35,7 +36,7 @@ static char *find_section(FILE *fp, const char *name)
{
char linebuf[BUFLEN];
char *n,*p,*t,*r;
-
+
while ((p = fgets(linebuf,BUFLEN,fp))) {
t = strtok(p," \t\n\r\f");
if ((t != NULL) && *t == '[') {
diff --git a/src/USER-COLVARS/ndx_group.h b/src/USER-COLVARS/ndx_group.h
index 479c442f6125fd48981725d6d87b6d9d78b89be2..cd3250a1d5c33be73c1d30c0682d4c9f3d6a6f98 100644
--- a/src/USER-COLVARS/ndx_group.h
+++ b/src/USER-COLVARS/ndx_group.h
@@ -1,4 +1,6 @@
-/* -*- c++ -*- ----------------------------------------------------------
+// -*- c++ -*-
+
+/* ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov