From 51a148df161cdb42a984d51add3613eb7bd7e96c Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 28 Oct 2014 19:53:17 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12650
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
lib/colvars/colvar.cpp | 132 +++++++++-------
lib/colvars/colvar.h | 11 +-
lib/colvars/colvaratoms.cpp | 11 ++
lib/colvars/colvarbias.cpp | 18 +++
lib/colvars/colvarbias.h | 9 ++
lib/colvars/colvarbias_abf.cpp | 145 ++++++++++++++++++
lib/colvars/colvarbias_abf.h | 23 ++-
lib/colvars/colvarbias_meta.cpp | 1 +
lib/colvars/colvarcomp.cpp | 13 +-
lib/colvars/colvarcomp.h | 116 --------------
lib/colvars/colvarcomp_distances.cpp | 15 +-
lib/colvars/colvargrid.h | 64 ++++++++
lib/colvars/colvarmodule.cpp | 218 ++++++++++++++++++++-------
lib/colvars/colvarmodule.h | 44 +++++-
lib/colvars/colvarproxy.h | 30 +++-
lib/colvars/colvarscript.cpp | 5 +-
lib/colvars/colvartypes.cpp | 40 ++---
lib/colvars/colvartypes.h | 11 +-
18 files changed, 632 insertions(+), 274 deletions(-)
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 96231ec89b..aafa60e1fb 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -53,11 +53,17 @@ colvar::colvar (std::string const &conf)
if (cvcp != NULL) { \
cvcs.push_back (cvcp); \
cvcp->check_keywords (def_conf, def_config_key); \
+ if (cvm::get_error()) { \
+ cvm::error("Error: in setting up component \"" \
+ def_config_key"\".\n"); \
+ return; \
+ } \
cvm::decrease_depth(); \
} else { \
- cvm::error ("Error: in allocating component \"" \
+ cvm::error("Error: in allocating component \"" \
def_config_key"\".\n", \
MEMORY_ERROR); \
+ return; \
} \
if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) { \
if ( (cvcp->function_type != std::string ("distance_z")) && \
@@ -70,6 +76,7 @@ colvar::colvar (std::string const &conf)
"Period: "+cvm::to_str(cvcp->period) + \
" wrapAround: "+cvm::to_str(cvcp->wrap_center), \
INPUT_ERROR); \
+ return; \
} \
} \
if ( ! cvcs.back()->name.size()) \
@@ -150,63 +157,64 @@ colvar::colvar (std::string const &conf)
"", colvarparse::parse_silent)) {
enable(task_scripted);
- cvm::log("This colvar is a scripted function.");
-
- std::string type_str;
- get_keyval (conf, "scriptedFunctionType", type_str, "scalar");
+ cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
- x.type(colvarvalue::type_notset);
- for (i = 0; i < colvarvalue::type_all; i++) {
- if (type_str == colvarvalue::type_keyword[i]) {
- x.type(colvarvalue::Type(i));
- break;
- }
- }
- if (x.type() == colvarvalue::type_notset) {
- cvm::error("Could not parse scripted colvar type.");
- return;
- }
- x_reported.type (x.type());
- cvm::log(std::string("Expecting colvar value of type ")
- + colvarvalue::type_desc[x.type()]);
+ // Only accept scalar scripted colvars
+ // might accept other types when the infrastructure is in place
+ // for derivatives of vectors wrt vectors
+ x.type(colvarvalue::type_scalar);
+ x_reported.type(x.type());
- // Build ordered list of component values that will be
- // passed to the script
+ // Sort array of cvcs based on values of componentExp
+ std::vector<cvc *> temp_vec;
for (i = 1; i <= cvcs.size(); i++) {
for (j = 0; j < cvcs.size(); j++) {
if (cvcs[j]->sup_np == int(i)) {
- sorted_cvc_values.push_back(cvcs[j]->p_value());
+ temp_vec.push_back(cvcs[j]);
break;
}
}
}
- if (sorted_cvc_values.size() != cvcs.size()) {
- cvm::error("Could not find order numbers for all components"
+ if (temp_vec.size() != cvcs.size()) {
+ cvm::error("Could not find order numbers for all components "
"in componentExp values.");
return;
}
- }
+ cvcs = temp_vec;
- // this is set false if any of the components has an exponent
- // different from 1 in the polynomial
- b_linear = true;
+ // Build ordered list of component values that will be
+ // passed to the script
+ for (j = 0; j < cvcs.size(); j++) {
+ sorted_cvc_values.push_back(cvcs[j]->p_value());
+ }
+ }
- // these will be set to false if any of the cvcs has them false
- b_inverse_gradients = true;
- b_Jacobian_force = true;
+ if (!tasks[task_scripted]) {
+ // this is set false if any of the components has an exponent
+ // different from 1 in the polynomial
+ b_linear = true;
+ // these will be set to false if any of the cvcs has them false
+ b_inverse_gradients = true;
+ b_Jacobian_force = true;
+ }
+ // Test whether this is a single-component variable
// Decide whether the colvar is periodic
// Used to wrap extended DOF if extendedLagrangian is on
- if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1
- && (cvcs[0])->sup_coeff == 1.0 ) {
- this->b_periodic = true;
- this->period = (cvcs[0])->period;
+ if (cvcs.size() == 1 && (cvcs[0])->sup_np == 1
+ && (cvcs[0])->sup_coeff == 1.0
+ && !tasks[task_scripted]) {
+
+ b_single_cvc = true;
+ b_periodic = (cvcs[0])->b_periodic;
+ period = (cvcs[0])->period;
// TODO write explicit wrap() function for colvars to allow for
// sup_coeff different from 1
// this->period = (cvcs[0])->period * (cvcs[0])->sup_coeff;
} else {
- this->b_periodic = false;
- this->period = 0.0;
+ b_single_cvc = false;
+ b_periodic = false;
+ period = 0.0;
}
// check the available features of each cvc
@@ -808,13 +816,14 @@ void colvar::calc()
// each atom group will take care of its own ref_pos_group, if defined
}
}
- if (tasks[task_output_velocity]) {
- for (i = 0; i < cvcs.size(); i++) {
- for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
- cvcs[i]->atom_groups[ig]->read_velocities();
- }
- }
- }
+//// Don't try to get atom velocities, as no back-end currently implements it
+// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
+// for (i = 0; i < cvcs.size(); i++) {
+// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
+// cvcs[i]->atom_groups[ig]->read_velocities();
+// }
+// }
+// }
if (tasks[task_system_force]) {
for (i = 0; i < cvcs.size(); i++) {
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
@@ -1153,8 +1162,11 @@ void colvar::communicate_forces()
for (i = 0; i < cvcs.size(); i++) {
cvm::increase_depth();
- // Note: we need a dot product here
- (cvcs[i])->apply_force (f * func_grads[i]);
+ // Force is scalar times colvarvalue (scalar or vector)
+ // Note: this can only handle scalar colvars (scalar values of f)
+ // A non-scalar colvar would need the gradient to be expressed
+ // as an order-2 tensor
+ (cvcs[i])->apply_force (f.real_value * func_grads[i]);
cvm::decrease_depth();
}
} else if (x.type() == colvarvalue::type_scalar) {
@@ -1223,30 +1235,38 @@ bool colvar::periodic_boundaries() const
cvm::real colvar::dist2 (colvarvalue const &x1,
colvarvalue const &x2) const
{
- return (cvcs[0])->dist2 (x1, x2);
+ if (b_single_cvc) {
+ return (cvcs[0])->dist2(x1, x2);
+ } else {
+ return x1.dist2(x2);
+ }
}
colvarvalue colvar::dist2_lgrad (colvarvalue const &x1,
colvarvalue const &x2) const
{
- return (cvcs[0])->dist2_lgrad (x1, x2);
+ if (b_single_cvc) {
+ return (cvcs[0])->dist2_lgrad (x1, x2);
+ } else {
+ return x1.dist2_grad(x2);
+ }
}
colvarvalue colvar::dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const
{
- return (cvcs[0])->dist2_rgrad (x1, x2);
-}
-
-cvm::real colvar::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return (cvcs[0])->compare (x1, x2);
+ if (b_single_cvc) {
+ return (cvcs[0])->dist2_rgrad (x1, x2);
+ } else {
+ return x2.dist2_grad(x1);
+ }
}
void colvar::wrap (colvarvalue &x) const
{
- (cvcs[0])->wrap (x);
+ if (b_single_cvc) {
+ (cvcs[0])->wrap (x);
+ }
return;
}
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 896613d7da..3d06cf8df5 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -83,6 +83,10 @@ public:
/// combination of \link cvc \endlink elements
bool b_linear;
+ /// \brief True if this \link colvar \endlink is equal to
+ /// its only constituent cvc
+ bool b_single_cvc;
+
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating inverse gradients
bool b_inverse_gradients;
@@ -344,13 +348,6 @@ public:
colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- /// \brief Use the internal metrics (as from \link cvc
- /// \endlink objects) to compare colvar values
- ///
- /// Handles correctly symmetries and periodic boundary conditions
- cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
-
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval
///
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index f56871d614..250d23acfa 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -138,6 +138,7 @@ int cvm::atom_group::parse (std::string const &conf,
for (size_t i = 0; i < atom_indexes.size(); i++) {
this->push_back (cvm::atom (atom_indexes[i]));
}
+ if (cvm::get_error()) return COLVARS_ERROR;
} else {
cvm::error ("Error: no numbers provided for \""
"atomNumbers\".\n", INPUT_ERROR);
@@ -166,6 +167,7 @@ int cvm::atom_group::parse (std::string const &conf,
for (size_t i = 0; i < index_groups_i->size(); i++) {
this->push_back (cvm::atom ((*index_groups_i)[i]));
}
+ if (cvm::get_error()) return COLVARS_ERROR;
}
}
@@ -185,6 +187,7 @@ int cvm::atom_group::parse (std::string const &conf,
for (int anum = initial; anum <= final; anum++) {
this->push_back (cvm::atom (anum));
}
+ if (cvm::get_error()) return COLVARS_ERROR;
range_conf = "";
continue;
}
@@ -235,6 +238,7 @@ int cvm::atom_group::parse (std::string const &conf,
for (int resid = initial; resid <= final; resid++) {
this->push_back (cvm::atom (resid, atom_name, psf_segid));
}
+ if (cvm::get_error()) return COLVARS_ERROR;
range_conf = "";
} else {
cvm::error ("Error: cannot parse definition for \""
@@ -272,6 +276,9 @@ int cvm::atom_group::parse (std::string const &conf,
}
}
+ // Catch any errors from all the initialization steps above
+ if (cvm::get_error()) return COLVARS_ERROR;
+
for (std::vector<cvm::atom>::iterator a1 = this->begin();
a1 != this->end(); ++a1) {
std::vector<cvm::atom>::iterator a2 = a1;
@@ -433,6 +440,10 @@ int cvm::atom_group::parse (std::string const &conf,
std::string (key)+"\".\n");
this->check_keywords (group_conf, key);
+ if (cvm::get_error()) {
+ cvm::error("Error setting up atom group \""+std::string (key)+"\".");
+ return COLVARS_ERROR;
+ }
cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
cvm::to_str (this->size())+" atoms initialized: total mass = "+
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index d88611b9b9..7c21c551c8 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -124,6 +124,24 @@ cvm::real colvarbias::energy_difference(std::string const &conf)
}
+// So far, these are only implemented in colvarsbias_abf
+int colvarbias::bin_num()
+{
+ cvm::error ("Error: bin_num() not implemented.\n");
+ return -1;
+}
+int colvarbias::current_bin()
+{
+ cvm::error ("Error: current_bin() not implemented.\n");
+ return -1;
+}
+int colvarbias::bin_count(int bin_index)
+{
+ cvm::error ("Error: bin_count() not implemented.\n");
+ return -1;
+}
+
+
std::ostream & colvarbias::write_traj_label (std::ostream &os)
{
os << " ";
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index c3624784d1..ad24f40f30 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -29,6 +29,15 @@ public:
/// Calculate change in energy from using alternate configuration
virtual cvm::real energy_difference(std::string const &conf);
+ /// Give the total number of bins for a given bias.
+ virtual int bin_num();
+ /// Calculate the bin index for a given bias.
+ virtual int current_bin();
+ //// Give the count at a given bin index.
+ virtual int bin_count(int bin_index);
+ //// Share information between replicas, whatever it may be.
+ virtual void replica_share() {};
+
/// Perform analysis tasks
virtual inline void analyse() {}
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index cf5825650c..812490013c 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -44,6 +44,18 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
get_keyval (conf, "historyFreq", history_freq, 0);
b_history_files = (history_freq > 0);
+ // shared ABF
+ get_keyval (conf, "shared", shared_on, false);
+ if (shared_on) {
+ if (!cvm::replica_enabled || cvm::replica_num() <= 1)
+ cvm::error ("Error: shared ABF requires more than one replica.");
+ else
+ cvm::log ("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n");
+
+ // If shared_freq is not set, we default to output_freq
+ get_keyval (conf, "sharedFreq", shared_freq, output_freq);
+ }
+
// ************* checking the associated colvars *******************
if (colvars.size() == 0) {
@@ -109,6 +121,15 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
gradients->samples = samples;
samples->has_parent_data = true;
+ // For shared ABF, we store a second set of grids.
+ // This used to be only if "shared" was defined,
+ // but now we allow calling share externally (e.g. from Tcl).
+ last_samples = new colvar_grid_count (colvars);
+ last_gradients = new colvar_grid_gradient (colvars);
+ last_gradients->samples = last_samples;
+ last_samples->has_parent_data = true;
+ shared_last_step = -1;
+
// If custom grids are provided, read them
if ( input_prefix.size() > 0 ) {
read_gradients_samples ();
@@ -130,6 +151,19 @@ colvarbias_abf::~colvarbias_abf()
gradients = NULL;
}
+ // shared ABF
+ // We used to only do this if "shared" was defined,
+ // but now we can call shared externally
+ if (last_samples) {
+ delete last_samples;
+ last_samples = NULL;
+ }
+
+ if (last_gradients) {
+ delete last_gradients;
+ last_gradients = NULL;
+ }
+
delete [] force;
if (cvm::n_abf_biases > 0)
@@ -225,13 +259,103 @@ cvm::real colvarbias_abf::update()
write_gradients_samples (output_prefix);
}
if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
+ cvm::log ("ABFHISTORYFILE "+cvm::to_str(cvm::step_absolute()));
// append to existing file only if cvm::step_absolute() > 0
// otherwise, backup and replace
write_gradients_samples (output_prefix + ".hist", (cvm::step_absolute() > 0));
}
+
+ if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
+ // Share gradients and samples for shared ABF.
+ replica_share();
+ }
+
+ // Prepare for the first sharing.
+ if (shared_last_step < 0) {
+ // Copy the current gradient and count values into last.
+ last_gradients->copy_grid(*gradients);
+ last_samples->copy_grid(*samples);
+ shared_last_step = cvm::step_absolute();
+ cvm::log ("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
+ }
+
return 0.0;
}
+void colvarbias_abf::replica_share () {
+ int p;
+
+ if( !cvm::replica_enabled() ) {
+ cvm::error ("Error: shared ABF: No replicas.\n");
+ return;
+ }
+ // We must have stored the last_gradients and last_samples.
+ if (shared_last_step < 0 ) {
+ cvm::error ("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n");
+ return;
+ }
+
+ // Share gradients for shared ABF.
+ cvm::log ("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) );
+
+ // Count of data items.
+ size_t data_n = gradients->raw_data_num();
+ size_t samp_start = data_n*sizeof(cvm::real);
+ size_t msg_total = data_n*sizeof(size_t) + samp_start;
+ char* msg_data = new char[msg_total];
+
+ if (cvm::replica_index() == 0) {
+ // Replica 0 collects the delta gradient and count from the others.
+ for (p = 1; p < cvm::replica_num(); p++) {
+ // Receive the deltas.
+ cvm::replica_comm_recv(msg_data, msg_total, p);
+
+ // Map the deltas from the others into the grids.
+ last_gradients->raw_data_in((cvm::real*)(&msg_data[0]));
+ last_samples->raw_data_in((size_t*)(&msg_data[samp_start]));
+
+ // Combine the delta gradient and count of the other replicas
+ // with Replica 0's current state (including its delta).
+ gradients->add_grid( *last_gradients );
+ samples->add_grid( *last_samples );
+ }
+
+ // Now we must send the combined gradient to the other replicas.
+ gradients->raw_data_out((cvm::real*)(&msg_data[0]));
+ samples->raw_data_out((size_t*)(&msg_data[samp_start]));
+ for (p = 1; p < cvm::replica_num(); p++) {
+ cvm::replica_comm_send(msg_data, msg_total, p);
+ }
+
+ } else {
+ // All other replicas send their delta gradient and count.
+ // Calculate the delta gradient and count.
+ last_gradients->delta_grid (*gradients);
+ last_samples->delta_grid (*samples);
+
+ // Cast the raw char data to the gradient and samples.
+ last_gradients->raw_data_out((cvm::real*)(&msg_data[0]));
+ last_samples->raw_data_out((size_t*)(&msg_data[samp_start]));
+ cvm::replica_comm_send(msg_data, msg_total, 0);
+
+ // We now receive the combined gradient from Replica 0.
+ cvm::replica_comm_recv(msg_data, msg_total, 0);
+ // We sync to the combined gradient computed by Replica 0.
+ gradients->raw_data_in((cvm::real*)(&msg_data[0]));
+ samples->raw_data_in((size_t*)(&msg_data[samp_start]));
+ }
+
+ // Without a barrier it's possible that one replica starts
+ // share 2 when other replicas haven't finished share 1.
+ cvm::replica_comm_barrier();
+ // Done syncing the replicas.
+ delete[] msg_data;
+
+ // Copy the current gradient and count values into last.
+ last_gradients->copy_grid(*gradients);
+ last_samples->copy_grid(*samples);
+ shared_last_step = cvm::step_absolute();
+}
void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool append)
{
@@ -268,6 +392,27 @@ void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool ap
return;
}
+
+// For Tcl implementation of selection rules.
+/// Give the total number of bins for a given bias.
+int colvarbias_abf::bin_num() {
+ return samples->number_of_points(0);
+}
+/// Calculate the bin index for a given bias.
+int colvarbias_abf::current_bin() {
+ return samples->current_bin_scalar(0);
+}
+/// Give the count at a given bin index.
+int colvarbias_abf::bin_count(int bin_index) {
+ if (bin_index < 0 || bin_index >= bin_num()) {
+ cvm::error ("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index));
+ return -1;
+ }
+ std::vector<int> ix(1,(int)bin_index);
+ return samples->value(ix);
+}
+
+
void colvarbias_abf::read_gradients_samples ()
{
std::string samples_in_name, gradients_in_name;
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index 6873a8f9b1..b89ca68dec 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -1,5 +1,3 @@
-/// -*- c++ -*-
-
/************************************************************************
* Headers for the ABF and histogram biases *
************************************************************************/
@@ -11,7 +9,6 @@
#include <list>
#include <sstream>
#include <iomanip>
-//#include <cmath>
#include "colvarbias.h"
#include "colvargrid.h"
@@ -62,8 +59,28 @@ private:
/// n-dim grid of number of samples
colvar_grid_count *samples;
+ // shared ABF
+ bool shared_on;
+ size_t shared_freq;
+ int shared_last_step;
+ // Share between replicas -- may be called independently of update
+ virtual void replica_share();
+
+ // Store the last set for shared ABF
+ colvar_grid_gradient *last_gradients;
+ colvar_grid_count *last_samples;
+
+ // For Tcl implementation of selection rules.
+ /// Give the total number of bins for a given bias.
+ virtual int bin_num();
+ /// Calculate the bin index for a given bias.
+ virtual int current_bin();
+ //// Give the count at a given bin index.
+ virtual int bin_count(int bin_index);
+
/// Write human-readable FE gradients and sample count
void write_gradients_samples (const std::string &prefix, bool append = false);
+ void write_last_gradients_samples (const std::string &prefix, bool append = false);
/// Read human-readable FE gradients and sample count (if not using restart)
void read_gradients_samples ();
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index b49aee87c0..5e4e002ded 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -723,6 +723,7 @@ void colvarbias_meta::calc_hills_force (size_t const &i,
}
break;
+ case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
for (h = h_first; h != h_last; h++) {
if (h->value() == 0.0) continue;
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index f333f5db3e..949256004e 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -47,11 +47,16 @@ void colvar::cvc::parse_group (std::string const &conf,
bool optional)
{
if (key_lookup (conf, group_key)) {
- group.parse (conf, group_key);
+ if (group.parse (conf, group_key) != COLVARS_OK) {
+ cvm::error ("Error parsing definition for atom group \""+
+ std::string (group_key)+"\".\n");
+ return;
+ }
} else {
if (! optional) {
- cvm::fatal_error ("Error: definition for atom group \""+
+ cvm::error ("Error: definition for atom group \""+
std::string (group_key)+"\" not found.\n");
+ return;
}
}
}
@@ -137,9 +142,7 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group)
21, 14)+"\n");
cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) /
- std::fabs (x_1 - x_0),
- 12, 5)+
- ".\n");
+ std::fabs (x_1 - x_0), 12, 5)+"\n");
}
}
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 253c1d8de7..735c2f5bc1 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -219,14 +219,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- /// \brief Return a positive number if x2>x1, zero if x2==x1,
- /// negative otherwise (can be redefined to transparently implement
- /// constraints, symmetries and periodicities) \b Note: \b it \b
- /// only \b works \b with \b scalar \b variables, otherwise raises
- /// an error
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
-
/// \brief Wrapp value (for periodic/symmetric cvcs)
virtual void wrap (colvarvalue &x) const;
@@ -300,18 +292,6 @@ inline colvarvalue colvar::cvc::dist2_rgrad (colvarvalue const &x1,
return x2.dist2_grad (x1);
}
-inline cvm::real colvar::cvc::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- if (this->type() == colvarvalue::type_scalar) {
- return cvm::real (x1 - x2);
- } else {
- cvm::error ("Error: you requested an operation which requires "
- "comparison between two non-scalar values.\n");
- return 0.0;
- }
-}
-
inline void colvar::cvc::wrap (colvarvalue &x) const
{
return;
@@ -356,8 +336,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -383,9 +361,6 @@ public:
/// Redefined to handle the box periodicity
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- /// Redefined to handle the box periodicity
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -408,8 +383,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -453,8 +426,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
/// \brief Redefined to make use of the user-provided period
virtual void wrap (colvarvalue &x) const;
};
@@ -485,8 +456,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -511,8 +480,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -541,8 +508,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -565,8 +530,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -592,8 +555,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -635,8 +596,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -685,8 +644,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -736,9 +693,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
- /// Redefined to handle the 2*PI periodicity
virtual void wrap (colvarvalue &x) const;
};
@@ -796,8 +750,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
/// \brief Colvar component: self-coordination number within a group
@@ -835,8 +787,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
/// \brief Colvar component: hydrogen bond, defined as the product of
@@ -872,8 +822,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -915,8 +863,6 @@ public:
// colvarvalue const &x2) const;
// virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
// colvarvalue const &x2) const;
-// virtual cvm::real compare (colvarvalue const &x1,
-// colvarvalue const &x2) const;
// };
@@ -958,8 +904,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
/// \brief Colvar component: dihedPC
@@ -988,8 +932,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
/// \brief Colvar component: orientation in space of an atom group,
@@ -1030,8 +972,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -1055,8 +995,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -1080,8 +1018,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -1108,8 +1044,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -1141,9 +1075,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
- /// Redefined to handle the 2*PI periodicity
virtual void wrap (colvarvalue &x) const;
};
@@ -1180,8 +1111,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
- virtual cvm::real compare (colvarvalue const &x1,
- colvarvalue const &x2) const;
};
@@ -1214,12 +1143,6 @@ public:
return this->dist2_lgrad (x2, x1); \
} \
\
- inline cvm::real colvar::TYPE::compare (colvarvalue const &x1, \
- colvarvalue const &x2) const \
- { \
- return this->dist2_lgrad (x1, x2); \
- } \
- \
simple_scalar_dist_functions (distance)
// NOTE: distance_z has explicit functions, see below
@@ -1268,12 +1191,6 @@ inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
-inline cvm::real colvar::dihedral::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return dist2_lgrad (x1, x2);
-}
-
inline void colvar::dihedral::wrap (colvarvalue &x) const
{
if ((x.real_value - wrap_center) >= 180.0) {
@@ -1313,12 +1230,6 @@ inline colvarvalue colvar::spin_angle::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
-inline cvm::real colvar::spin_angle::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return dist2_lgrad (x1, x2);
-}
-
inline void colvar::spin_angle::wrap (colvarvalue &x) const
{
if ((x.real_value - wrap_center) >= 180.0) {
@@ -1370,12 +1281,6 @@ inline colvarvalue colvar::distance_z::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
-inline cvm::real colvar::distance_z::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- return dist2_lgrad (x1, x2);
-}
-
inline void colvar::distance_z::wrap (colvarvalue &x) const
{
if (! this->b_periodic) {
@@ -1412,13 +1317,6 @@ inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1,
return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
}
-inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::error ("Error: cannot compare() two distance vectors.\n");
- return 0.0;
-}
-
inline cvm::real colvar::distance_dir::dist2 (colvarvalue const &x1,
colvarvalue const &x2) const
{
@@ -1437,13 +1335,6 @@ inline colvarvalue colvar::distance_dir::dist2_rgrad (colvarvalue const &x1,
return colvarvalue ((x2.rvector_value - x1.rvector_value), colvarvalue::type_unitvector);
}
-inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::error ("Error: cannot compare() two distance directions.\n");
- return 0.0;
-}
-
// distance between quaternions
inline cvm::real colvar::orientation::dist2 (colvarvalue const &x1,
@@ -1464,12 +1355,5 @@ inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1,
return x2.quaternion_value.dist2_grad (x1);
}
-inline cvm::real colvar::orientation::compare (colvarvalue const &x1,
- colvarvalue const &x2) const
-{
- cvm::error ("Error: cannot compare() two quaternions.\n");
- return 0.0;
-}
-
#endif
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 3e1163db4b..43bfebf188 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -725,7 +725,7 @@ colvar::rmsd::rmsd (std::string const &conf)
return;
}
- if (atoms.ref_pos_group != NULL) {
+ if (atoms.ref_pos_group != NULL && b_Jacobian_derivative) {
cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
"Jacobian derivatives of the RMSD will not be calculated.\n");
b_Jacobian_derivative = false;
@@ -788,10 +788,21 @@ colvar::rmsd::rmsd (std::string const &conf)
atoms.ref_pos = ref_pos;
atoms.center_ref_pos();
+ cvm::log ("This is a standard minimum RMSD, derivatives of the optimal rotation "
+ "will not be computed as they cancel out in the gradients.");
+ atoms.b_fit_gradients = false;
+ }
+
+ if (atoms.b_rotate) {
+ // TODO: finer-grained control of this would require exposing a
+ // "request_Jacobian_derivative()" method to the colvar, and the same
+ // from the colvar to biases
+ // TODO: this should not be enabled here anyway, as it is not specific of the
+ // component - instead it should be decided in a generic way by the atom group
+
// request the calculation of the derivatives of the rotation defined by the atom group
atoms.rot.request_group1_gradients (atoms.size());
// request derivatives of optimal rotation wrt reference coordinates for Jacobian:
- // this is only required for ABF, but we do both groups here for better caching
atoms.rot.request_group2_gradients (atoms.size());
}
}
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index c3a61dcaa0..7cc255ee74 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -362,6 +362,70 @@ public:
has_data = true;
}
+ /// \brief Get the change from this to other_grid
+ /// and store the result in this.
+ /// this_grid := other_grid - this_grid
+ /// Grids must have the same dimensions.
+ void delta_grid (colvar_grid<T> const &other_grid)
+ {
+
+ if (other_grid.multiplicity() != this->multiplicity()) {
+ cvm::error ("Error: trying to subtract two grids with "
+ "different multiplicity.\n");
+ return;
+ }
+
+ if (other_grid.data.size() != this->data.size()) {
+ cvm::error ("Error: trying to subtract two grids with "
+ "different size.\n");
+ return;
+ }
+
+ for (size_t i = 0; i < data.size(); i++) {
+ data[i] = other_grid.data[i] - data[i];
+ }
+ has_data = true;
+ }
+
+ /// \brief Copy data from another grid of the same type, AND
+ /// identical definition (boundaries, widths)
+ /// Added for shared ABF.
+ void copy_grid (colvar_grid<T> const &other_grid)
+ {
+ if (other_grid.multiplicity() != this->multiplicity()) {
+ cvm::error ("Error: trying to copy two grids with "
+ "different multiplicity.\n");
+ return;
+ }
+
+ if (other_grid.data.size() != this->data.size()) {
+ cvm::error ("Error: trying to copy two grids with "
+ "different size.\n");
+ return;
+ }
+
+
+ for (size_t i = 0; i < data.size(); i++) {
+ data[i] = other_grid.data[i];
+ }
+ has_data = true;
+ }
+
+ /// \brief Extract the grid data as they are represented in memory.
+ /// Put the results in "out_data".
+ void raw_data_out (T* out_data) const
+ {
+ for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
+ }
+ /// \brief Input the data as they are represented in memory.
+ void raw_data_in (const T* in_data)
+ {
+ for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
+ has_data = true;
+ }
+ /// \brief Size of the data as they are represented in memory.
+ size_t raw_data_num() const { return data.size(); }
+
/// \brief Get the binned value indexed by ix, or the first of them
/// if the multiplicity is larger than 1
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 3b85273a60..03f9ba704c 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -99,7 +99,7 @@ int colvarmodule::config (std::string &conf)
// parse global options
error_code |= parse_global_params (conf);
- if (error_code != COLVARS_OK) {
+ if (error_code != COLVARS_OK || cvm::get_error()) {
set_error_bits(INPUT_ERROR);
return COLVARS_ERROR;
}
@@ -107,13 +107,23 @@ int colvarmodule::config (std::string &conf)
// parse the options for collective variables
error_code |= parse_colvars (conf);
+ if (error_code != COLVARS_OK || cvm::get_error()) {
+ set_error_bits(INPUT_ERROR);
+ return COLVARS_ERROR;
+ }
+
// parse the options for biases
error_code |= parse_biases (conf);
+ if (error_code != COLVARS_OK || cvm::get_error()) {
+ set_error_bits(INPUT_ERROR);
+ return COLVARS_ERROR;
+ }
+
// done parsing known keywords, check that all keywords found were valid ones
error_code |= parse->check_keywords (conf, "colvarmodule");
- if (error_code != COLVARS_OK) {
+ if (error_code != COLVARS_OK || cvm::get_error()) {
set_error_bits(INPUT_ERROR);
return COLVARS_ERROR;
}
@@ -179,12 +189,11 @@ int colvarmodule::parse_colvars (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
colvars.push_back (new colvar (colvar_conf));
- if (cvm::get_error()) {
- delete colvars.back();
- colvars.pop_back();
- return COLVARS_ERROR;
- }
- if ((colvars.back())->check_keywords (colvar_conf, "colvar") != COLVARS_OK) {
+ if (cvm::get_error() ||
+ ((colvars.back())->check_keywords (colvar_conf, "colvar") != COLVARS_OK)) {
+ cvm::log("Error while constructing colvar number " +
+ cvm::to_str(colvars.size()) + " : deleting.");
+ delete colvars.back(); // the colvar destructor updates the colvars array
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -192,6 +201,7 @@ int colvarmodule::parse_colvars (std::string const &conf)
cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
+ cvm::decrease_depth();
colvar_conf = "";
}
@@ -208,6 +218,17 @@ int colvarmodule::parse_colvars (std::string const &conf)
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
+bool colvarmodule::check_new_bias(std::string &conf, char const *key)
+{
+ if (cvm::get_error() ||
+ (biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
+ cvm::log("Error while constructing bias number " +
+ cvm::to_str(biases.size()) + " : deleting.\n");
+ delete biases.back(); // the bias destructor updates the biases array
+ return true;
+ }
+ return false;
+}
int colvarmodule::parse_biases (std::string const &conf)
{
@@ -223,12 +244,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_abf (abf_conf, "abf"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (abf_conf, "abf") != COLVARS_OK) {
+ if (cvm::check_new_bias(abf_conf, "abf")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -250,12 +266,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_restraint_harmonic (harm_conf, "harmonic"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (harm_conf, "harmonic") != COLVARS_OK) {
+ if (cvm::check_new_bias(harm_conf, "harmonic")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -277,12 +288,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_restraint_linear (lin_conf, "linear"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (lin_conf, "linear") != COLVARS_OK) {
+ if (cvm::check_new_bias(lin_conf, "linear")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -304,12 +310,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_alb (alb_conf, "ALB"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (alb_conf, "ALB") != COLVARS_OK) {
+ if (cvm::check_new_bias(alb_conf, "ALB")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -332,12 +333,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_histogram (histo_conf, "histogram"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (histo_conf, "histogram") != COLVARS_OK) {
+ if (cvm::check_new_bias(histo_conf, "histogram")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -359,12 +355,7 @@ int colvarmodule::parse_biases (std::string const &conf)
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_meta (meta_conf, "metadynamics"));
- if (cvm::get_error()) {
- delete biases.back();
- biases.pop_back();
- return COLVARS_ERROR;
- }
- if ((biases.back())->check_keywords (meta_conf, "metadynamics") != COLVARS_OK) {
+ if (cvm::check_new_bias(meta_conf, "metadynamics")) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
@@ -444,7 +435,6 @@ std::string colvarmodule::read_colvar(std::string const &name)
return ss.str();
}
-
cvm::real colvarmodule::energy_difference (std::string const &bias_name,
std::string const &conf)
{
@@ -458,6 +448,125 @@ cvm::real colvarmodule::energy_difference (std::string const &bias_name,
return energy_diff;
}
+// For shared ABF
+cvm::real colvarmodule::read_width(std::string const &name)
+{
+ cvm::increase_depth();
+ int found = 0;
+ real width = -1.0;
+ for (std::vector<colvar *>::iterator cvi = colvars.begin();
+ cvi != colvars.end();
+ cvi++) {
+ if ( (*cvi)->name == name ) {
+ ++found;
+ width = (*cvi)->width;
+ }
+ }
+ if (found < 1) {
+ cvm::error ("Error: bias not found.\n");
+ } else if (found > 1) {
+ cvm::error ("Error: duplicate bias name.\n");
+ } else {
+ cvm::decrease_depth();
+ }
+ return width;
+}
+
+size_t colvarmodule::bias_current_bin (std::string const &bias_name)
+{
+ cvm::increase_depth();
+ int found = 0;
+ size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
+
+ for (std::vector<colvarbias *>::iterator bi = biases.begin();
+ bi != biases.end();
+ bi++) {
+ if ( (*bi)->name == bias_name ) {
+ ++found;
+ ret = (*bi)->current_bin ();
+ }
+ }
+ if (found < 1) {
+ cvm::error ("Error: bias not found.\n");
+ } else if (found > 1) {
+ cvm::error ("Error: duplicate bias name.\n");
+ } else {
+ cvm::decrease_depth();
+ }
+ return ret;
+}
+
+size_t colvarmodule::bias_bin_num (std::string const &bias_name)
+{
+ cvm::increase_depth();
+ int found = 0;
+ size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
+
+ for (std::vector<colvarbias *>::iterator bi = biases.begin();
+ bi != biases.end();
+ bi++) {
+ if ( (*bi)->name == bias_name ) {
+ ++found;
+ ret = (*bi)->bin_num ();
+ }
+ }
+ if (found < 1) {
+ cvm::error ("Error: bias not found.\n");
+ } else if (found > 1) {
+ cvm::error ("Error: duplicate bias name.\n");
+ } else {
+ cvm::decrease_depth();
+ }
+ return ret;
+}
+
+size_t colvarmodule::bias_bin_count (std::string const &bias_name, size_t bin_index)
+{
+ cvm::increase_depth();
+ int found = 0;
+ size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
+
+ for (std::vector<colvarbias *>::iterator bi = biases.begin();
+ bi != biases.end();
+ bi++) {
+ if ( (*bi)->name == bias_name ) {
+ ++found;
+ ret = (*bi)->bin_count (bin_index);
+ }
+ }
+ if (found < 1) {
+ cvm::error ("Error: bias not found.\n");
+ } else if (found > 1) {
+ cvm::error ("Error: duplicate bias name.\n");
+ } else {
+ cvm::decrease_depth();
+ }
+ return ret;
+}
+
+void colvarmodule::bias_share (std::string const &bias_name)
+{
+ cvm::increase_depth();
+ int found = 0;
+
+ for (std::vector<colvarbias *>::iterator bi = biases.begin();
+ bi != biases.end();
+ bi++) {
+ if ( (*bi)->name == bias_name ) {
+ ++found;
+ (*bi)->replica_share ();
+ }
+ }
+ if (found < 1) {
+ cvm::error ("Error: bias not found.\n");
+ return;
+ }
+ if (found > 1) {
+ cvm::error ("Error: duplicate bias name.\n");
+ return;
+ }
+ cvm::decrease_depth();
+}
int colvarmodule::calc() {
@@ -675,21 +784,20 @@ colvarmodule::~colvarmodule()
int colvarmodule::reset()
{
- if (cvm::debug())
- cvm::log ("colvars::reset() called.\n");
- for (std::vector<colvar *>::iterator cvi = colvars.begin();
- cvi != colvars.end();
+ cvm::log("Resetting the Collective Variables Module.\n");
+ // Iterate backwards because we are deleting the elements as we go
+ for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
+ cvi != colvars.rend();
cvi++) {
- delete *cvi;
- cvi--;
+ delete *cvi; // the colvar destructor updates the colvars array
}
colvars.clear();
- for (std::vector<colvarbias *>::iterator bi = biases.begin();
- bi != biases.end();
+ // Iterate backwards because we are deleting the elements as we go
+ for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
+ bi != biases.rend();
bi++) {
- delete *bi;
- bi--;
+ delete *bi; // the bias destructor updates the biases array
}
biases.clear();
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index fab6ae5802..ba1a4937e9 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2014-10-07"
+#define COLVARS_VERSION "2014-10-21"
#endif
#ifndef COLVARS_DEBUG
@@ -219,6 +219,9 @@ public:
/// Parse and initialize collective variable biases
int parse_biases (std::string const &conf);
+ /// Test error condition and keyword parsing
+ /// on error, delete new bias
+ bool check_new_bias(std::string &conf, char const *key);
// "Setup" functions (change internal data based on related data
// from the proxy that may change during program execution)
@@ -270,6 +273,17 @@ public:
/// currently works for harmonic (force constant and/or centers)
real energy_difference (std::string const &bias_name, std::string const &conf);
+ /// Give the bin width in the units of the colvar.
+ real read_width(std::string const &name);
+ /// Give the total number of bins for a given bias.
+ size_t bias_bin_num(std::string const &bias_name);
+ /// Calculate the bin index for a given bias.
+ size_t bias_current_bin(std::string const &bias_name);
+ //// Give the count at a given bin index.
+ size_t bias_bin_count(std::string const &bias_name, size_t bin_index);
+ //// Share among replicas.
+ void bias_share(std::string const &bias_name);
+
/// Calculate collective variables and biases
int calc();
@@ -346,6 +360,13 @@ public:
/// Print a message to the main log and exit normally
static void exit (std::string const &message);
+ // Replica exchange commands.
+ static bool replica_enabled();
+ static int replica_index();
+ static int replica_num();
+ static void replica_comm_barrier();
+ static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
+ static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
/// \brief Get the distance between two atomic positions with pbcs handled
/// correctly
@@ -533,6 +554,27 @@ inline cvm::real cvm::dt()
return proxy->dt();
}
+// Replica exchange commands
+inline bool cvm::replica_enabled() {
+ return proxy->replica_enabled();
+}
+inline int cvm::replica_index() {
+ return proxy->replica_index();
+}
+inline int cvm::replica_num() {
+ return proxy->replica_num();
+}
+inline void cvm::replica_comm_barrier() {
+ return proxy->replica_comm_barrier();
+}
+inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
+ return proxy->replica_comm_recv(msg_data,buf_len,src_rep);
+}
+inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
+ return proxy->replica_comm_send(msg_data,msg_len,dest_rep);
+}
+
+
inline void cvm::request_system_force()
{
proxy->request_system_force (true);
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 89c229dde3..ec458d4867 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -5,7 +5,7 @@
#ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2014-09-19"
+#define COLVARPROXY_VERSION "2014-10-21"
#endif
@@ -68,8 +68,31 @@ public:
virtual int frame (int) { return COLVARS_NOT_IMPLEMENTED; }
- // **************** SIMULATION PARAMETERS ****************
+ // Replica exchange commands:
+
+ /// \brief Indicate if multi-replica support is available and active
+ virtual bool replica_enabled() { return false; }
+
+ /// \brief Index of this replica
+ virtual int replica_index() { return 0; }
+ /// \brief Total number of replica
+ virtual int replica_num() { return 1; }
+
+ /// \brief Synchronize replica
+ virtual void replica_comm_barrier() {}
+
+ /// \brief Receive data from other replica
+ virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
+ return COLVARS_NOT_IMPLEMENTED;
+ }
+
+ /// \brief Send data to other replica
+ virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
+ return COLVARS_NOT_IMPLEMENTED;
+ }
+
+ // **************** SIMULATION PARAMETERS ****************
/// \brief Prefix to be used for input files (restarts, not
/// configuration)
@@ -93,7 +116,7 @@ public:
return output_prefix_str;
}
- /// \brief Restarts will be fritten each time this number of steps has passed
+ /// \brief Restarts will be written each time this number of steps has passed
virtual size_t restart_frequency() = 0;
@@ -212,4 +235,3 @@ inline cvm::real colvarproxy::position_dist2 (cvm::atom_pos const &pos1,
}
#endif
-
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 5b6c333f7a..13dffd35d1 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -214,10 +214,7 @@ int colvarscript::proc_colvar (int argc, char const *argv[]) {
}
if (subcmd == "update") {
- // note: this is not sufficient for a newly created colvar
- // as atom coordinates may not be properly loaded
- // a full CVM update is required
- // or otherwise updating atom positions
+ cv->calc();
cv->update();
result = cvm::to_str(cv->value(), 0, cvm::cv_prec);
return COLVARSCRIPT_OK;
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 53dc1c1c28..b24be512a1 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -350,22 +350,22 @@ void colvarmodule::rotation::calc_optimal_rotation
// derivative of the S matrix
ds_1.reset();
- ds_1[0][0] = cvm::rvector ( a2x, a2y, a2z);
- ds_1[1][0] = cvm::rvector ( 0.0, a2z, -a2y);
+ ds_1[0][0].set( a2x, a2y, a2z);
+ ds_1[1][0].set( 0.0, a2z, -a2y);
ds_1[0][1] = ds_1[1][0];
- ds_1[2][0] = cvm::rvector (-a2z, 0.0, a2x);
+ ds_1[2][0].set(-a2z, 0.0, a2x);
ds_1[0][2] = ds_1[2][0];
- ds_1[3][0] = cvm::rvector ( a2y, -a2x, 0.0);
+ ds_1[3][0].set( a2y, -a2x, 0.0);
ds_1[0][3] = ds_1[3][0];
- ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z);
- ds_1[2][1] = cvm::rvector ( a2y, a2x, 0.0);
+ ds_1[1][1].set( a2x, -a2y, -a2z);
+ ds_1[2][1].set( a2y, a2x, 0.0);
ds_1[1][2] = ds_1[2][1];
- ds_1[3][1] = cvm::rvector ( a2z, 0.0, a2x);
+ ds_1[3][1].set( a2z, 0.0, a2x);
ds_1[1][3] = ds_1[3][1];
- ds_1[2][2] = cvm::rvector (-a2x, a2y, -a2z);
- ds_1[3][2] = cvm::rvector ( 0.0, a2z, a2y);
+ ds_1[2][2].set(-a2x, a2y, -a2z);
+ ds_1[3][2].set( 0.0, a2z, a2y);
ds_1[2][3] = ds_1[3][2];
- ds_1[3][3] = cvm::rvector (-a2x, -a2y, a2z);
+ ds_1[3][3].set(-a2x, -a2y, a2z);
cvm::rvector &dl0_1 = dL0_1[ia];
vector1d<cvm::rvector, 4> &dq0_1 = dQ0_1[ia];
@@ -404,22 +404,22 @@ void colvarmodule::rotation::calc_optimal_rotation
matrix2d<cvm::rvector, 4, 4> &ds_2 = dS_2[ia];
ds_2.reset();
- ds_2[0][0] = cvm::rvector ( a1x, a1y, a1z);
- ds_2[1][0] = cvm::rvector ( 0.0, -a1z, a1y);
+ ds_2[0][0].set( a1x, a1y, a1z);
+ ds_2[1][0].set( 0.0, -a1z, a1y);
ds_2[0][1] = ds_2[1][0];
- ds_2[2][0] = cvm::rvector ( a1z, 0.0, -a1x);
+ ds_2[2][0].set( a1z, 0.0, -a1x);
ds_2[0][2] = ds_2[2][0];
- ds_2[3][0] = cvm::rvector (-a1y, a1x, 0.0);
+ ds_2[3][0].set(-a1y, a1x, 0.0);
ds_2[0][3] = ds_2[3][0];
- ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z);
- ds_2[2][1] = cvm::rvector ( a1y, a1x, 0.0);
+ ds_2[1][1].set( a1x, -a1y, -a1z);
+ ds_2[2][1].set( a1y, a1x, 0.0);
ds_2[1][2] = ds_2[2][1];
- ds_2[3][1] = cvm::rvector ( a1z, 0.0, a1x);
+ ds_2[3][1].set( a1z, 0.0, a1x);
ds_2[1][3] = ds_2[3][1];
- ds_2[2][2] = cvm::rvector (-a1x, a1y, -a1z);
- ds_2[3][2] = cvm::rvector ( 0.0, a1z, a1y);
+ ds_2[2][2].set(-a1x, a1y, -a1z);
+ ds_2[3][2].set( 0.0, a1z, a1y);
ds_2[2][3] = ds_2[3][2];
- ds_2[3][3] = cvm::rvector (-a1x, -a1y, a1z);
+ ds_2[3][3].set(-a1x, -a1y, a1z);
cvm::rvector &dl0_2 = dL0_2[ia];
vector1d<cvm::rvector, 4> &dq0_2 = dQ0_2[ia];
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 0d97740622..18a6de3cfb 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -37,10 +37,19 @@ public:
{}
/// \brief Set all components to a scalar value
- inline void set (cvm::real const &value = 0.0) {
+ inline void set (cvm::real const &value) {
x = y = z = value;
}
+ /// \brief Assign all components
+ inline void set (cvm::real const &x_i,
+ cvm::real const &y_i,
+ cvm::real const &z_i) {
+ x = x_i;
+ y = y_i;
+ z = z_i;
+ }
+
/// \brief Set all components to zero
inline void reset() {
x = y = z = 0.0;
--
GitLab