diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index daa1393269b953d569aff9846a135480dcdd42c7..2d28758819a9a7a0dac700cd4a8239f51257b8e1 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index ce76b3b9eb2c8bd5f37b85af332611fa7ca963b3..5a4e8b617a9007fff02e87233b8f8dfebb1e7a25 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -195,7 +195,7 @@ int colvar::init(std::string const &conf)
// - it is homogeneous
// - all cvcs are periodic
// - all cvcs have the same period
- if (cvcs[0]->b_periodic) { // TODO make this a CVC feature
+ if (is_enabled(f_cv_homogeneous) && 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++) {
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 20dad2771b6c25afef9c4667ca23f24cfe61d46d..32c329460d06cfa5ceba43e245aa08a0d0810100 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -11,8 +11,6 @@
#define COLVAR_H
#include <iostream>
-#include <iomanip>
-#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
@@ -60,10 +58,13 @@ public:
/// \brief Current actual value (not extended DOF)
colvarvalue const & actual_value() const;
-
+
+ /// \brief Current running average (if calculated as set by analysis flag)
+ colvarvalue const & run_ave() const;
+
/// \brief Force constant of the spring
cvm::real const & force_constant() const;
-
+
/// \brief Current velocity (previously set by calc() or by read_traj())
colvarvalue const & velocity() const;
@@ -516,7 +517,7 @@ public:
// collective variable component base class
class cvc;
- // currently available collective variable components
+ // list of available collective variable components
// scalar colvar components
class distance;
@@ -611,12 +612,15 @@ inline colvarvalue const & colvar::value() const
return x_reported;
}
-
inline colvarvalue const & colvar::actual_value() const
{
return x;
}
+inline colvarvalue const & colvar::run_ave() const
+{
+ return runave;
+}
inline colvarvalue const & colvar::velocity() const
{
diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h
index 7fc7f870a10932d176525bb3c6c42d574a572fb1..36ed938119d41232dfdbdc42589cd35f8f67a735 100644
--- a/lib/colvars/colvar_UIestimator.h
+++ b/lib/colvars/colvar_UIestimator.h
@@ -45,7 +45,7 @@ namespace UIestimator {
this->width = width;
this->dimension = lowerboundary.size();
this->y_size = y_size; // keep in mind the internal (spare) matrix is stored in diagonal form
- this->y_total_size = int(pow(double(y_size), dimension) + EPSILON);
+ this->y_total_size = int(std::pow(double(y_size), double(dimension)) + EPSILON);
// the range of the matrix is [lowerboundary, upperboundary]
x_total_size = 1;
@@ -121,7 +121,7 @@ namespace UIestimator {
int index = 0;
for (i = 0; i < dimension; i++) {
if (i + 1 < dimension)
- index += temp[i] * int(pow(double(y_size), dimension - i - 1) + EPSILON);
+ index += temp[i] * int(std::pow(double(y_size), double(dimension - i - 1)) + EPSILON);
else
index += temp[i];
}
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index d2a0f0a807d2727500d96e6a19a56eaa0af23f39..1be6f42e594f86adf16a5b4702a631a6bef9c1e7 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -8,9 +8,11 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvarparse.h"
#include "colvaratoms.h"
+
cvm::atom::atom()
{
index = -1;
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 71c587e23084e516c6b66358e97327c8d404822f..0dda6ab792d8ae20afcf6a9aeaa615675f332d1b 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -11,6 +11,7 @@
#define COLVARATOMS_H
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvarparse.h"
#include "colvardeps.h"
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 301e83e73015a91f14ff75a82c045954083d4074..29620fbee8f103a032b5698c52ca841b2f122cbc 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -8,6 +8,7 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvarvalue.h"
#include "colvarbias.h"
#include "colvargrid.h"
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index e4aea8eb86e7dcb6712fc6dc2c771adeeb6909f6..b3b5b3eb1654df12992b2fafb206d82851d3788b 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -8,6 +8,7 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias_abf.h"
@@ -18,16 +19,18 @@ colvarbias_abf::colvarbias_abf(char const *key)
b_CZAR_estimator(false),
system_force(NULL),
gradients(NULL),
+ pmf(NULL),
samples(NULL),
z_gradients(NULL),
z_samples(NULL),
czar_gradients(NULL),
+ czar_pmf(NULL),
last_gradients(NULL),
- last_samples(NULL)
+ last_samples(NULL),
+ pabf_freq(0)
{
}
-
int colvarbias_abf::init(std::string const &conf)
{
colvarbias::init(conf);
@@ -91,7 +94,7 @@ int colvarbias_abf::init(std::string const &conf)
// ************* checking the associated colvars *******************
- if (colvars.size() == 0) {
+ if (num_variables() == 0) {
cvm::error("Error: no collective variables specified for the ABF bias.\n");
return COLVARS_ERROR;
}
@@ -102,7 +105,8 @@ int colvarbias_abf::init(std::string const &conf)
}
bool b_extended = false;
- for (size_t i = 0; i < colvars.size(); i++) {
+ size_t i;
+ for (i = 0; i < num_variables(); i++) {
if (colvars[i]->value().type() != colvarvalue::type_scalar) {
cvm::error("Error: ABF bias can only use scalar-type variables.\n");
@@ -132,10 +136,10 @@ int colvarbias_abf::init(std::string const &conf)
}
if (get_keyval(conf, "maxForce", max_force)) {
- if (max_force.size() != colvars.size()) {
+ if (max_force.size() != num_variables()) {
cvm::error("Error: Number of parameters to maxForce does not match number of colvars.");
}
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
if (max_force[i] < 0.0) {
cvm::error("Error: maxForce should be non-negative.");
}
@@ -145,9 +149,9 @@ int colvarbias_abf::init(std::string const &conf)
cap_force = false;
}
- bin.assign(colvars.size(), 0);
- force_bin.assign(colvars.size(), 0);
- system_force = new cvm::real [colvars.size()];
+ bin.assign(num_variables(), 0);
+ force_bin.assign(num_variables(), 0);
+ system_force = new cvm::real [num_variables()];
// Construct empty grids based on the colvars
if (cvm::debug()) {
@@ -159,14 +163,14 @@ int colvarbias_abf::init(std::string const &conf)
gradients->samples = samples;
samples->has_parent_data = true;
- // Data for eABF z-based estimator
- if (b_extended) {
+ // Data for eAB F z-based estimator
+ if ( b_extended ) {
get_keyval(conf, "CZARestimator", b_CZAR_estimator, true);
// CZAR output files for stratified eABF
get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
colvarparse::parse_silent);
- z_bin.assign(colvars.size(), 0);
+ z_bin.assign(num_variables(), 0);
z_samples = new colvar_grid_count(colvars);
z_samples->request_actual_value();
z_gradients = new colvar_grid_gradient(colvars);
@@ -176,6 +180,27 @@ int colvarbias_abf::init(std::string const &conf)
czar_gradients = new colvar_grid_gradient(colvars);
}
+ // For now, we integrate on-the-fly iff the grid is < 3D
+ if ( num_variables() <= 3 ) {
+ pmf = new integrate_potential(colvars, gradients);
+ if ( b_CZAR_estimator ) {
+ czar_pmf = new integrate_potential(colvars, czar_gradients);
+ }
+ get_keyval(conf, "integrate", b_integrate, true); // Integrate for output
+ if ( num_variables() > 1 ) {
+ // Projected ABF
+ get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0);
+ // Parameters for integrating initial (and final) gradient data
+ get_keyval(conf, "integrateInitSteps", integrate_initial_steps, 1e4);
+ get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6);
+ // for updating the integrated PMF on the fly
+ get_keyval(conf, "integrateSteps", integrate_steps, 100);
+ get_keyval(conf, "integrateTol", integrate_tol, 1e-4);
+ }
+ } else {
+ b_integrate = false;
+ }
+
// 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).
@@ -188,6 +213,8 @@ int colvarbias_abf::init(std::string const &conf)
// If custom grids are provided, read them
if ( input_prefix.size() > 0 ) {
read_gradients_samples();
+ // Update divergence to account for input data
+ pmf->set_div();
}
// if extendedLangrangian is on, then call UI estimator
@@ -202,7 +229,7 @@ int colvarbias_abf::init(std::string const &conf)
bool UI_restart = (input_prefix.size() > 0);
- for (size_t i = 0; i < colvars.size(); i++)
+ for (i = 0; i < num_variables(); i++)
{
UI_lowerboundary.push_back(colvars[i]->lower_boundary);
UI_upperboundary.push_back(colvars[i]->upper_boundary);
@@ -238,6 +265,11 @@ colvarbias_abf::~colvarbias_abf()
gradients = NULL;
}
+ if (pmf) {
+ delete pmf;
+ pmf = NULL;
+ }
+
if (z_samples) {
delete z_samples;
z_samples = NULL;
@@ -253,6 +285,11 @@ colvarbias_abf::~colvarbias_abf()
czar_gradients = NULL;
}
+ if (czar_pmf) {
+ delete czar_pmf;
+ czar_pmf = NULL;
+ }
+
// shared ABF
// We used to only do this if "shared" was defined,
// but now we can call shared externally
@@ -278,44 +315,48 @@ colvarbias_abf::~colvarbias_abf()
int colvarbias_abf::update()
{
- if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
+ int iter;
- if (cvm::step_relative() == 0) {
-
- // At first timestep, do only:
- // initialization stuff (file operations relying on n_abf_biases
- // compute current value of colvars
-
- for (size_t i = 0; i < colvars.size(); i++) {
- bin[i] = samples->current_bin_scalar(i);
- }
+ if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
- } else {
+ size_t i;
+ for (i = 0; i < num_variables(); i++) {
+ bin[i] = samples->current_bin_scalar(i);
+ }
+ if (cvm::proxy->total_forces_same_step()) {
+ // e.g. in LAMMPS, total forces are current
+ force_bin = bin;
+ }
- for (size_t i = 0; i < colvars.size(); i++) {
- bin[i] = samples->current_bin_scalar(i);
- }
+ if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
- if ( update_bias && samples->index_ok(force_bin) ) {
- // Only if requested and within bounds of the grid...
+ if (update_bias) {
+// if (b_adiabatic_reweighting) {
+// // Update gradients non-locally based on conditional distribution of
+// // fictitious variable TODO
+//
+// } else
+ if (samples->index_ok(force_bin)) {
+ // Only if requested and within bounds of the grid...
- for (size_t i = 0; i < colvars.size(); i++) {
- // get total forces (lagging by 1 timestep) from colvars
- // and subtract previous ABF force if necessary
- update_system_force(i);
- }
- if (cvm::proxy->total_forces_same_step()) {
- // e.g. in LAMMPS, total forces are current
- force_bin = bin;
+ for (i = 0; i < num_variables(); i++) {
+ // get total forces (lagging by 1 timestep) from colvars
+ // and subtract previous ABF force if necessary
+ update_system_force(i);
+ }
+ gradients->acc_force(force_bin, system_force);
+ if ( b_integrate ) {
+ pmf->update_div_neighbors(force_bin);
+ }
}
- gradients->acc_force(force_bin, system_force);
}
+
if ( z_gradients && update_bias ) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
z_bin[i] = z_samples->current_bin_scalar(i);
}
if ( z_samples->index_ok(z_bin) ) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
// If we are outside the range of xi, the force has not been obtained above
// the function is just an accessor, so cheap to call again anyway
update_system_force(i);
@@ -323,6 +364,14 @@ int colvarbias_abf::update()
z_gradients->acc_force(z_bin, system_force);
}
}
+
+ if ( b_integrate ) {
+ if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
+ cvm::real err;
+ iter = pmf->integrate(integrate_steps, integrate_tol, err);
+ pmf->set_zero_minimum(); // TODO: do this only when necessary
+ }
+ }
}
if (!cvm::proxy->total_forces_same_step()) {
@@ -332,14 +381,14 @@ int colvarbias_abf::update()
}
// Reset biasing forces from previous timestep
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvar_forces[i].reset();
}
// Compute and apply the new bias, if applicable
if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
- size_t count = samples->value(bin);
+ cvm::real count = samples->value(bin);
cvm::real fact = 1.0;
// Factor that ensures smooth introduction of the force
@@ -348,21 +397,34 @@ int colvarbias_abf::update()
(cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
}
- const cvm::real * grad = &(gradients->value(bin));
+ std::vector<cvm::real> grad(num_variables());
+ if ( pabf_freq ) {
+ // In projected ABF, the force is the PMF gradient estimate
+ pmf->vector_gradient_finite_diff(bin, grad);
+ } else {
+ // Normal ABF
+ gradients->vector_value(bin, grad);
+ }
+
+// if ( b_adiabatic_reweighting) {
+// // Average of force according to conditional distribution of fictitious variable
+// // need freshly integrated PMF, gradient TODO
+// } else
if ( fact != 0.0 ) {
- if ( (colvars.size() == 1) && colvars[0]->periodic_boundaries() ) {
+ if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) {
// Enforce a zero-mean bias on periodic, 1D coordinates
// in other words: boundary condition is that the biasing potential is periodic
- colvar_forces[0].real_value = fact * (grad[0] / cvm::real(count) - gradients->average());
+ // This is enforced naturally if using integrated PMF
+ colvar_forces[0].real_value = fact * (grad[0] - gradients->average ());
} else {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
// subtracting the mean force (opposite of the FE gradient) means adding the gradient
- colvar_forces[i].real_value = fact * grad[i] / cvm::real(count);
+ colvar_forces[i].real_value = fact * grad[i];
}
}
if (cap_force) {
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) {
colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]);
}
@@ -407,9 +469,9 @@ int colvarbias_abf::update()
// update UI estimator every step
if (b_UI_estimator)
{
- std::vector<double> x(colvars.size(),0);
- std::vector<double> y(colvars.size(),0);
- for (size_t i = 0; i < colvars.size(); i++)
+ std::vector<double> x(num_variables(),0);
+ std::vector<double> y(num_variables(),0);
+ for (size_t i = 0; i < num_variables(); i++)
{
x[i] = colvars[i]->actual_value();
y[i] = colvars[i]->value();
@@ -509,26 +571,60 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) {
cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
+ return;
}
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
+ // In dimension higher than 2, dx is easier to handle and visualize
+ if (num_variables() > 2) {
+ std::string samples_dx_out_name = prefix + ".count.dx";
+ std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
+ if (!samples_os) {
+ cvm::error("Error opening samples file " + samples_dx_out_name + " for writing");
+ return;
+ }
+ samples->write_opendx(*samples_dx_os);
+ *samples_dx_os << std::endl;
+ cvm::proxy->close_output_stream(samples_dx_out_name);
+ }
+
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) {
cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
+ return;
}
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
- if (colvars.size() == 1) {
- // Do numerical integration and output a PMF
+ if (b_integrate) {
+ // Do numerical integration (to high precision) and output a PMF
+ cvm::real err;
+ pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
+ pmf->set_zero_minimum();
+
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) {
cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
+ return;
+ }
+ pmf->write_multicol(*pmf_os);
+
+ // In dimension higher than 2, dx is easier to handle and visualize
+ if (num_variables() > 2) {
+ std::string pmf_dx_out_name = prefix + ".pmf.dx";
+ std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
+ if (!pmf_dx_os) {
+ cvm::error("Error opening pmf file " + pmf_dx_out_name + " for writing");
+ return;
+ }
+ pmf->write_opendx(*pmf_dx_os);
+ *pmf_dx_os << std::endl;
+ cvm::proxy->close_output_stream(pmf_dx_out_name);
}
- gradients->write_1D_integral(*pmf_os);
+
*pmf_os << std::endl;
cvm::proxy->close_output_stream(pmf_out_name);
}
@@ -542,6 +638,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) {
cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
+ return;
}
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
@@ -553,6 +650,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) {
cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
+ return;
}
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
@@ -563,8 +661,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
- - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n),
- n);
+ - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n);
}
}
@@ -574,17 +671,39 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) {
cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
+ return;
}
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
- if (colvars.size() == 1) {
- // Do numerical integration and output a PMF
+ if (b_integrate) {
+ // Do numerical integration (to high precision) and output a PMF
+ cvm::real err;
+ czar_pmf->set_div();
+ czar_pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
+ czar_pmf->set_zero_minimum();
+
std::string czar_pmf_out_name = prefix + ".czar.pmf";
- std::ostream *czar_pmf_os =
- cvm::proxy->output_stream(czar_pmf_out_name, mode);
- if (!czar_pmf_os) cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
- czar_gradients->write_1D_integral(*czar_pmf_os);
+ std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
+ if (!czar_pmf_os) {
+ cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
+ return;
+ }
+ czar_pmf->write_multicol(*czar_pmf_os);
+
+ // In dimension higher than 2, dx is easier to handle and visualize
+ if (num_variables() > 2) {
+ std::string czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
+ std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
+ if (!czar_pmf_dx_os) {
+ cvm::error("Error opening CZAR pmf file " + czar_pmf_dx_out_name + " for writing");
+ return;
+ }
+ czar_pmf->write_opendx(*czar_pmf_dx_os);
+ *czar_pmf_dx_os << std::endl;
+ cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
+ }
+
*czar_pmf_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_out_name);
}
@@ -708,6 +827,10 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
if (! gradients->read_raw(is)) {
return is;
}
+ if (b_integrate) {
+ // Update divergence to account for restart data
+ pmf->set_div();
+ }
if (b_CZAR_estimator) {
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index 1defe72268c487e82a536d2f5e441e1067897272..0260401292ecc40ee38c5fce3c282f359603591f 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -40,28 +40,44 @@ private:
/// Base filename(s) for reading previous gradient data (replaces data from restart file)
std::vector<std::string> input_prefix;
- bool update_bias;
- bool hide_Jacobian;
- size_t full_samples;
- size_t min_samples;
+ bool update_bias;
+ bool hide_Jacobian;
+ bool b_integrate;
+
+ size_t full_samples;
+ size_t min_samples;
/// frequency for updating output files
- int output_freq;
+ int output_freq;
/// Write combined files with a history of all output data?
- bool b_history_files;
+ bool b_history_files;
/// Write CZAR output file for stratified eABF (.zgrad)
- bool b_czar_window_file;
- size_t history_freq;
+ bool b_czar_window_file;
+ size_t history_freq;
/// Umbrella Integration estimator of free energy from eABF
UIestimator::UIestimator eabf_UI;
- // Run UI estimator?
- bool b_UI_estimator;
- // Run CZAR estimator?
- bool b_CZAR_estimator;
-
- /// Cap applied biasing force?
+ /// Run UI estimator?
+ bool b_UI_estimator;
+ /// Run CZAR estimator?
+ bool b_CZAR_estimator;
+
+ /// Frequency for updating pABF PMF (if zero, pABF is not used)
+ int pabf_freq;
+ /// Max number of CG iterations for integrating PMF at startup and for file output
+ int integrate_initial_steps;
+ /// Tolerance for integrating PMF at startup and for file output
+ cvm::real integrate_initial_tol;
+ /// Max number of CG iterations for integrating PMF at on-the-fly pABF updates
+ int integrate_steps;
+ /// Tolerance for integrating PMF at on-the-fly pABF updates
+ cvm::real integrate_tol;
+
+ /// Cap the biasing force to be applied?
bool cap_force;
std::vector<cvm::real> max_force;
+ // Frequency for updating 2D gradients
+ int integrate_freq;
+
// Internal data and methods
std::vector<int> bin, force_bin, z_bin;
@@ -71,12 +87,16 @@ private:
colvar_grid_gradient *gradients;
/// n-dim grid of number of samples
colvar_grid_count *samples;
+ /// n-dim grid of pmf (dimension 1 to 3)
+ integrate_potential *pmf;
/// n-dim grid: average force on "real" coordinate for eABF z-based estimator
colvar_grid_gradient *z_gradients;
/// n-dim grid of number of samples on "real" coordinate for eABF z-based estimator
colvar_grid_count *z_samples;
/// n-dim grid contining CZAR estimator of "real" free energy gradients
colvar_grid_gradient *czar_gradients;
+ /// n-dim grid of CZAR pmf (dimension 1 to 3)
+ integrate_potential *czar_pmf;
inline int update_system_force(size_t i)
{
@@ -96,9 +116,9 @@ private:
}
// shared ABF
- bool shared_on;
- size_t shared_freq;
- int shared_last_step;
+ bool shared_on;
+ size_t shared_freq;
+ int shared_last_step;
// Share between replicas -- may be called independently of update
virtual int replica_share();
@@ -114,12 +134,12 @@ private:
//// 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);
+ /// Write human-readable FE gradients and sample count, and DX file in dim > 2
+ 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();
+ void read_gradients_samples();
std::istream& read_state_data(std::istream&);
std::ostream& write_state_data(std::ostream&);
diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp
index 124a15c5da3f7b59e1e1da4ef5d1d65dc0f60c67..187ecc363a9da559555b7eb296703167f1d26389 100644
--- a/lib/colvars/colvarbias_alb.cpp
+++ b/lib/colvars/colvarbias_alb.cpp
@@ -7,13 +7,11 @@
// 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>
+#include <cstdlib>
#include "colvarmodule.h"
-#include "colvarbias_alb.h"
#include "colvarbias.h"
+#include "colvarbias_alb.h"
#ifdef _MSC_VER
#if _MSC_VER <= 1700
@@ -45,22 +43,22 @@ int colvarbias_alb::init(std::string const &conf)
size_t i;
// get the initial restraint centers
- colvar_centers.resize(colvars.size());
+ colvar_centers.resize(num_variables());
- means.resize(colvars.size());
- ssd.resize(colvars.size()); //sum of squares of differences from mean
+ means.resize(num_variables());
+ ssd.resize(num_variables()); //sum of squares of differences from mean
//setup force vectors
- max_coupling_range.resize(colvars.size());
- max_coupling_rate.resize(colvars.size());
- coupling_accum.resize(colvars.size());
- set_coupling.resize(colvars.size());
- current_coupling.resize(colvars.size());
- coupling_rate.resize(colvars.size());
+ max_coupling_range.resize(num_variables());
+ max_coupling_rate.resize(num_variables());
+ coupling_accum.resize(num_variables());
+ set_coupling.resize(num_variables());
+ current_coupling.resize(num_variables());
+ coupling_rate.resize(num_variables());
enable(f_cvb_apply_force);
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvar_centers[i].type(colvars[i]->value());
//zero moments
means[i] = ssd[i] = 0;
@@ -70,7 +68,7 @@ int colvarbias_alb::init(std::string const &conf)
}
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvar_centers[i].apply_constraints();
}
} else {
@@ -78,7 +76,7 @@ int colvarbias_alb::init(std::string const &conf)
cvm::fatal_error("Error: must define the initial centers of adaptive linear bias .\n");
}
- if (colvar_centers.size() != colvars.size())
+ if (colvar_centers.size() != num_variables())
cvm::fatal_error("Error: number of centers does not match "
"that of collective variables.\n");
@@ -100,17 +98,17 @@ int colvarbias_alb::init(std::string const &conf)
//initial guess
if (!get_keyval(conf, "forceConstant", set_coupling, set_coupling))
- for (i =0 ; i < colvars.size(); i++)
+ for (i =0 ; i < num_variables(); i++)
set_coupling[i] = 0.;
//how we're going to increase to that point
- for (i = 0; i < colvars.size(); i++)
+ for (i = 0; i < num_variables(); i++)
coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq;
if (!get_keyval(conf, "forceRange", max_coupling_range, max_coupling_range)) {
//set to default
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
if (cvm::temperature() > 0)
max_coupling_range[i] = 3 * cvm::temperature() * cvm::boltzmann();
else
@@ -120,7 +118,7 @@ int colvarbias_alb::init(std::string const &conf)
if (!get_keyval(conf, "rateMax", max_coupling_rate, max_coupling_rate)) {
//set to default
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
max_coupling_rate[i] = max_coupling_range[i] / (10 * update_freq);
}
}
@@ -151,7 +149,7 @@ int colvarbias_alb::update()
// Force and energy calculation
bool finished_equil_flag = 1;
cvm::real delta;
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(current_coupling[i], colvars[i]->width),
colvars[i],
colvar_centers[i]);
@@ -168,7 +166,9 @@ int colvarbias_alb::update()
} else {
//check if we've reached the setpoint
- if (coupling_rate[i] == 0 || pow(current_coupling[i] - set_coupling[i],2) < pow(coupling_rate[i],2)) {
+ cvm::real const coupling_diff = current_coupling[i] - set_coupling[i];
+ if ((coupling_rate[i] == 0) ||
+ ((coupling_diff*coupling_diff) < (coupling_rate[i]*coupling_rate[i]))) {
finished_equil_flag &= 1; //we continue equilibrating as long as we haven't reached all the set points
}
else {
@@ -209,7 +209,7 @@ int colvarbias_alb::update()
cvm::real temp;
//reset means and sum of squares of differences
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
temp = 2. * (means[i] / (static_cast<cvm::real> (colvar_centers[i])) - 1) * ssd[i] / (update_calls - 1);
@@ -222,7 +222,7 @@ int colvarbias_alb::update()
ssd[i] = 0;
//stochastic if we do that update or not
- if (colvars.size() == 1 || rand() < RAND_MAX / ((int) colvars.size())) {
+ if (num_variables() == 1 || rand() < RAND_MAX / ((int) num_variables())) {
coupling_accum[i] += step_size * step_size;
current_coupling[i] = set_coupling[i];
set_coupling[i] += max_coupling_range[i] / sqrt(coupling_accum[i]) * step_size;
@@ -284,37 +284,37 @@ std::string const colvarbias_alb::get_state_params() const
std::ostringstream os;
os << " setCoupling ";
size_t i;
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << set_coupling[i] << "\n";
}
os << " currentCoupling ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << current_coupling[i] << "\n";
}
os << " maxCouplingRange ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << max_coupling_range[i] << "\n";
}
os << " couplingRate ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << coupling_rate[i] << "\n";
}
os << " couplingAccum ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << coupling_accum[i] << "\n";
}
os << " mean ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << means[i] << "\n";
}
os << " ssd ";
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
os << std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << ssd[i] << "\n";
}
@@ -350,7 +350,7 @@ std::ostream & colvarbias_alb::write_traj_label(std::ostream &os)
}
if (b_output_centers)
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
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);
@@ -378,7 +378,7 @@ std::ostream & colvarbias_alb::write_traj(std::ostream &os)
if (b_output_centers)
- for (size_t i = 0; i < colvars.size(); i++) {
+ for (size_t i = 0; i < num_variables(); i++) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers[i];
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index 0722e6384dc7047e5f6eb0bdf0c623537cfaeb53..329b1d9dc08d918295f5ec89a0519a4a292168a0 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -8,10 +8,10 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias_histogram.h"
-/// Histogram "bias" constructor
colvarbias_histogram::colvarbias_histogram(char const *key)
: colvarbias(key),
@@ -44,7 +44,7 @@ int colvarbias_histogram::init(std::string const &conf)
get_keyval(conf, "gatherVectorColvars", colvar_array, colvar_array);
if (colvar_array) {
- for (i = 0; i < colvars.size(); i++) { // should be all vector
+ for (i = 0; i < num_variables(); i++) { // should be all vector
if (colvars[i]->value().type() != colvarvalue::type_vector) {
cvm::error("Error: used gatherVectorColvars with non-vector colvar.\n", INPUT_ERROR);
return INPUT_ERROR;
@@ -63,7 +63,7 @@ int colvarbias_histogram::init(std::string const &conf)
}
}
} else {
- for (i = 0; i < colvars.size(); i++) { // should be all scalar
+ for (i = 0; i < num_variables(); i++) { // should be all scalar
if (colvars[i]->value().type() != colvarvalue::type_scalar) {
cvm::error("Error: only scalar colvars are supported when gatherVectorColvars is off.\n", INPUT_ERROR);
return INPUT_ERROR;
@@ -77,7 +77,7 @@ int colvarbias_histogram::init(std::string const &conf)
get_keyval(conf, "weights", weights, weights);
}
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
colvars[i]->enable(f_cv_grid);
}
@@ -116,7 +116,7 @@ int colvarbias_histogram::update()
}
// assign a valid bin size
- bin.assign(colvars.size(), 0);
+ bin.assign(num_variables(), 0);
if (out_name.size() == 0) {
// At the first timestep, we need to assign out_name since
@@ -137,7 +137,7 @@ int colvarbias_histogram::update()
if (colvar_array_size == 0) {
// update indices for scalar values
size_t i;
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
bin[i] = grid->current_bin_scalar(i);
}
@@ -148,7 +148,7 @@ int colvarbias_histogram::update()
// update indices for vector/array values
size_t iv, i;
for (iv = 0; iv < colvar_array_size; iv++) {
- for (i = 0; i < colvars.size(); i++) {
+ for (i = 0; i < num_variables(); i++) {
bin[i] = grid->current_bin_scalar(i, iv);
}
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index b0d154dfc9ebf591f49b15ba517f4bcf7ba10317..f3ae3631a02f2478cb62a7f9a84e7535926501c4 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -27,7 +27,8 @@
#define PATHSEP "/"
#endif
-
+#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias_meta.h"
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index 23534f56eb43c1cd0a2e257ed8581449026d0c8f..4ed1a95f942fe664d336b7dfb170f1ea1c3a62a8 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -7,7 +7,10 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
+#include <cmath>
+
#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvarvalue.h"
#include "colvarbias_restraint.h"
@@ -150,13 +153,14 @@ colvarbias_restraint_k::colvarbias_restraint_k(char const *key)
: colvarbias(key), colvarbias_ti(key), colvarbias_restraint(key)
{
force_k = -1.0;
+ check_positive_k = true;
}
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) {
+ if (check_positive_k && (force_k < 0.0)) {
cvm::error("Error: undefined or invalid force constant.\n", INPUT_ERROR);
return INPUT_ERROR;
}
@@ -177,6 +181,7 @@ colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key)
target_nstages = 0;
target_nsteps = 0;
stage = 0;
+ acc_work = 0.0;
b_chg_centers = false;
b_chg_force_k = false;
}
@@ -203,6 +208,14 @@ int colvarbias_restraint_moving::init(std::string const &conf)
cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", INPUT_ERROR);
return cvm::get_error();
}
+
+ get_keyval_feature(this, conf, "outputAccumulatedWork",
+ f_cvb_output_acc_work,
+ is_enabled(f_cvb_output_acc_work));
+ if (is_enabled(f_cvb_output_acc_work) && (target_nstages > 0)) {
+ return cvm::error("Error: outputAccumulatedWork and targetNumStages "
+ "are incompatible.\n", INPUT_ERROR);
+ }
}
return COLVARS_OK;
@@ -246,8 +259,6 @@ colvarbias_restraint_centers_moving::colvarbias_restraint_centers_moving(char co
{
b_chg_centers = false;
b_output_centers = false;
- b_output_acc_work = false;
- acc_work = 0.0;
}
@@ -288,9 +299,6 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf)
0.5);
}
- get_keyval(conf, "outputAccumulatedWork", b_output_acc_work,
- b_output_acc_work); // TODO this conflicts with stages
-
} else {
target_centers.clear();
}
@@ -382,12 +390,14 @@ int colvarbias_restraint_centers_moving::update()
int colvarbias_restraint_centers_moving::update_acc_work()
{
- if (b_output_acc_work) {
- if ((cvm::step_relative() > 0) &&
- (cvm::step_absolute() <= target_nsteps)) {
- for (size_t i = 0; i < num_variables(); i++) {
- // project forces on the calculated increments at this step
- acc_work += colvar_forces[i] * centers_incr[i];
+ if (b_chg_centers) {
+ if (is_enabled(f_cvb_output_acc_work)) {
+ if ((cvm::step_relative() > 0) &&
+ (cvm::step_absolute() <= target_nsteps)) {
+ for (size_t i = 0; i < num_variables(); i++) {
+ // project forces on the calculated increments at this step
+ acc_work += colvar_forces[i] * centers_incr[i];
+ }
}
}
}
@@ -410,7 +420,7 @@ std::string const colvarbias_restraint_centers_moving::get_state_params() const
}
os << "\n";
- if (b_output_acc_work) {
+ if (is_enabled(f_cvb_output_acc_work)) {
os << "accumulatedWork "
<< std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
<< acc_work << "\n";
@@ -429,7 +439,7 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con
// 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 (b_output_acc_work) {
+ if (is_enabled(f_cvb_output_acc_work)) {
if (!get_keyval(conf, "accumulatedWork", acc_work))
cvm::error("Error: accumulatedWork is missing from the restart.\n");
}
@@ -449,7 +459,7 @@ std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostrea
}
}
- if (b_output_acc_work) {
+ if (b_chg_centers && is_enabled(f_cvb_output_acc_work)) {
os << " W_"
<< cvm::wrap_string(this->name, cvm::en_width-2);
}
@@ -468,7 +478,7 @@ std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
}
}
- if (b_output_acc_work) {
+ if (b_chg_centers && is_enabled(f_cvb_output_acc_work)) {
os << " "
<< std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
<< acc_work;
@@ -488,10 +498,11 @@ colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key)
{
b_chg_force_k = false;
target_equil_steps = 0;
- target_force_k = 0.0;
- starting_force_k = 0.0;
+ target_force_k = -1.0;
+ starting_force_k = -1.0;
force_k_exp = 1.0;
restraint_FE = 0.0;
+ force_k_incr = 0.0;
}
@@ -569,14 +580,13 @@ int colvarbias_restraint_k_moving::update()
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;
+ // Derivative of energy with respect to force_k
+ cvm::real dU_dk = 0.0;
for (size_t i = 0; i < num_variables(); i++) {
- dist_sq += d_restraint_potential_dk(i);
+ dU_dk += 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;
+ restraint_FE += force_k_exp * std::pow(lambda, force_k_exp - 1.0)
+ * (target_force_k - starting_force_k) * dU_dk;
}
// Finish current stage...
@@ -607,10 +617,13 @@ int colvarbias_restraint_k_moving::update()
} else if (cvm::step_absolute() <= target_nsteps) {
+
// update force constant (slow growth)
lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
+ cvm::real const force_k_old = force_k;
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow(lambda, force_k_exp);
+ force_k_incr = force_k - force_k_old;
}
}
@@ -618,6 +631,23 @@ int colvarbias_restraint_k_moving::update()
}
+int colvarbias_restraint_k_moving::update_acc_work()
+{
+ if (b_chg_force_k) {
+ if (is_enabled(f_cvb_output_acc_work)) {
+ if (cvm::step_relative() > 0) {
+ cvm::real dU_dk = 0.0;
+ for (size_t i = 0; i < num_variables(); i++) {
+ dU_dk += d_restraint_potential_dk(i);
+ }
+ acc_work += dU_dk * force_k_incr;
+ }
+ }
+ }
+ return COLVARS_OK;
+}
+
+
std::string const colvarbias_restraint_k_moving::get_state_params() const
{
std::ostringstream os;
@@ -626,6 +656,12 @@ std::string const colvarbias_restraint_k_moving::get_state_params() const
os << "forceConstant "
<< std::setprecision(cvm::en_prec)
<< std::setw(cvm::en_width) << force_k << "\n";
+
+ if (is_enabled(f_cvb_output_acc_work)) {
+ os << "accumulatedWork "
+ << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+ << acc_work << "\n";
+ }
}
return os.str();
}
@@ -639,6 +675,10 @@ int colvarbias_restraint_k_moving::set_state_params(std::string const &conf)
// 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");
+ if (is_enabled(f_cvb_output_acc_work)) {
+ if (!get_keyval(conf, "accumulatedWork", acc_work))
+ cvm::error("Error: accumulatedWork is missing from the restart.\n");
+ }
}
return COLVARS_OK;
@@ -647,12 +687,21 @@ int colvarbias_restraint_k_moving::set_state_params(std::string const &conf)
std::ostream & colvarbias_restraint_k_moving::write_traj_label(std::ostream &os)
{
+ if (b_chg_force_k && is_enabled(f_cvb_output_acc_work)) {
+ os << " W_"
+ << cvm::wrap_string(this->name, cvm::en_width-2);
+ }
return os;
}
std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os)
{
+ if (b_chg_force_k && is_enabled(f_cvb_output_acc_work)) {
+ os << " "
+ << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+ << acc_work;
+ }
return os;
}
@@ -765,6 +814,7 @@ int colvarbias_restraint_harmonic::update()
// update accumulated work using the current forces
error_code |= colvarbias_restraint_centers_moving::update_acc_work();
+ error_code |= colvarbias_restraint_k_moving::update_acc_work();
return error_code;
}
@@ -876,8 +926,8 @@ colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char co
colvarbias_restraint_moving(key),
colvarbias_restraint_k_moving(key)
{
- lower_wall_k = 0.0;
- upper_wall_k = 0.0;
+ lower_wall_k = -1.0;
+ upper_wall_k = -1.0;
}
@@ -887,26 +937,6 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
colvarbias_restraint_moving::init(conf);
colvarbias_restraint_k_moving::init(conf);
- get_keyval(conf, "lowerWallConstant", lower_wall_k,
- (lower_wall_k > 0.0) ? lower_wall_k : force_k);
- get_keyval(conf, "upperWallConstant", upper_wall_k,
- (upper_wall_k > 0.0) ? upper_wall_k : force_k);
-
- if (lower_wall_k * upper_wall_k > 0.0) {
- for (size_t i = 0; i < num_variables(); i++) {
- if (variables(i)->width != 1.0)
- cvm::log("The lower and upper wall force constants for colvar \""+
- variables(i)->name+
- "\" will be rescaled to "+
- cvm::to_str(lower_wall_k /
- (variables(i)->width * variables(i)->width))+
- " and "+
- cvm::to_str(upper_wall_k /
- (variables(i)->width * variables(i)->width))+
- " according to the specified width.\n");
- }
- }
-
enable(f_cvb_scalar_variables);
size_t i;
@@ -942,16 +972,23 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
}
if ((lower_walls.size() == 0) && (upper_walls.size() == 0)) {
- cvm::error("Error: no walls provided.\n", INPUT_ERROR);
- return INPUT_ERROR;
+ return cvm::error("Error: no walls provided.\n", INPUT_ERROR);
+ }
+
+ if (lower_walls.size() > 0) {
+ get_keyval(conf, "lowerWallConstant", lower_wall_k,
+ (lower_wall_k > 0.0) ? lower_wall_k : force_k);
+ }
+ if (upper_walls.size() > 0) {
+ get_keyval(conf, "upperWallConstant", upper_wall_k,
+ (upper_wall_k > 0.0) ? upper_wall_k : force_k);
}
if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) {
for (i = 0; i < num_variables(); i++) {
if (variables(i)->is_enabled(f_cv_periodic)) {
- cvm::error("Error: at least one variable is periodic, "
- "both walls must be provided.\n", INPUT_ERROR);
- return INPUT_ERROR;
+ return cvm::error("Error: at least one variable is periodic, "
+ "both walls must be provided.\n", INPUT_ERROR);
}
}
}
@@ -972,19 +1009,49 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
INPUT_ERROR);
return INPUT_ERROR;
}
- force_k = lower_wall_k * upper_wall_k;
- // transform the two constants to relative values
+ force_k = std::sqrt(lower_wall_k * upper_wall_k);
+ // transform the two constants to relative values using gemetric mean as ref
+ // to preserve force_k if provided as single parameter
// (allow changing both via force_k)
lower_wall_k /= force_k;
upper_wall_k /= force_k;
+ } else {
+ // If only one wall is defined, need to rescale as well
+ if (lower_walls.size() > 0) {
+ force_k = lower_wall_k;
+ lower_wall_k = 1.0;
+ }
+ if (upper_walls.size() > 0) {
+ force_k = upper_wall_k;
+ upper_wall_k = 1.0;
+ }
}
- for (i = 0; i < num_variables(); i++) {
- if (variables(i)->width != 1.0)
- cvm::log("The force constant for colvar \""+variables(i)->name+
- "\" will be rescaled to "+
- cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+
- " according to the specified width.\n");
+ // Initialize starting value of the force constant (in case it's changing)
+ starting_force_k = force_k;
+
+ if (lower_walls.size() > 0) {
+ for (i = 0; i < num_variables(); i++) {
+ if (variables(i)->width != 1.0)
+ cvm::log("The lower wall force constant for colvar \""+
+ variables(i)->name+
+ "\" will be rescaled to "+
+ cvm::to_str(lower_wall_k * force_k /
+ (variables(i)->width * variables(i)->width))+
+ " according to the specified width.\n");
+ }
+ }
+
+ if (upper_walls.size() > 0) {
+ for (i = 0; i < num_variables(); i++) {
+ if (variables(i)->width != 1.0)
+ cvm::log("The upper wall force constant for colvar \""+
+ variables(i)->name+
+ "\" will be rescaled to "+
+ cvm::to_str(upper_wall_k * force_k /
+ (variables(i)->width * variables(i)->width))+
+ " according to the specified width.\n");
+ }
}
return COLVARS_OK;
@@ -1001,6 +1068,8 @@ int colvarbias_restraint_harmonic_walls::update()
error_code |= colvarbias_restraint::update();
+ error_code |= colvarbias_restraint_k_moving::update_acc_work();
+
return error_code;
}
@@ -1134,6 +1203,7 @@ colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
colvarbias_restraint_centers_moving(key),
colvarbias_restraint_k_moving(key)
{
+ check_positive_k = false;
}
@@ -1177,6 +1247,7 @@ int colvarbias_restraint_linear::update()
// update accumulated work using the current forces
error_code |= colvarbias_restraint_centers_moving::update_acc_work();
+ error_code |= colvarbias_restraint_k_moving::update_acc_work();
return error_code;
}
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index b10649cab112f8a10bfedaff0c10823cc7a43f04..3ee999c2624cfe0cfd9997a2ae096106631444ec 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -89,8 +89,12 @@ public:
virtual int change_configuration(std::string const &conf);
protected:
+
/// \brief Restraint force constant
cvm::real force_k;
+
+ /// \brief Whether the force constant should be positive
+ bool check_positive_k;
};
@@ -129,6 +133,9 @@ protected:
/// \brief Number of steps required to reach the target force constant
/// or restraint centers
long target_nsteps;
+
+ /// \brief Accumulated work (computed when outputAccumulatedWork == true)
+ cvm::real acc_work;
};
@@ -157,8 +164,7 @@ protected:
/// \brief Initial value of the restraint centers
std::vector<colvarvalue> initial_centers;
- /// \brief Amplitude of the restraint centers' increment at each step
- /// towards the new values (calculated from target_nsteps)
+ /// \brief Increment of the restraint centers at each step
std::vector<colvarvalue> centers_incr;
/// \brief Update the centers by interpolating between initial and target
@@ -167,12 +173,6 @@ protected:
/// Whether to write the current restraint centers to the trajectory file
bool b_output_centers;
- /// Whether to write the current accumulated work to the trajectory file
- bool b_output_acc_work;
-
- /// \brief Accumulated work
- cvm::real acc_work;
-
/// Update the accumulated work
int update_acc_work();
};
@@ -212,6 +212,12 @@ protected:
/// \brief Equilibration steps for restraint FE calculation through TI
cvm::real target_equil_steps;
+
+ /// \brief Increment of the force constant at each step
+ cvm::real force_k_incr;
+
+ /// Update the accumulated work
+ int update_acc_work();
};
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index b94d798be9f0a3cd49abc27bb8cdaf41f88e8c87..52078a3a3059cb91315ac9f52e62b2f06e1e55f3 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -20,52 +20,48 @@
// simple_scalar_dist_functions (derived_class)
-#include <fstream>
-#include <cmath>
-
-
#include "colvarmodule.h"
#include "colvar.h"
#include "colvaratoms.h"
-/// \brief Colvar component (base class); most implementations of
-/// \link cvc \endlink utilize one or more \link
-/// colvarmodule::atom \endlink or \link colvarmodule::atom_group
-/// \endlink objects to access atoms.
+/// \brief Colvar component (base class for collective variables)
///
/// A \link cvc \endlink object (or an object of a
-/// cvc-derived class) specifies how to calculate a collective
-/// variable, its gradients and other related physical quantities
-/// which do not depend only on the numeric value (the \link colvar
-/// \endlink class already serves this purpose).
+/// cvc-derived class) implements the calculation of a collective
+/// variable, its gradients and any other related physical quantities
+/// that depend on microscopic degrees of freedom.
+///
+/// No restriction is set to what kind of calculation a \link cvc \endlink
+/// object performs (usually an analytical function of atomic coordinates).
+/// The only constraints are that: \par
///
-/// No restriction is set to what kind of calculation a \link
-/// cvc \endlink object performs (usually calculate an
-/// analytical function of atomic coordinates). The only constraint
-/// is that the value calculated is implemented as a \link colvarvalue
-/// \endlink object. This serves to provide a unique way to calculate
-/// scalar and non-scalar collective variables, and specify if and how
-/// they can be combined together by the parent \link colvar \endlink
-/// object.
+/// - The value is calculated by the \link calc_value() \endlink
+/// method, and is an object of \link colvarvalue \endlink class. This
+/// provides a transparent way to treat scalar and non-scalar variables
+/// alike, and allows an automatic selection of the applicable algorithms.
+///
+/// - The object provides an implementation \link apply_force() \endlink to
+/// apply forces to atoms. Typically, one or more \link cvm::atom_group
+/// \endlink objects are used, but this is not a requirement for as long as
+/// the \link cvc \endlink object communicates with the simulation program.
///
/// <b> If you wish to implement a new collective variable component, you
/// should write your own class by inheriting directly from \link
-/// cvc \endlink, or one of its derived classes (for instance,
-/// \link distance \endlink is frequently used, because it provides
+/// colvar::cvc \endlink, or one of its derived classes (for instance,
+/// \link colvar::distance \endlink is frequently used, because it provides
/// useful data and function members for any colvar based on two
-/// atom groups). The steps are: \par
-/// 1. add the name of this class under colvar.h \par
-/// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par
-/// 3. declare the class in colvarcomp.h \par
-/// 4. implement the class in one of the files colvarcomp_*.C
+/// atom groups).</b>
+///
+/// The steps are: \par
+/// 1. Declare the new class as a derivative of \link colvar::cvc \endlink
+/// in the file \link colvarcomp.h \endlink
+/// 2. Implement the new class in a file named colvarcomp_<something>.cpp
+/// 3. Declare the name of the new class inside the \link colvar \endlink class
+/// in \link colvar.h \endlink (see "list of available components")
+/// 4. Add a call for the new class in colvar::init_components()
+//// (file: colvar.cpp)
///
-/// </b>
-/// The cvm::atom and cvm::atom_group classes are available to
-/// transparently communicate with the simulation program. However,
-/// they are not strictly needed, as long as all the degrees of
-/// freedom associated to the cvc are properly evolved by a simple
-/// call to e.g. apply_force().
class colvar::cvc
: public colvarparse, public colvardeps
@@ -155,7 +151,7 @@ public:
/// \brief Calculate the atomic gradients, to be reused later in
/// order to apply forces
- virtual void calc_gradients() = 0;
+ virtual void calc_gradients() {}
/// \brief Calculate the atomic fit gradients
void calc_fit_gradients();
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index ce8055843f93d389117f38926e81d97c66bc7954..9911f4c87edef9c17605403f21a0ffc1a45826e9 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -581,6 +581,12 @@ colvar::distance_inv::distance_inv(std::string const &conf)
}
}
+ if (is_enabled(f_cvc_debug_gradient)) {
+ cvm::log("Warning: debugGradients will not give correct results "
+ "for distanceInv, because its value and gradients are computed "
+ "simultaneously.\n");
+ }
+
x.type(colvarvalue::type_scalar);
}
@@ -601,11 +607,9 @@ void colvar::distance_inv::calc_value()
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
cvm::rvector const dv = ai2->pos - ai1->pos;
cvm::real const d2 = dv.norm2();
- cvm::real dinv = 1.0;
- for (int ne = 0; ne < exponent/2; ne++)
- dinv *= 1.0/d2;
+ cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
x.real_value += dinv;
- cvm::rvector const dsumddv = -(cvm::real(exponent)) * dinv/d2 * dv;
+ cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
ai1->grad += -1.0 * dsumddv;
ai2->grad += dsumddv;
}
@@ -615,11 +619,9 @@ void colvar::distance_inv::calc_value()
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos);
cvm::real const d2 = dv.norm2();
- cvm::real dinv = 1.0;
- for (int ne = 0; ne < exponent/2; ne++)
- dinv *= 1.0/d2;
+ cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
x.real_value += dinv;
- cvm::rvector const dsumddv = -(cvm::real(exponent)) * dinv/d2 * dv;
+ cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
ai1->grad += -1.0 * dsumddv;
ai2->grad += dsumddv;
}
@@ -627,13 +629,11 @@ void colvar::distance_inv::calc_value()
}
x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
- x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));
-}
+ 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());
+ cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) *
+ cvm::integer_power(x.real_value, exponent+1) /
+ cvm::real(group1->size() * group2->size());
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
ai1->grad *= dxdsum;
}
@@ -643,6 +643,11 @@ void colvar::distance_inv::calc_gradients()
}
+void colvar::distance_inv::calc_gradients()
+{
+}
+
+
void colvar::distance_inv::apply_force(colvarvalue const &force)
{
if (!group1->noforce)
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index ac906e7be755fc382d1cb383730522da38860092..a058ad55c2a4911324c9cd07b4ea7bb1e1afb147 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -8,10 +8,16 @@
// Colvars repository at GitHub.
+#include "colvarmodule.h"
+#include "colvarproxy.h"
#include "colvardeps.h"
+
colvardeps::colvardeps()
- : time_step_factor (1) {}
+{
+ time_step_factor = 1;
+}
+
colvardeps::~colvardeps() {
size_t i;
@@ -416,6 +422,9 @@ void colvardeps::init_cvb_requires() {
init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
f_req_children(f_cvb_get_total_force, f_cv_total_force);
+ init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user);
+ f_req_self(f_cvb_output_acc_work, f_cvb_apply_force);
+
init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index bd892fbca8cef746bb4dfad818b1e4577507a735..940eefb01b1a83c5793f9f7f5de68229b7bfe0ad 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -225,6 +225,8 @@ public:
f_cvb_apply_force,
/// \brief requires total forces
f_cvb_get_total_force,
+ /// \brief whether this bias should record the accumulated work
+ f_cvb_output_acc_work,
/// \brief depends on simulation history
f_cvb_history_dependent,
/// \brief depends on time
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 9016e2c23a0ce573f8711075ee6ecd93b694aa5c..1ac4aae133c3a225aa8d4a8bb2dbdd1c3ac381cb 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -14,6 +14,7 @@
#include "colvarcomp.h"
#include "colvargrid.h"
+#include <ctime>
colvar_grid_count::colvar_grid_count()
: colvar_grid<size_t>()
@@ -22,43 +23,37 @@ colvar_grid_count::colvar_grid_count()
}
colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
- size_t const &def_count)
+ size_t const &def_count)
: colvar_grid<size_t>(nx_i, def_count, 1)
{}
colvar_grid_count::colvar_grid_count(std::vector<colvar *> &colvars,
- size_t const &def_count)
- : colvar_grid<size_t>(colvars, def_count, 1)
+ size_t const &def_count,
+ bool margin)
+ : colvar_grid<size_t>(colvars, def_count, 1, margin)
{}
colvar_grid_scalar::colvar_grid_scalar()
- : colvar_grid<cvm::real>(), samples(NULL), grad(NULL)
+ : colvar_grid<cvm::real>(), samples(NULL)
{}
colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
- : colvar_grid<cvm::real>(g), samples(NULL), grad(NULL)
+ : colvar_grid<cvm::real>(g), samples(NULL)
{
- grad = new cvm::real[nd];
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
- : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL), grad(NULL)
+ : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
{
- grad = new cvm::real[nd];
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
- : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL), grad(NULL)
+ : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
{
- grad = new cvm::real[nd];
}
colvar_grid_scalar::~colvar_grid_scalar()
{
- if (grad) {
- delete [] grad;
- grad = NULL;
- }
}
cvm::real colvar_grid_scalar::maximum_value() const
@@ -143,18 +138,18 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
os << "# xi A(xi)\n";
- if ( cv.size() != 1 ) {
+ if (cv.size() != 1) {
cvm::error("Cannot write integral for multi-dimensional gradient grids.");
return;
}
integral = 0.0;
- int_vals.push_back( 0.0 );
+ int_vals.push_back(0.0);
min = 0.0;
// correction for periodic colvars, so that the PMF is periodic
cvm::real corr;
- if ( periodic[0] ) {
+ if (periodic[0]) {
corr = average();
} else {
corr = 0.0;
@@ -171,7 +166,7 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
}
if ( integral < min ) min = integral;
- int_vals.push_back( integral );
+ int_vals.push_back(integral);
}
bin = 0.0;
@@ -192,3 +187,670 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
+integrate_potential::integrate_potential(std::vector<colvar *> &colvars, colvar_grid_gradient * gradients)
+ : colvar_grid_scalar(colvars, true),
+ gradients(gradients)
+{
+ // parent class colvar_grid_scalar is constructed with margin option set to true
+ // hence PMF grid is wider than gradient grid if non-PBC
+
+ if (nd > 1) {
+ divergence.resize(nt);
+
+ // Compute inverse of Laplacian diagonal for Jacobi preconditioning
+ // For now all code related to preconditioning is commented out
+ // until a method better than Jacobi is implemented
+// cvm::log("Preparing inverse diagonal for preconditioning...");
+// inv_lap_diag.resize(nt);
+// std::vector<cvm::real> id(nt), lap_col(nt);
+// for (int i = 0; i < nt; i++) {
+// if (i % (nt / 100) == 0)
+// cvm::log(cvm::to_str(i));
+// id[i] = 1.;
+// atimes(id, lap_col);
+// id[i] = 0.;
+// inv_lap_diag[i] = 1. / lap_col[i];
+// }
+// cvm::log("Done.");
+ }
+}
+
+
+int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::real & err)
+{
+ int iter = 0;
+
+ if (nd == 1) {
+
+ cvm::real sum = 0.0;
+ cvm::real corr;
+ if ( periodic[0] ) {
+ corr = gradients->average(); // Enforce PBC by subtracting average gradient
+ } else {
+ corr = 0.0;
+ }
+
+ std::vector<int> ix;
+ // Iterate over valid indices in gradient grid
+ for (ix = new_index(); gradients->index_ok(ix); incr(ix)) {
+ set_value(ix, sum);
+ sum += (gradients->value_output(ix) - corr) * widths[0];
+ }
+ if (index_ok(ix)) {
+ // This will happen if non-periodic: then PMF grid has one extra bin wrt gradient grid
+ set_value(ix, sum);
+ }
+
+ } else if (nd <= 3) {
+
+ nr_linbcg_sym(divergence, data, tol, itmax, iter, err);
+ cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err));
+
+ } else {
+ cvm::error("Cannot integrate PMF in dimension > 3\n");
+ }
+
+ return iter;
+}
+
+
+void integrate_potential::set_div()
+{
+ if (nd == 1) return;
+ for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
+ update_div_local(ix);
+ }
+}
+
+
+void integrate_potential::update_div_neighbors(const std::vector<int> &ix0)
+{
+ std::vector<int> ix(ix0);
+ int i, j, k;
+
+ // If not periodic, expanded grid ensures that neighbors of ix0 are valid grid points
+ if (nd == 1) {
+ return;
+
+ } else if (nd == 2) {
+
+ update_div_local(ix);
+ ix[0]++; wrap(ix);
+ update_div_local(ix);
+ ix[1]++; wrap(ix);
+ update_div_local(ix);
+ ix[0]--; wrap(ix);
+ update_div_local(ix);
+
+ } else if (nd == 3) {
+
+ for (i = 0; i<2; i++) {
+ ix[1] = ix0[1];
+ for (j = 0; j<2; j++) {
+ ix[2] = ix0[2];
+ for (k = 0; k<2; k++) {
+ wrap(ix);
+ update_div_local(ix);
+ ix[2]++;
+ }
+ ix[1]++;
+ }
+ ix[0]++;
+ }
+ }
+}
+
+void integrate_potential::get_grad(cvm::real * g, std::vector<int> &ix)
+{
+ size_t count, i;
+ bool edge = gradients->wrap_edge(ix); // Detect edge if non-PBC
+
+ if (gradients->samples) {
+ count = gradients->samples->value(ix);
+ } else {
+ count = 1;
+ }
+
+ if (!edge && count) {
+ cvm::real const *grad = &(gradients->value(ix));
+ cvm::real const fact = 1.0 / count;
+ for ( i = 0; i<nd; i++ ) {
+ g[i] = fact * grad[i];
+ }
+ } else {
+ for ( i = 0; i<nd; i++ ) {
+ g[i] = 0.0;
+ }
+ }
+}
+
+void integrate_potential::update_div_local(const std::vector<int> &ix0)
+{
+ const int linear_index = address(ix0);
+ int i, j, k;
+ std::vector<int> ix = ix0;
+ const cvm::real * g;
+
+ if (nd == 2) {
+ // gradients at grid points surrounding the current scalar grid point
+ cvm::real g00[2], g01[2], g10[2], g11[2];
+
+ get_grad(g11, ix);
+ ix[0] = ix0[0] - 1;
+ get_grad(g01, ix);
+ ix[1] = ix0[1] - 1;
+ get_grad(g00, ix);
+ ix[0] = ix0[0];
+ get_grad(g10, ix);
+
+ divergence[linear_index] = ((g10[0]-g00[0] + g11[0]-g01[0]) / widths[0]
+ + (g01[1]-g00[1] + g11[1]-g10[1]) / widths[1]) * 0.5;
+ } else if (nd == 3) {
+ cvm::real gc[24]; // stores 3d gradients in 8 contiguous bins
+ int index = 0;
+
+ ix[0] = ix0[0] - 1;
+ for (i = 0; i<2; i++) {
+ ix[1] = ix0[1] - 1;
+ for (j = 0; j<2; j++) {
+ ix[2] = ix0[2] - 1;
+ for (k = 0; k<2; k++) {
+ get_grad(gc + index, ix);
+ index += 3;
+ ix[2]++;
+ }
+ ix[1]++;
+ }
+ ix[0]++;
+ }
+
+ divergence[linear_index] =
+ ((gc[3*4]-gc[0] + gc[3*5]-gc[3*1] + gc[3*6]-gc[3*2] + gc[3*7]-gc[3*3])
+ / widths[0]
+ + (gc[3*2+1]-gc[0+1] + gc[3*3+1]-gc[3*1+1] + gc[3*6+1]-gc[3*4+1] + gc[3*7+1]-gc[3*5+1])
+ / widths[1]
+ + (gc[3*1+2]-gc[0+2] + gc[3*3+2]-gc[3*2+2] + gc[3*5+2]-gc[3*4+2] + gc[3*7+2]-gc[3*6+2])
+ / widths[2]) * 0.25;
+ }
+}
+
+
+/// Multiplication by sparse matrix representing Laplacian
+/// NOTE: Laplacian must be symmetric for solving with CG
+void integrate_potential::atimes(const std::vector<cvm::real> &A, std::vector<cvm::real> &LA)
+{
+ if (nd == 2) {
+ // DIMENSION 2
+
+ size_t index, index2;
+ int i, j;
+ cvm::real fact;
+ const cvm::real ffx = 1.0 / (widths[0] * widths[0]);
+ const cvm::real ffy = 1.0 / (widths[1] * widths[1]);
+ const int h = nx[1];
+ const int w = nx[0];
+ // offsets for 4 reference points of the Laplacian stencil
+ int xm = -h;
+ int xp = h;
+ int ym = -1;
+ int yp = 1;
+
+ // NOTE on performance: this version is slightly sub-optimal because
+ // it contains two double loops on the core of the array (for x and y terms)
+ // The slightly faster version is in commit 0254cb5a2958cb2e135f268371c4b45fad34866b
+ // yet it is much uglier, and probably horrible to extend to dimension 3
+ // All terms in the matrix are assigned (=) during the x loops, then updated (+=)
+ // with the y (and z) contributions
+
+
+ // All x components except on x edges
+ index = h; // Skip first column
+
+ // Halve the term on y edges (if any) to preserve symmetry of the Laplacian matrix
+ // (Long Chen, Finite Difference Methods, UCI, 2017)
+ fact = periodic[1] ? 1.0 : 0.5;
+
+ for (i=1; i<w-1; i++) {
+ // Full range of j, but factor may change on y edges (j == 0 and j == h-1)
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ for (j=1; j<h-1; j++) {
+ LA[index] = ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ }
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ }
+ // Edges along x (x components only)
+ index = 0; // Follows left edge
+ index2 = h * (w - 1); // Follows right edge
+ if (periodic[0]) {
+ xm = h * (w - 1);
+ xp = h;
+ fact = periodic[1] ? 1.0 : 0.5;
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ for (j=1; j<h-1; j++) {
+ LA[index] = ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ }
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ } else {
+ xm = -h;
+ xp = h;
+ fact = periodic[1] ? 1.0 : 0.5; // Halve in corners in full PBC only
+ // lower corner, "j == 0"
+ LA[index] = fact * ffx * (A[index + xp] - A[index]);
+ LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
+ index++;
+ index2++;
+ for (j=1; j<h-1; j++) {
+ // x gradient (+ y term of laplacian, calculated below)
+ LA[index] = ffx * (A[index + xp] - A[index]);
+ LA[index2] = ffx * (A[index2 + xm] - A[index2]);
+ index++;
+ index2++;
+ }
+ // upper corner, j == h-1
+ LA[index] = fact * ffx * (A[index + xp] - A[index]);
+ LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
+ }
+
+ // Now adding all y components
+ // All y components except on y edges
+ index = 1; // Skip first element (in first row)
+
+ fact = periodic[0] ? 1.0 : 0.5; // for i == 0
+ for (i=0; i<w; i++) {
+ // Factor of 1/2 on x edges if non-periodic
+ if (i == 1) fact = 1.0;
+ if (i == w - 1) fact = periodic[0] ? 1.0 : 0.5;
+ for (j=1; j<h-1; j++) {
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ index++;
+ }
+ index += 2; // skip the edges and move to next column
+ }
+ // Edges along y (y components only)
+ index = 0; // Follows bottom edge
+ index2 = h - 1; // Follows top edge
+ if (periodic[1]) {
+ fact = periodic[0] ? 1.0 : 0.5;
+ ym = h - 1;
+ yp = 1;
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ index += h;
+ index2 += h;
+ for (i=1; i<w-1; i++) {
+ LA[index] += ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ index += h;
+ index2 += h;
+ }
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ } else {
+ ym = -1;
+ yp = 1;
+ fact = periodic[0] ? 1.0 : 0.5; // Halve in corners in full PBC only
+ // Left corner
+ LA[index] += fact * ffy * (A[index + yp] - A[index]);
+ LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
+ index += h;
+ index2 += h;
+ for (i=1; i<w-1; i++) {
+ // y gradient (+ x term of laplacian, calculated above)
+ LA[index] += ffy * (A[index + yp] - A[index]);
+ LA[index2] += ffy * (A[index2 + ym] - A[index2]);
+ index += h;
+ index2 += h;
+ }
+ // Right corner
+ LA[index] += fact * ffy * (A[index + yp] - A[index]);
+ LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
+ }
+
+ } else if (nd == 3) {
+ // DIMENSION 3
+
+ int i, j, k;
+ size_t index, index2;
+ cvm::real fact = 1.0;
+ const cvm::real ffx = 1.0 / (widths[0] * widths[0]);
+ const cvm::real ffy = 1.0 / (widths[1] * widths[1]);
+ const cvm::real ffz = 1.0 / (widths[2] * widths[2]);
+ const int h = nx[2]; // height
+ const int d = nx[1]; // depth
+ const int w = nx[0]; // width
+ // offsets for 6 reference points of the Laplacian stencil
+ int xm = -d * h;
+ int xp = d * h;
+ int ym = -h;
+ int yp = h;
+ int zm = -1;
+ int zp = 1;
+
+ cvm::real factx = periodic[0] ? 1 : 0.5; // factor to be applied on x edges
+ cvm::real facty = periodic[1] ? 1 : 0.5; // same for y
+ cvm::real factz = periodic[2] ? 1 : 0.5; // same for z
+ cvm::real ifactx = 1 / factx;
+ cvm::real ifacty = 1 / facty;
+ cvm::real ifactz = 1 / factz;
+
+ // All x components except on x edges
+ index = d * h; // Skip left slab
+ fact = facty * factz;
+ for (i=1; i<w-1; i++) {
+ for (j=0; j<d; j++) { // full range of y
+ if (j == 1) fact *= ifacty;
+ if (j == d-1) fact *= facty;
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) { // full range of z
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ }
+ fact *= factz;
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ index++;
+ }
+ }
+ // Edges along x (x components only)
+ index = 0; // Follows left slab
+ index2 = d * h * (w - 1); // Follows right slab
+ if (periodic[0]) {
+ xm = d * h * (w - 1);
+ xp = d * h;
+ fact = facty * factz;
+ for (j=0; j<d; j++) {
+ if (j == 1) fact *= ifacty;
+ if (j == d-1) fact *= facty;
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) {
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ }
+ fact *= factz;
+ LA[index] = fact * ffx * (A[index + xm] + A[index + xp] - 2.0 * A[index]);
+ LA[index2] = fact * ffx * (A[index2 - xp] + A[index2 - xm] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ }
+ } else {
+ xm = -d * h;
+ xp = d * h;
+ fact = facty * factz;
+ for (j=0; j<d; j++) {
+ if (j == 1) fact *= ifacty;
+ if (j == d-1) fact *= facty;
+ LA[index] = fact * ffx * (A[index + xp] - A[index]);
+ LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
+ index++;
+ index2++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) {
+ // x gradient (+ y, z terms of laplacian, calculated below)
+ LA[index] = fact * ffx * (A[index + xp] - A[index]);
+ LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
+ index++;
+ index2++;
+ }
+ fact *= factz;
+ LA[index] = fact * ffx * (A[index + xp] - A[index]);
+ LA[index2] = fact * ffx * (A[index2 + xm] - A[index2]);
+ index++;
+ index2++;
+ }
+ }
+
+ // Now adding all y components
+ // All y components except on y edges
+ index = h; // Skip first column (in front slab)
+ fact = factx * factz;
+ for (i=0; i<w; i++) { // full range of x
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ for (j=1; j<d-1; j++) {
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ index++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) {
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ index++;
+ }
+ fact *= factz;
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ index++;
+ }
+ index += 2 * h; // skip columns in front and back slabs
+ }
+ // Edges along y (y components only)
+ index = 0; // Follows front slab
+ index2 = h * (d - 1); // Follows back slab
+ if (periodic[1]) {
+ ym = h * (d - 1);
+ yp = h;
+ fact = factx * factz;
+ for (i=0; i<w; i++) {
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) {
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ }
+ fact *= factz;
+ LA[index] += fact * ffy * (A[index + ym] + A[index + yp] - 2.0 * A[index]);
+ LA[index2] += fact * ffy * (A[index2 - yp] + A[index2 - ym] - 2.0 * A[index2]);
+ index++;
+ index2++;
+ index += h * (d - 1);
+ index2 += h * (d - 1);
+ }
+ } else {
+ ym = -h;
+ yp = h;
+ fact = factx * factz;
+ for (i=0; i<w; i++) {
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ LA[index] += fact * ffy * (A[index + yp] - A[index]);
+ LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
+ index++;
+ index2++;
+ fact *= ifactz;
+ for (k=1; k<h-1; k++) {
+ // y gradient (+ x, z terms of laplacian, calculated above and below)
+ LA[index] += fact * ffy * (A[index + yp] - A[index]);
+ LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
+ index++;
+ index2++;
+ }
+ fact *= factz;
+ LA[index] += fact * ffy * (A[index + yp] - A[index]);
+ LA[index2] += fact * ffy * (A[index2 + ym] - A[index2]);
+ index++;
+ index2++;
+ index += h * (d - 1);
+ index2 += h * (d - 1);
+ }
+ }
+
+ // Now adding all z components
+ // All z components except on z edges
+ index = 1; // Skip first element (in bottom slab)
+ fact = factx * facty;
+ for (i=0; i<w; i++) { // full range of x
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ for (k=1; k<h-1; k++) {
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ index++;
+ }
+ fact *= ifacty;
+ index += 2; // skip edge slabs
+ for (j=1; j<d-1; j++) { // full range of y
+ for (k=1; k<h-1; k++) {
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ index++;
+ }
+ index += 2; // skip edge slabs
+ }
+ fact *= facty;
+ for (k=1; k<h-1; k++) {
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ index++;
+ }
+ index += 2; // skip edge slabs
+ }
+ // Edges along z (z components onlz)
+ index = 0; // Follows bottom slab
+ index2 = h - 1; // Follows top slab
+ if (periodic[2]) {
+ zm = h - 1;
+ zp = 1;
+ fact = factx * facty;
+ for (i=0; i<w; i++) {
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
+ index += h;
+ index2 += h;
+ fact *= ifacty;
+ for (j=1; j<d-1; j++) {
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
+ index += h;
+ index2 += h;
+ }
+ fact *= facty;
+ LA[index] += fact * ffz * (A[index + zm] + A[index + zp] - 2.0 * A[index]);
+ LA[index2] += fact * ffz * (A[index2 - zp] + A[index2 - zm] - 2.0 * A[index2]);
+ index += h;
+ index2 += h;
+ }
+ } else {
+ zm = -1;
+ zp = 1;
+ fact = factx * facty;
+ for (i=0; i<w; i++) {
+ if (i == 1) fact *= ifactx;
+ if (i == w-1) fact *= factx;
+ LA[index] += fact * ffz * (A[index + zp] - A[index]);
+ LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
+ index += h;
+ index2 += h;
+ fact *= ifacty;
+ for (j=1; j<d-1; j++) {
+ // z gradient (+ x, y terms of laplacian, calculated above)
+ LA[index] += fact * ffz * (A[index + zp] - A[index]);
+ LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
+ index += h;
+ index2 += h;
+ }
+ fact *= facty;
+ LA[index] += fact * ffz * (A[index + zp] - A[index]);
+ LA[index2] += fact * ffz * (A[index2 + zm] - A[index2]);
+ index += h;
+ index2 += h;
+ }
+ }
+ }
+}
+
+
+/*
+/// Inversion of preconditioner matrix (e.g. diagonal of the Laplacian)
+void integrate_potential::asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x)
+{
+ for (size_t i=0; i<nt; i++) {
+ x[i] = b[i] * inv_lap_diag[i]; // Jacobi preconditioner - little benefit in tests so far
+ }
+ return;
+}*/
+
+
+// b : RHS of equation
+// x : initial guess for the solution; output is solution
+// itol : convergence criterion
+void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x, const cvm::real &tol,
+ const int itmax, int &iter, cvm::real &err)
+{
+ cvm::real ak,akden,bk,bkden,bknum,bnrm;
+ const cvm::real EPS=1.0e-14;
+ int j;
+ std::vector<cvm::real> p(nt), r(nt), z(nt);
+
+ iter=0;
+ atimes(x,r);
+ for (j=0;j<nt;j++) {
+ r[j]=b[j]-r[j];
+ }
+ bnrm=l2norm(b);
+ if (bnrm < EPS) {
+ return; // Target is zero, will break relative error calc
+ }
+// asolve(r,z); // precon
+ bkden = 1.0;
+ while (iter < itmax) {
+ ++iter;
+ for (bknum=0.0,j=0;j<nt;j++) {
+ bknum += r[j]*r[j]; // precon: z[j]*r[j]
+ }
+ if (iter == 1) {
+ for (j=0;j<nt;j++) {
+ p[j] = r[j]; // precon: p[j] = z[j]
+ }
+ } else {
+ bk=bknum/bkden;
+ for (j=0;j<nt;j++) {
+ p[j] = bk*p[j] + r[j]; // precon: bk*p[j] + z[j]
+ }
+ }
+ bkden = bknum;
+ atimes(p,z);
+ for (akden=0.0,j=0;j<nt;j++) {
+ akden += z[j]*p[j];
+ }
+ ak = bknum/akden;
+ for (j=0;j<nt;j++) {
+ x[j] += ak*p[j];
+ r[j] -= ak*z[j];
+ }
+// asolve(r,z); // precon
+ err = l2norm(r)/bnrm;
+ if (cvm::debug())
+ std::cout << "iter=" << std::setw(4) << iter+1 << std::setw(12) << err << std::endl;
+ if (err <= tol)
+ break;
+ }
+}
+
+cvm::real integrate_potential::l2norm(const std::vector<cvm::real> &x)
+{
+ size_t i;
+ cvm::real sum = 0.0;
+ for (i=0;i<x.size();i++)
+ sum += x[i]*x[i];
+ return sqrt(sum);
+}
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index a01104dba819c2745a9d63d434b8e6832c9d85b7..9a0fe4c8ecb805045a6e19caa4f47aefc9c37434 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -97,8 +97,8 @@ public:
/// Whether this grid has been filled with data or is still empty
bool has_data;
- /// Return the number of colvars
- inline size_t number_of_colvars() const
+ /// Return the number of colvar objects
+ inline size_t num_variables() const
{
return nd;
}
@@ -374,6 +374,20 @@ public:
}
}
+ /// Wrap an index vector around periodic boundary conditions
+ /// or detects edges if non-periodic
+ inline bool wrap_edge(std::vector<int> & ix) const
+ {
+ bool edge = false;
+ for (size_t i = 0; i < nd; i++) {
+ if (periodic[i]) {
+ ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
+ } else if (ix[i] == -1 || ix[i] == nx[i]) {
+ edge = true;
+ }
+ }
+ return edge;
+ }
/// \brief Report the bin corresponding to the current value of variable i
inline int current_bin_scalar(int const i) const
@@ -445,6 +459,12 @@ public:
has_data = true;
}
+ /// Set the value at the point with linear address i (for speed)
+ inline void set_value(size_t i, T const &t)
+ {
+ data[i] = t;
+ }
+
/// \brief Get the change from this to other_grid
/// and store the result in this.
/// this_grid := other_grid - this_grid
@@ -518,6 +538,11 @@ public:
return data[this->address(ix) + imult];
}
+ /// \brief Get the binned value indexed by linear address i
+ inline T const & value(size_t i) const
+ {
+ return data[i];
+ }
/// \brief Add a constant to all elements (fast loop)
inline void add_constant(T const &t)
@@ -1110,20 +1135,20 @@ public:
// write the header
os << "object 1 class gridpositions counts";
size_t icv;
- for (icv = 0; icv < number_of_colvars(); icv++) {
+ for (icv = 0; icv < num_variables(); icv++) {
os << " " << number_of_points(icv);
}
os << "\n";
os << "origin";
- for (icv = 0; icv < number_of_colvars(); icv++) {
+ for (icv = 0; icv < num_variables(); icv++) {
os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]);
}
os << "\n";
- for (icv = 0; icv < number_of_colvars(); icv++) {
+ for (icv = 0; icv < num_variables(); icv++) {
os << "delta";
- for (size_t icv2 = 0; icv2 < number_of_colvars(); icv2++) {
+ for (size_t icv2 = 0; icv2 < num_variables(); icv2++) {
if (icv == icv2) os << " " << widths[icv];
else os << " " << 0.0;
}
@@ -1131,7 +1156,7 @@ public:
}
os << "object 2 class gridconnections counts";
- for (icv = 0; icv < number_of_colvars(); icv++) {
+ for (icv = 0; icv < num_variables(); icv++) {
os << " " << number_of_points(icv);
}
os << "\n";
@@ -1167,7 +1192,8 @@ public:
/// Constructor from a vector of colvars
colvar_grid_count(std::vector<colvar *> &colvars,
- size_t const &def_count = 0);
+ size_t const &def_count = 0,
+ bool margin = false);
/// Increment the counter at given position
inline void incr_count(std::vector<int> const &ix)
@@ -1210,12 +1236,13 @@ public:
int A0, A1, A2;
std::vector<int> ix = ix0;
+ // TODO this can be rewritten more concisely with wrap_edge()
if (periodic[n]) {
ix[n]--; wrap(ix);
- A0 = data[address(ix)];
+ A0 = value(ix);
ix = ix0;
ix[n]++; wrap(ix);
- A1 = data[address(ix)];
+ A1 = value(ix);
if (A0 * A1 == 0) {
return 0.; // can't handle empty bins
} else {
@@ -1224,10 +1251,10 @@ public:
}
} else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
ix[n]--;
- A0 = data[address(ix)];
+ A0 = value(ix);
ix = ix0;
ix[n]++;
- A1 = data[address(ix)];
+ A1 = value(ix);
if (A0 * A1 == 0) {
return 0.; // can't handle empty bins
} else {
@@ -1238,9 +1265,9 @@ public:
// 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)];
+ A0 = value(ix);
+ ix[n] += increment; A1 = value(ix);
+ ix[n] += increment; A2 = value(ix);
if (A0 * A1 * A2 == 0) {
return 0.; // can't handle empty bins
} else {
@@ -1249,6 +1276,49 @@ public:
}
}
}
+
+ /// \brief Return the gradient of discrete count from finite differences
+ /// on the *same* grid for dimension n
+ inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
+ int n = 0)
+ {
+ int A0, A1, A2;
+ std::vector<int> ix = ix0;
+
+ // FIXME this can be rewritten more concisely with wrap_edge()
+ if (periodic[n]) {
+ ix[n]--; wrap(ix);
+ A0 = value(ix);
+ ix = ix0;
+ ix[n]++; wrap(ix);
+ A1 = value(ix);
+ if (A0 * A1 == 0) {
+ return 0.; // can't handle empty bins
+ } else {
+ return cvm::real(A1 - A0) / (widths[n] * 2.);
+ }
+ } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
+ ix[n]--;
+ A0 = value(ix);
+ ix = ix0;
+ ix[n]++;
+ A1 = value(ix);
+ if (A0 * A1 == 0) {
+ return 0.; // can't handle empty bins
+ } else {
+ return cvm::real(A1 - A0) / (widths[n] * 2.);
+ }
+ } else {
+ // edge: use 2nd order derivative
+ int increment = (ix[n] == 0 ? 1 : -1);
+ // move right from left edge, or the other way around
+ A0 = value(ix);
+ ix[n] += increment; A1 = value(ix);
+ ix[n] += increment; A2 = value(ix);
+ return (-1.5 * cvm::real(A0) + 2. * cvm::real(A1)
+ - 0.5 * cvm::real(A2)) * increment / widths[n];
+ }
+ }
};
@@ -1289,27 +1359,57 @@ public:
has_data = true;
}
- /// Return the gradient of the scalar field from finite differences
- inline const cvm::real * gradient_finite_diff( const std::vector<int> &ix0 )
+ /// \brief Return the gradient of the scalar field from finite differences
+ /// Input coordinates are those of gradient grid, shifted wrt scalar grid
+ /// Should not be called on edges of scalar grid, provided the latter has margins
+ /// wrt gradient grid
+ inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad)
{
cvm::real A0, A1;
std::vector<int> ix;
- if (nd != 2) {
- cvm::error("Finite differences available in dimension 2 only.");
- return grad;
- }
- for (unsigned int n = 0; n < nd; n++) {
+ size_t i, j, k, n;
+
+ if (nd == 2) {
+ for (n = 0; n < 2; n++) {
+ ix = ix0;
+ A0 = value(ix);
+ ix[n]++; wrap(ix);
+ A1 = value(ix);
+ ix[1-n]++; wrap(ix);
+ A1 += value(ix);
+ ix[n]--; wrap(ix);
+ A0 += value(ix);
+ grad[n] = 0.5 * (A1 - A0) / widths[n];
+ }
+ } else if (nd == 3) {
+
+ cvm::real p[8]; // potential values within cube, indexed in binary (4 i + 2 j + k)
ix = ix0;
- A0 = data[address(ix)];
- ix[n]++; wrap(ix);
- A1 = data[address(ix)];
- ix[1-n]++; wrap(ix);
- A1 += data[address(ix)];
- ix[n]--; wrap(ix);
- A0 += data[address(ix)];
- grad[n] = 0.5 * (A1 - A0) / widths[n];
+ int index = 0;
+ for (i = 0; i<2; i++) {
+ ix[1] = ix0[1];
+ for (j = 0; j<2; j++) {
+ ix[2] = ix0[2];
+ for (k = 0; k<2; k++) {
+ wrap(ix);
+ p[index++] = value(ix);
+ ix[2]++;
+ }
+ ix[1]++;
+ }
+ ix[0]++;
+ }
+
+ // The following would be easier to read using binary literals
+ // 100 101 110 111 000 001 010 011
+ grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0];
+ // 010 011 110 111 000 001 100 101
+ grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[0];
+ // 001 011 101 111 000 010 100 110
+ grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[0];
+ } else {
+ cvm::error("Finite differences available in dimension 2 and 3 only.");
}
- return grad;
}
/// \brief Return the value of the function at ix divided by its
@@ -1373,10 +1473,6 @@ public:
/// \brief Assuming that the map is a normalized probability density,
/// calculates the entropy (uses widths if they are defined)
cvm::real entropy() const;
-
-private:
- // gradient
- cvm::real * grad;
};
@@ -1390,6 +1486,10 @@ public:
/// should be divided
colvar_grid_count *samples;
+ /// \brief Provide the floating point weights by which each binned value
+ /// should be divided (alternate to samples, only one should be non-null)
+ colvar_grid_scalar *weights;
+
/// Default constructor
colvar_grid_gradient();
@@ -1403,6 +1503,29 @@ public:
/// Constructor from a vector of colvars
colvar_grid_gradient(std::vector<colvar *> &colvars);
+ /// \brief Get a vector with the binned value(s) indexed by ix, normalized if applicable
+ inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const
+ {
+ cvm::real const * p = &value(ix);
+ if (samples) {
+ int count = samples->value(ix);
+ if (count) {
+ cvm::real invcount = 1.0 / count;
+ for (size_t i = 0; i < mult; i++) {
+ v[i] = invcount * p[i];
+ }
+ } else {
+ for (size_t i = 0; i < mult; i++) {
+ v[i] = 0.0;
+ }
+ }
+ } else {
+ for (size_t i = 0; i < mult; i++) {
+ v[i] = p[i];
+ }
+ }
+ }
+
/// \brief Accumulate the value
inline void acc_value(std::vector<int> const &ix, std::vector<colvarvalue> const &values) {
for (size_t imult = 0; imult < mult; imult++) {
@@ -1412,23 +1535,25 @@ public:
samples->incr_count(ix);
}
- /// \brief Accumulate the gradient
- inline void acc_grad(std::vector<int> const &ix, cvm::real const *grads) {
+ /// \brief Accumulate the gradient based on the force (i.e. sums the
+ /// opposite of the force)
+ inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
for (size_t imult = 0; imult < mult; imult++) {
- data[address(ix) + imult] += grads[imult];
+ data[address(ix) + imult] -= forces[imult];
}
if (samples)
samples->incr_count(ix);
}
/// \brief Accumulate the gradient based on the force (i.e. sums the
- /// opposite of the force)
- inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
+ /// opposite of the force) with a non-integer weight
+ inline void acc_force_weighted(std::vector<int> const &ix,
+ cvm::real const *forces,
+ cvm::real weight) {
for (size_t imult = 0; imult < mult; imult++) {
- data[address(ix) + imult] -= forces[imult];
+ data[address(ix) + imult] -= forces[imult] * weight;
}
- if (samples)
- samples->incr_count(ix);
+ weights->acc_value(ix, weight);
}
/// \brief Return the value of the function at ix divided by its
@@ -1498,5 +1623,70 @@ public:
};
+
+/// Integrate (1D, 2D or 3D) gradients
+
+class integrate_potential : public colvar_grid_scalar
+{
+ public:
+
+ integrate_potential();
+
+ virtual ~integrate_potential()
+ {}
+
+ /// Constructor from a vector of colvars + gradient grid
+ integrate_potential (std::vector<colvar *> &colvars, colvar_grid_gradient * gradients);
+
+ /// \brief Calculate potential from divergence (in 2D); return number of steps
+ int integrate (const int itmax, const cvm::real & tol, cvm::real & err);
+
+ /// \brief Update matrix containing divergence and boundary conditions
+ /// based on new gradient point value, in neighboring bins
+ void update_div_neighbors(const std::vector<int> &ix);
+
+ /// \brief Set matrix containing divergence and boundary conditions
+ /// based on complete gradient grid
+ void set_div();
+
+ /// \brief Add constant to potential so that its minimum value is zero
+ /// Useful e.g. for output
+ inline void set_zero_minimum() {
+ add_constant(-1.0 * minimum_value());
+ }
+
+ protected:
+
+ // Reference to gradient grid
+ colvar_grid_gradient *gradients;
+
+ /// Array holding divergence + boundary terms (modified Neumann) if not periodic
+ std::vector<cvm::real> divergence;
+
+// std::vector<cvm::real> inv_lap_diag; // Inverse of the diagonal of the Laplacian; for conditioning
+
+ /// \brief Update matrix containing divergence and boundary conditions
+ /// called by update_div_neighbors
+ void update_div_local(const std::vector<int> &ix);
+
+ /// Obtain the gradient vector at given location ix, if available
+ /// or zero if it is on the edge of the gradient grid
+ /// ix gets wrapped in PBC
+ void get_grad(cvm::real * g, std::vector<int> &ix);
+
+ /// \brief Solve linear system based on CG, valid for symmetric matrices only
+ void nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x,
+ const cvm::real &tol, const int itmax, int &iter, cvm::real &err);
+
+ /// l2 norm of a vector
+ cvm::real l2norm(const std::vector<cvm::real> &x);
+
+ /// Multiplication by sparse matrix representing Lagrangian (or its transpose)
+ void atimes(const std::vector<cvm::real> &x, std::vector<cvm::real> &r);
+
+// /// Inversion of preconditioner matrix
+// void asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x);
+};
+
#endif
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 200c2d6848f10cf6a52ab70745c6049c9387dff7..9898e2d5e8264d96ba79d12aac5e70bf9dbfd67d 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -24,6 +24,7 @@
#include "colvaratoms.h"
#include "colvarcomp.h"
+
colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
depth_s = 0;
@@ -417,10 +418,10 @@ int colvarmodule::parse_biases(std::string const &conf)
"Please ensure that their forces do not counteract each other.\n");
}
- if (biases.size() || use_scripted_forces) {
+ if (num_biases() || use_scripted_forces) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables biases initialized, "+
- cvm::to_str(biases.size())+" in total.\n");
+ cvm::to_str(num_biases())+" in total.\n");
} else {
if (!use_scripted_forces) {
cvm::log("No collective variables biases were defined.\n");
@@ -431,12 +432,37 @@ int colvarmodule::parse_biases(std::string const &conf)
}
+int colvarmodule::num_variables() const
+{
+ return colvars.size();
+}
+
+
+int colvarmodule::num_variables_feature(int feature_id) const
+{
+ size_t n = 0;
+ for (std::vector<colvar *>::const_iterator cvi = colvars.begin();
+ cvi != colvars.end();
+ cvi++) {
+ if ((*cvi)->is_enabled(feature_id)) {
+ n++;
+ }
+ }
+ return n;
+}
+
+
+int colvarmodule::num_biases() const
+{
+ return biases.size();
+}
+
+
int colvarmodule::num_biases_feature(int feature_id) const
{
- colvarmodule *cv = cvm::main();
size_t n = 0;
- for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
- bi != cv->biases.end();
+ for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
+ bi != biases.end();
bi++) {
if ((*bi)->is_enabled(feature_id)) {
n++;
@@ -448,10 +474,9 @@ int colvarmodule::num_biases_feature(int feature_id) const
int colvarmodule::num_biases_type(std::string const &type) const
{
- colvarmodule *cv = cvm::main();
size_t n = 0;
- for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
- bi != cv->biases.end();
+ for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
+ bi != biases.end();
bi++) {
if ((*bi)->bias_type == type) {
n++;
@@ -465,7 +490,7 @@ std::vector<std::string> const colvarmodule::time_dependent_biases() const
{
size_t i;
std::vector<std::string> biases_names;
- for (i = 0; i < biases.size(); i++) {
+ for (i = 0; i < num_biases(); i++) {
if (biases[i]->is_enabled(colvardeps::f_cvb_apply_force) &&
biases[i]->is_enabled(colvardeps::f_cvb_active) &&
(biases[i]->is_enabled(colvardeps::f_cvb_history_dependent) ||
@@ -790,7 +815,7 @@ int colvarmodule::calc_biases()
{
// update the biases and communicate their forces to the collective
// variables
- if (cvm::debug() && biases.size())
+ if (cvm::debug() && num_biases())
cvm::log("Updating collective variable biases.\n");
std::vector<colvarbias *>::iterator bi;
@@ -852,7 +877,7 @@ int colvarmodule::update_colvar_forces()
std::vector<colvarbias *>::iterator bi;
// sum the forces from all biases for each collective variable
- if (cvm::debug() && biases.size())
+ if (cvm::debug() && num_biases())
cvm::log("Collecting forces from all biases.\n");
cvm::increase_depth();
for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
@@ -1073,8 +1098,6 @@ int colvarmodule::reset()
int colvarmodule::setup_input()
{
- if (this->size() == 0) return cvm::get_error();
-
std::string restart_in_name("");
// read the restart configuration, if available
@@ -1107,14 +1130,12 @@ int colvarmodule::setup_input()
}
}
- return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+ return cvm::get_error();
}
int colvarmodule::setup_output()
{
- if (this->size() == 0) return cvm::get_error();
-
int error_code = COLVARS_OK;
// output state file (restart)
@@ -1123,7 +1144,8 @@ int colvarmodule::setup_output()
std::string("");
if (restart_out_name.size()) {
- cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
+ cvm::log("The restart output state file will be \""+
+ restart_out_name+"\".\n");
}
output_prefix() = proxy->output_prefix();
@@ -1154,7 +1176,7 @@ int colvarmodule::setup_output()
set_error_bits(FILE_ERROR);
}
- return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+ return cvm::get_error();
}
@@ -1738,6 +1760,89 @@ int cvm::load_coords_xyz(char const *filename,
}
+
+// Wrappers to proxy functions: these may go in the future
+
+cvm::real cvm::unit_angstrom()
+{
+ return proxy->unit_angstrom();
+}
+
+
+cvm::real cvm::boltzmann()
+{
+ return proxy->boltzmann();
+}
+
+
+cvm::real cvm::temperature()
+{
+ return proxy->temperature();
+}
+
+
+cvm::real cvm::dt()
+{
+ return proxy->dt();
+}
+
+
+void cvm::request_total_force()
+{
+ proxy->request_total_force(true);
+}
+
+cvm::rvector cvm::position_distance(atom_pos const &pos1,
+ atom_pos const &pos2)
+{
+ return proxy->position_distance(pos1, pos2);
+}
+
+cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
+ cvm::atom_pos const &pos2)
+{
+ return proxy->position_dist2(pos1, pos2);
+}
+
+cvm::real cvm::rand_gaussian(void)
+{
+ return proxy->rand_gaussian();
+}
+
+
+
+bool cvm::replica_enabled()
+{
+ return proxy->replica_enabled();
+}
+
+int cvm::replica_index()
+{
+ return proxy->replica_index();
+}
+
+int cvm::replica_num()
+{
+ return proxy->replica_num();
+}
+
+void cvm::replica_comm_barrier()
+{
+ return proxy->replica_comm_barrier();
+}
+
+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);
+}
+
+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);
+}
+
+
+
// shared pointer to the proxy object
colvarproxy *colvarmodule::proxy = NULL;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 14e5d56701f75730c7bb1ce07764c89edd0b665b..64a98b77a02d59eae598a3f223214fc103018263 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -39,16 +39,14 @@ You can browse the class hierarchy or the list of source files.
#define FILE_ERROR (1<<4)
#define MEMORY_ERROR (1<<5)
#define FATAL_ERROR (1<<6) // Should be set, or not, together with other bits
-#define DELETE_COLVARS (1<<7) // Instruct the caller to delete cvm
+//#define DELETE_COLVARS (1<<7) // Instruct the caller to delete cvm
#define COLVARS_NO_SUCH_FRAME (1<<8) // Cannot load the requested frame
#include <iostream>
#include <iomanip>
-#include <string>
-#include <cstring>
-#include <sstream>
#include <fstream>
-#include <cmath>
+#include <sstream>
+#include <string>
#include <vector>
#include <list>
@@ -84,12 +82,18 @@ public:
/// Defining an abstract real number allows to switch precision
typedef double real;
- /// Override std::pow with a product for n positive integer
- static inline real integer_power(real x, int n)
+ /// Override std::pow with a product for n integer
+ static inline real integer_power(real const &x, int const n)
{
- real result = 1.0;
- for (int i = 0; i < n; i++) result *= x;
- return result;
+ // Original code: math_special.h in LAMMPS
+ double yy, ww;
+ if (x == 0.0) return 0.0;
+ int nn = (n > 0) ? n : -n;
+ ww = x;
+ for (yy = 1.0; nn != 0; nn >>= 1, ww *=ww) {
+ if (nn & 1) yy *= ww;
+ }
+ return (n > 0) ? yy : 1.0/yy;
}
/// Residue identifier
@@ -301,13 +305,23 @@ private:
public:
+ /// Return how many variables are defined
+ int num_variables() const;
+
+ /// Return how many variables have this feature enabled
+ int num_variables_feature(int feature_id) const;
+
+ /// Return how many biases are defined
+ int num_biases() const;
+
/// Return how many biases have this feature enabled
int num_biases_feature(int feature_id) const;
- /// Return how many biases are defined with this type
+ /// Return how many biases of this type are defined
int num_biases_type(std::string const &type) const;
- /// Return the names of time-dependent biases with forces enabled
+ /// Return the names of time-dependent biases with forces enabled (ABF,
+ /// metadynamics, etc)
std::vector<std::string> const time_dependent_biases() const;
private:
@@ -602,16 +616,14 @@ public:
typedef colvarmodule cvm;
-#include "colvartypes.h"
-
std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);
template<typename T> std::string cvm::to_str(T const &x,
- size_t const &width,
- size_t const &prec) {
+ size_t const &width,
+ size_t const &prec) {
std::ostringstream os;
if (width) os.width(width);
if (prec) {
@@ -622,9 +634,10 @@ template<typename T> std::string cvm::to_str(T const &x,
return os.str();
}
+
template<typename T> std::string cvm::to_str(std::vector<T> const &x,
- size_t const &width,
- size_t const &prec) {
+ size_t const &width,
+ size_t const &prec) {
if (!x.size()) return std::string("");
std::ostringstream os;
if (prec) {
@@ -645,70 +658,4 @@ template<typename T> std::string cvm::to_str(std::vector<T> const &x,
}
-#include "colvarproxy.h"
-
-
-inline cvm::real cvm::unit_angstrom()
-{
- return proxy->unit_angstrom();
-}
-
-inline cvm::real cvm::boltzmann()
-{
- return proxy->boltzmann();
-}
-
-inline cvm::real cvm::temperature()
-{
- return proxy->temperature();
-}
-
-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_total_force()
-{
- proxy->request_total_force(true);
-}
-
-inline cvm::rvector cvm::position_distance(atom_pos const &pos1,
- atom_pos const &pos2)
-{
- return proxy->position_distance(pos1, pos2);
-}
-
-inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
- cvm::atom_pos const &pos2)
-{
- return proxy->position_dist2(pos1, pos2);
-}
-
-inline cvm::real cvm::rand_gaussian(void)
-{
- return proxy->rand_gaussian();
-}
-
#endif
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index 9f333b7b76bc029f142fa5443d16c3d73249b77a..43c8c69c040e6b6fdb14a31249f39a293b997d8a 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -553,7 +553,8 @@ bool colvarparse::key_lookup(std::string const &conf,
size_t *save_pos)
{
if (cvm::debug()) {
- cvm::log("Looking for the keyword \""+std::string(key_in)+"\" and its value.\n");
+ cvm::log("Looking for the keyword \""+std::string(key_in)+
+ "\" and its value.\n");
}
// add this keyword to the register (in its camelCase version)
diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp
index 8160144c6bab1e11d36f4d1b7712e75bce90e6c8..8767d5f459090854cf19ad011f579e752c6d2230 100644
--- a/lib/colvars/colvarproxy.cpp
+++ b/lib/colvars/colvarproxy.cpp
@@ -14,6 +14,11 @@
#include <omp.h>
#endif
+#if defined(NAMD_TCL) || defined(VMDTCL)
+#define COLVARS_TCL
+#include <tcl.h>
+#endif
+
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarscript.h"
@@ -420,8 +425,10 @@ colvarproxy_script::colvarproxy_script()
colvarproxy_script::~colvarproxy_script() {}
-char *colvarproxy_script::script_obj_to_str(unsigned char *obj)
+char const *colvarproxy_script::script_obj_to_str(unsigned char *obj)
{
+ cvm::error("Error: trying to print a script object without a scripting "
+ "language interface.\n", BUG_ERROR);
return reinterpret_cast<char *>(obj);
}
@@ -451,6 +458,140 @@ int colvarproxy_script::run_colvar_gradient_callback(
+colvarproxy_tcl::colvarproxy_tcl()
+{
+ _tcl_interp = NULL;
+}
+
+
+colvarproxy_tcl::~colvarproxy_tcl()
+{
+}
+
+
+void colvarproxy_tcl::init_tcl_pointers()
+{
+ cvm::error("Error: Tcl support is currently unavailable "
+ "outside NAMD or VMD.\n", COLVARS_NOT_IMPLEMENTED);
+}
+
+
+char const *colvarproxy_tcl::tcl_obj_to_str(unsigned char *obj)
+{
+#if defined(COLVARS_TCL)
+ return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
+#else
+ return NULL;
+#endif
+}
+
+
+int colvarproxy_tcl::tcl_run_force_callback()
+{
+#if defined(COLVARS_TCL)
+ Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+ std::string cmd = std::string("calc_colvar_forces ")
+ + cvm::to_str(cvm::step_absolute());
+ int err = Tcl_Eval(tcl_interp, cmd.c_str());
+ if (err != TCL_OK) {
+ cvm::log(std::string("Error while executing calc_colvar_forces:\n"));
+ cvm::error(Tcl_GetStringResult(tcl_interp));
+ return COLVARS_ERROR;
+ }
+ return cvm::get_error();
+#else
+ return COLVARS_NOT_IMPLEMENTED;
+#endif
+}
+
+
+int colvarproxy_tcl::tcl_run_colvar_callback(
+ std::string const &name,
+ std::vector<const colvarvalue *> const &cvc_values,
+ colvarvalue &value)
+{
+#if defined(COLVARS_TCL)
+
+ Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+ size_t i;
+ std::string cmd = std::string("calc_") + name;
+ for (i = 0; i < cvc_values.size(); i++) {
+ cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
+ std::string("}");
+ }
+ int err = Tcl_Eval(tcl_interp, cmd.c_str());
+ const char *result = Tcl_GetStringResult(tcl_interp);
+ if (err != TCL_OK) {
+ return cvm::error(std::string("Error while executing ")
+ + cmd + std::string(":\n") +
+ std::string(Tcl_GetStringResult(tcl_interp)), COLVARS_ERROR);
+ }
+ std::istringstream is(result);
+ if (value.from_simple_string(is.str()) != COLVARS_OK) {
+ cvm::log("Error parsing colvar value from script:");
+ cvm::error(result);
+ return COLVARS_ERROR;
+ }
+ return cvm::get_error();
+
+#else
+
+ return COLVARS_NOT_IMPLEMENTED;
+
+#endif
+}
+
+
+int colvarproxy_tcl::tcl_run_colvar_gradient_callback(
+ std::string const &name,
+ std::vector<const colvarvalue *> const &cvc_values,
+ std::vector<cvm::matrix2d<cvm::real> > &gradient)
+{
+#if defined(COLVARS_TCL)
+
+ Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+ size_t i;
+ std::string cmd = std::string("calc_") + name + "_gradient";
+ for (i = 0; i < cvc_values.size(); i++) {
+ cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
+ std::string("}");
+ }
+ int err = Tcl_Eval(tcl_interp, cmd.c_str());
+ if (err != TCL_OK) {
+ return cvm::error(std::string("Error while executing ")
+ + cmd + std::string(":\n") +
+ std::string(Tcl_GetStringResult(tcl_interp)), COLVARS_ERROR);
+ }
+ Tcl_Obj **list;
+ int n;
+ Tcl_ListObjGetElements(tcl_interp, Tcl_GetObjResult(tcl_interp),
+ &n, &list);
+ if (n != int(gradient.size())) {
+ cvm::error("Error parsing list of gradient values from script: found "
+ + cvm::to_str(n) + " values instead of " +
+ cvm::to_str(gradient.size()));
+ return COLVARS_ERROR;
+ }
+ for (i = 0; i < gradient.size(); i++) {
+ std::istringstream is(Tcl_GetString(list[i]));
+ if (gradient[i].from_simple_string(is.str()) != COLVARS_OK) {
+ cvm::log("Gradient matrix size: " + cvm::to_str(gradient[i].size()));
+ cvm::log("Gradient string: " + cvm::to_str(Tcl_GetString(list[i])));
+ cvm::error("Error parsing gradient value from script", COLVARS_ERROR);
+ return COLVARS_ERROR;
+ }
+ }
+
+ return cvm::get_error();
+
+#else
+
+ return COLVARS_NOT_IMPLEMENTED;
+
+#endif
+}
+
+
colvarproxy_io::colvarproxy_io() {}
@@ -541,6 +682,7 @@ colvarproxy::colvarproxy()
{
colvars = NULL;
b_simulation_running = true;
+ b_delete_requested = false;
}
@@ -556,6 +698,14 @@ int colvarproxy::reset()
}
+int colvarproxy::request_deletion()
+{
+ return cvm::error("Error: \"delete\" command is only available in VMD; "
+ "please use \"reset\" instead.\n",
+ COLVARS_NOT_IMPLEMENTED);
+}
+
+
int colvarproxy::setup()
{
return COLVARS_OK;
@@ -579,13 +729,3 @@ size_t colvarproxy::restart_frequency()
return 0;
}
-
-
-
-
-
-
-
-
-
-
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index e51ddfbe3bfc16a9d9a77264b86912a2cd8aa34f..bf290482704af3204518c0ffae7daf5ad5d88c0b 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -415,7 +415,7 @@ public:
};
-/// Method for scripting language interface (Tcl or Python)
+/// Methods for scripting language interface (Tcl or Python)
class colvarproxy_script {
public:
@@ -427,7 +427,7 @@ public:
virtual ~colvarproxy_script();
/// Convert a script object (Tcl or Python call argument) to a C string
- virtual char *script_obj_to_str(unsigned char *obj);
+ virtual char const *script_obj_to_str(unsigned char *obj);
/// Pointer to the scripting interface object
/// (does not need to be allocated in a new interface)
@@ -454,6 +454,46 @@ public:
};
+/// Methods for using Tcl within Colvars
+class colvarproxy_tcl {
+
+public:
+
+ /// Constructor
+ colvarproxy_tcl();
+
+ /// Destructor
+ virtual ~colvarproxy_tcl();
+
+ /// Is Tcl available? (trigger initialization if needed)
+ int tcl_available();
+
+ /// Tcl implementation of script_obj_to_str()
+ char const *tcl_obj_to_str(unsigned char *obj);
+
+ /// Run a user-defined colvar forces script
+ int tcl_run_force_callback();
+
+ int tcl_run_colvar_callback(
+ std::string const &name,
+ std::vector<const colvarvalue *> const &cvcs,
+ colvarvalue &value);
+
+ int tcl_run_colvar_gradient_callback(
+ std::string const &name,
+ std::vector<const colvarvalue *> const &cvcs,
+ std::vector<cvm::matrix2d<cvm::real> > &gradient);
+
+protected:
+
+ /// Pointer to Tcl interpreter object
+ void *_tcl_interp;
+
+ /// Set Tcl pointers
+ virtual void init_tcl_pointers();
+};
+
+
/// Methods for data input/output
class colvarproxy_io {
@@ -540,6 +580,7 @@ class colvarproxy
public colvarproxy_smp,
public colvarproxy_replicas,
public colvarproxy_script,
+ public colvarproxy_tcl,
public colvarproxy_io
{
@@ -554,6 +595,15 @@ public:
/// Destructor
virtual ~colvarproxy();
+ /// Request deallocation of the module (currently only implemented by VMD)
+ virtual int request_deletion();
+
+ /// Whether deallocation was requested
+ inline bool delete_requested()
+ {
+ return b_delete_requested;
+ }
+
/// \brief Reset proxy state, e.g. requested atoms
virtual int reset();
@@ -591,6 +641,9 @@ protected:
/// Whether a simulation is running (warn against irrecovarable errors)
bool b_simulation_running;
+ /// Whether the entire module should be deallocated by the host engine
+ bool b_delete_requested;
+
};
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index a92a776f8a3abc88700794228bcf4a1e0ad8da51..dc4b8bd07e420d83289b99f458ed2bf935ab210e 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2017-10-20"
+#define COLVARS_VERSION "2018-01-17"
// 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
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 9570acd8327252a99e6394f9bd29d6f6a0bd2f4e..0977496b9e1c252e52c5670f44321796cfafd889 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -74,7 +74,9 @@ int colvarscript::run(int objc, unsigned char *const objv[])
}
if (objc < 2) {
- return exec_command(cv_help, NULL, objc, objv);
+ set_str_result("No commands given: use \"cv help\" "
+ "for a list of commands.");
+ return COLVARSCRIPT_ERROR;
}
std::string const cmd(obj_to_str(objv[1]));
@@ -123,8 +125,7 @@ int colvarscript::run(int objc, unsigned char *const objv[])
if (cmd == "delete") {
// Note: the delete bit may be ignored by some backends
// it is mostly useful in VMD
- colvars->set_error_bits(DELETE_COLVARS);
- return COLVARS_OK;
+ return proxy->request_deletion();
}
if (cmd == "update") {
@@ -272,6 +273,11 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
return COLVARS_OK;
}
+ if (subcmd == "run_ave") {
+ result = (cv->run_ave()).to_simple_string();
+ return COLVARS_OK;
+ }
+
if (subcmd == "width") {
result = cvm::to_str(cv->width, 0, cvm::cv_prec);
return COLVARS_OK;
diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h
index 45ecea867f17c7fdf2b5bbad780ac860287e8836..79f77dad4e8284d4a7c0c771f0e464626a8f828d 100644
--- a/src/USER-COLVARS/colvarproxy_lammps_version.h
+++ b/src/USER-COLVARS/colvarproxy_lammps_version.h
@@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2017-10-20"
+#define COLVARPROXY_VERSION "2017-12-01"
// 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