diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index ad15752107cc006471cc4b3e7a400aa1658f26e2..ef1ed3b54577643d29722a9a176cf5b683d75006 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/Makefile.common b/lib/colvars/Makefile.common
index 5e2f0f268e93944e6b8ce349a21e7754f10c45b7..e3fa4662e638ef1668a21259f7fb5c021c599e9a 100644
--- a/lib/colvars/Makefile.common
+++ b/lib/colvars/Makefile.common
@@ -59,8 +59,7 @@ LEPTON_OBJS = \
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o) $(LEPTON_OBJS)
%.o: %.cpp
- $(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) \
- -Ilepton/include -DLEPTON -c -o $@ $<
+ $(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) -Ilepton/include -DLEPTON -c -o $@ $<
$(COLVARS_LIB): Makefile.deps $(COLVARS_OBJS)
$(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS) $(LEPTON_OBJS)
diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps
index fec7b2be123606b013910c273cdeb9e8b08540e8..e0c768dbc9b042b3ae14731b6319929af9a20705 100644
--- a/lib/colvars/Makefile.deps
+++ b/lib/colvars/Makefile.deps
@@ -14,7 +14,7 @@ $(COLVARS_OBJ_DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
- colvarbias_abf.h colvarbias.h colvargrid.h
+ colvarbias_abf.h colvarbias.h colvargrid.h colvar_UIestimator.h
$(COLVARS_OBJ_DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarbias_alb.h colvar.h colvarparse.h colvardeps.h \
@@ -39,7 +39,8 @@ $(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
- lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h
+ lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
+ colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_histogram.o: colvarbias_histogram.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h \
@@ -198,9 +199,9 @@ $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
- colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
- colvarbias_histogram.h colvarbias_meta.h colvarbias_restraint.h \
- colvarscript.h colvaratoms.h
+ colvarbias.h colvarbias_abf.h colvargrid.h colvar_UIestimator.h \
+ colvarbias_alb.h colvarbias_histogram.h colvarbias_meta.h \
+ colvarbias_restraint.h colvarscript.h colvaratoms.h colvarcomp.h
$(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index d23bd852aaf8be1fe81c41039c9d577af87c7384..d1a398ce4dad6392b4587bffb2508f9ae12b0d19 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1008,6 +1008,8 @@ int colvar::calc()
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
+ colvarproxy *proxy = cvm::main()->proxy;
+
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\", components "+
@@ -1018,14 +1020,18 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
return error_code;
}
- if (cvm::step_relative() > 0) {
- // Total force depends on Jacobian derivative from previous timestep
+ if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
+ // Use Jacobian derivative from previous timestep
error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
}
// atom coordinates are updated by the next line
error_code |= calc_cvc_values(first_cvc, num_cvcs);
error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
+ if (proxy->total_forces_same_step()){
+ // Use Jacobian derivative from this timestep
+ error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
+ }
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\".\n");
@@ -1043,6 +1049,7 @@ int colvar::collect_cvc_data()
if (cvm::step_relative() > 0) {
// Total force depends on Jacobian derivative from previous timestep
+ // collect_cvc_total_forces() uses the previous value of jd
error_code |= collect_cvc_total_forces();
}
error_code |= collect_cvc_values();
@@ -1471,9 +1478,15 @@ cvm::real colvar::update_forces_energy()
// Coupling force is a slow force, to be applied to atomic coords impulse-style
f *= cvm::real(time_step_factor);
- // The total force acting on the extended variable is f_ext
- // This will be used in the next timestep
- ft_reported = f_ext;
+ if (is_enabled(f_cv_subtract_applied_force)) {
+ // Report a "system" force without the biases on this colvar
+ // that is, just the spring force
+ ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+ } else {
+ // The total force acting on the extended variable is f_ext
+ // This will be used in the next timestep
+ ft_reported = f_ext;
+ }
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index dfa9e093a537b797284909e2d7a468571cb8cb98..20dad2771b6c25afef9c4667ca23f24cfe61d46d 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -60,7 +60,10 @@ public:
/// \brief Current actual value (not extended DOF)
colvarvalue const & actual_value() 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;
@@ -96,6 +99,12 @@ public:
{
return cv_features;
}
+ static void delete_features() {
+ for (size_t i=0; i < cv_features.size(); i++) {
+ delete cv_features[i];
+ }
+ cv_features.clear();
+ }
/// Implements possible actions to be carried out
/// when a given feature is enabled
@@ -592,6 +601,10 @@ public:
}
};
+inline cvm::real const & colvar::force_constant() const
+{
+ return ext_force_k;
+}
inline colvarvalue const & colvar::value() const
{
diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fb250b71de01bbb9ffee981209dc14f98bf8b52
--- /dev/null
+++ b/lib/colvars/colvar_UIestimator.h
@@ -0,0 +1,736 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+#ifndef COLVAR_UIESTIMATOR_H
+#define COLVAR_UIESTIMATOR_H
+
+#include <cmath>
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <string>
+
+#include <typeinfo>
+
+// only for colvar module!
+// when integrated into other code, just remove this line and "...cvm::backup_file(...)"
+#include "colvarmodule.h"
+
+namespace UIestimator {
+ const int Y_SIZE = 21; // defines the range of extended CV with respect to a given CV
+ // For example, CV=10, width=1, Y_SIZE=21, then eCV=[0-20], having a size of 21
+ const int HALF_Y_SIZE = 10;
+ const int EXTENDED_X_SIZE = HALF_Y_SIZE;
+ const double EPSILON = 0.000001; // for comparison of float numbers
+
+ class n_matrix { // Stores the distribution matrix of n(x,y)
+
+ public:
+ n_matrix() {}
+ n_matrix(const std::vector<double> & lowerboundary, // lowerboundary of x
+ const std::vector<double> & upperboundary, // upperboundary of
+ const std::vector<double> & width, // width of x
+ const int y_size) { // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
+
+ int i;
+
+ this->lowerboundary = lowerboundary;
+ this->upperboundary = upperboundary;
+ 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(y_size, dimension) + EPSILON);
+
+ // the range of the matrix is [lowerboundary, upperboundary]
+ x_total_size = 1;
+ for (i = 0; i < dimension; i++) {
+ x_size.push_back(int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON));
+ x_total_size *= x_size[i];
+ }
+
+ // initialize the internal matrix
+ matrix.reserve(x_total_size);
+ for (i = 0; i < x_total_size; i++) {
+ matrix.push_back(std::vector<int>(y_total_size, 0));
+ }
+
+ temp.resize(dimension);
+ }
+
+ int inline get_value(const std::vector<double> & x, const std::vector<double> & y) {
+ return matrix[convert_x(x)][convert_y(x, y)];
+ }
+
+ void inline set_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
+ matrix[convert_x(x)][convert_y(x,y)] = value;
+ }
+
+ void inline increase_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
+ matrix[convert_x(x)][convert_y(x,y)] += value;
+ }
+
+ private:
+ std::vector<double> lowerboundary;
+ std::vector<double> upperboundary;
+ std::vector<double> width;
+ int dimension;
+ std::vector<int> x_size; // the size of x in each dimension
+ int x_total_size; // the size of x of the internal matrix
+ int y_size; // the size of y in each dimension
+ int y_total_size; // the size of y of the internal matrix
+
+ std::vector<std::vector<int> > matrix; // the internal matrix
+
+ std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
+
+ int i, j;
+
+ int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
+ for (i = 0; i < dimension; i++) {
+ temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON);
+ }
+
+ int index = 0;
+ for (i = 0; i < dimension; i++) {
+ if (i + 1 < dimension) {
+ int x_temp = 1;
+ for (j = i + 1; j < dimension; j++)
+ x_temp *= x_size[j];
+ index += temp[i] * x_temp;
+ }
+ else
+ index += temp[i];
+ }
+ return index;
+ }
+
+ int convert_y(const std::vector<double> & x, const std::vector<double> & y) { // convert real y value to its interal index
+
+ int i;
+
+ for (i = 0; i < dimension; i++) {
+ temp[i] = round((round(y[i] / width[i] + EPSILON) - round(x[i] / width[i] + EPSILON)) + (y_size - 1) / 2 + EPSILON);
+ }
+
+ int index = 0;
+ for (i = 0; i < dimension; i++) {
+ if (i + 1 < dimension)
+ index += temp[i] * int(pow(y_size, dimension - i - 1) + EPSILON);
+ else
+ index += temp[i];
+ }
+ return index;
+ }
+
+ double round(double r) {
+ return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
+ }
+ };
+
+ // vector, store the sum_x, sum_x_square, count_y
+ template <typename T>
+ class n_vector {
+
+ public:
+ n_vector() {}
+ n_vector(const std::vector<double> & lowerboundary, // lowerboundary of x
+ const std::vector<double> & upperboundary, // upperboundary of
+ const std::vector<double> & width, // width of x
+ const int y_size, // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
+ const T & default_value) { // the default value of T
+
+ this->width = width;
+ this->dimension = lowerboundary.size();
+
+ x_total_size = 1;
+ for (int i = 0; i < dimension; i++) {
+ this->lowerboundary.push_back(lowerboundary[i] - (y_size - 1) / 2 * width[i] - EPSILON);
+ this->upperboundary.push_back(upperboundary[i] + (y_size - 1) / 2 * width[i] + EPSILON);
+
+ x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + EPSILON));
+ x_total_size *= x_size[i];
+ }
+
+ // initialize the internal vector
+ vector.resize(x_total_size, default_value);
+
+ temp.resize(dimension);
+ }
+
+ const T inline get_value(const std::vector<double> & x) {
+ return vector[convert_x(x)];
+ }
+
+ void inline set_value(const std::vector<double> & x, const T value) {
+ vector[convert_x(x)] = value;
+ }
+
+ void inline increase_value(const std::vector<double> & x, const T value) {
+ vector[convert_x(x)] += value;
+ }
+ private:
+ std::vector<double> lowerboundary;
+ std::vector<double> upperboundary;
+ std::vector<double> width;
+ int dimension;
+ std::vector<int> x_size; // the size of x in each dimension
+ int x_total_size; // the size of x of the internal matrix
+
+ std::vector<T> vector; // the internal vector
+
+ std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
+
+ int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
+
+ int i, j;
+
+ for (i = 0; i < dimension; i++) {
+ temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON);
+ }
+
+ int index = 0;
+ for (i = 0; i < dimension; i++) {
+ if (i + 1 < dimension) {
+ int x_temp = 1;
+ for (j = i + 1; j < dimension; j++)
+ x_temp *= x_size[j];
+ index += temp[i] * x_temp;
+ }
+ else
+ index += temp[i];
+ }
+ return index;
+ }
+ };
+
+ class UIestimator { // the implemension of UI estimator
+
+ public:
+ UIestimator() {}
+
+ //called when (re)start an eabf simulation
+ UIestimator(const std::vector<double> & lowerboundary,
+ const std::vector<double> & upperboundary,
+ const std::vector<double> & width,
+ const std::vector<double> & krestr, // force constant in eABF
+ const std::string & output_filename, // the prefix of output files
+ const int output_freq,
+ const bool restart, // whether restart from a .count and a .grad file
+ const std::vector<std::string> & input_filename, // the prefixes of input files
+ const double temperature) {
+
+ // initialize variables
+ this->lowerboundary = lowerboundary;
+ this->upperboundary = upperboundary;
+ this->width = width;
+ this->krestr = krestr;
+ this->output_filename = output_filename;
+ this->output_freq = output_freq;
+ this->restart = restart;
+ this->input_filename = input_filename;
+ this->temperature = temperature;
+
+ int i, j;
+
+ dimension = lowerboundary.size();
+
+ for (i = 0; i < dimension; i++) {
+ sum_x.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
+ sum_x_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
+
+ x_av.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
+ sigma_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
+ }
+
+ count_y = n_vector<int>(lowerboundary, upperboundary, width, Y_SIZE, 0);
+ distribution_x_y = n_matrix(lowerboundary, upperboundary, width, Y_SIZE);
+
+ grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
+ count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
+
+ written = false;
+ written_1D = false;
+
+ if (dimension == 1) {
+ std::vector<double> upperboundary_temp = upperboundary;
+ upperboundary_temp[0] = upperboundary[0] + width[0];
+ oneD_pmf = n_vector<double>(lowerboundary, upperboundary_temp, width, 1, 0.0);
+ }
+
+ if (restart == true) {
+ input_grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
+ input_count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
+
+ // initialize input_Grad and input_count
+ // the loop_flag is a n-dimensional vector, increae from lowerboundary to upperboundary when looping
+ std::vector<double> loop_flag(dimension, 0);
+ for (i = 0; i < dimension; i++) {
+ loop_flag[i] = lowerboundary[i];
+ }
+
+ i = 0;
+ while (i >= 0) {
+ for (j = 0; j < dimension; j++) {
+ input_grad.set_value(loop_flag, std::vector<double>(dimension,0));
+ }
+ input_count.set_value(loop_flag, 0);
+
+ // iterate over any dimensions
+ i = dimension - 1;
+ while (i >= 0) {
+ loop_flag[i] += width[i];
+ if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) {
+ loop_flag[i] = lowerboundary[i];
+ i--;
+ }
+ else
+ break;
+ }
+ }
+ read_inputfiles(input_filename);
+ }
+ }
+
+ ~UIestimator() {}
+
+ // called from MD engine every step
+ bool update(const int step, std::vector<double> x, std::vector<double> y) {
+
+ int i;
+
+ if (step % output_freq == 0) {
+ calc_pmf();
+ write_files();
+ //write_interal_data();
+ }
+
+ for (i = 0; i < dimension; i++) {
+ // for dihedral RC, it is possible that x = 179 and y = -179, should correct it
+ // may have problem, need to fix
+ if (x[i] > 150 && y[i] < -150) {
+ y[i] += 360;
+ }
+ if (x[i] < -150 && y[i] > 150) {
+ y[i] -= 360;
+ }
+
+ if (x[i] < lowerboundary[i] - EXTENDED_X_SIZE * width[i] + EPSILON || x[i] > upperboundary[i] + EXTENDED_X_SIZE * width[i] - EPSILON \
+ || y[i] - x[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - x[i] > HALF_Y_SIZE * width[i] - EPSILON \
+ || y[i] - lowerboundary[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - upperboundary[i] > HALF_Y_SIZE * width[i] - EPSILON)
+ return false;
+ }
+
+ for (i = 0; i < dimension; i++) {
+ sum_x[i].increase_value(y, x[i]);
+ sum_x_square[i].increase_value(y, x[i] * x[i]);
+ }
+ count_y.increase_value(y, 1);
+
+ for (i = 0; i < dimension; i++) {
+ // adapt colvars precision
+ if (x[i] < lowerboundary[i] + EPSILON || x[i] > upperboundary[i] - EPSILON)
+ return false;
+ }
+ distribution_x_y.increase_value(x, y, 1);
+
+ return true;
+ }
+
+ // update the output_filename
+ void update_output_filename(const std::string& filename) {
+ output_filename = filename;
+ }
+
+ private:
+ std::vector<n_vector<double> > sum_x; // the sum of x in each y bin
+ std::vector<n_vector<double> > sum_x_square; // the sum of x in each y bin
+ n_vector<int> count_y; // the distribution of y
+ n_matrix distribution_x_y; // the distribution of <x, y> pair
+
+ int dimension;
+
+ std::vector<double> lowerboundary;
+ std::vector<double> upperboundary;
+ std::vector<double> width;
+ std::vector<double> krestr;
+ std::string output_filename;
+ int output_freq;
+ bool restart;
+ std::vector<std::string> input_filename;
+ double temperature;
+
+ n_vector<std::vector<double> > grad;
+ n_vector<int> count;
+
+ n_vector<double> oneD_pmf;
+
+ n_vector<std::vector<double> > input_grad;
+ n_vector<int> input_count;
+
+ // used in double integration
+ std::vector<n_vector<double> > x_av;
+ std::vector<n_vector<double> > sigma_square;
+
+ bool written;
+ bool written_1D;
+
+ // calculate gradients from the internal variables
+ void calc_pmf() {
+ int norm;
+ int i, j, k;
+
+ std::vector<double> loop_flag(dimension, 0);
+ for (i = 0; i < dimension; i++) {
+ loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
+ }
+
+ i = 0;
+ while (i >= 0) {
+ norm = count_y.get_value(loop_flag) > 0 ? count_y.get_value(loop_flag) : 1;
+ for (j = 0; j < dimension; j++) {
+ x_av[j].set_value(loop_flag, sum_x[j].get_value(loop_flag) / norm);
+ sigma_square[j].set_value(loop_flag, sum_x_square[j].get_value(loop_flag) / norm - x_av[j].get_value(loop_flag) * x_av[j].get_value(loop_flag));
+ }
+
+ // iterate over any dimensions
+ i = dimension - 1;
+ while (i >= 0) {
+ loop_flag[i] += width[i];
+ if (loop_flag[i] > upperboundary[i] + HALF_Y_SIZE * width[i] - width[i] + EPSILON) {
+ loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
+ i--;
+ }
+ else
+ break;
+ }
+ }
+
+ // double integration
+ std::vector<double> av(dimension, 0);
+ std::vector<double> diff_av(dimension, 0);
+
+ std::vector<double> loop_flag_x(dimension, 0);
+ std::vector<double> loop_flag_y(dimension, 0);
+ for (i = 0; i < dimension; i++) {
+ loop_flag_x[i] = lowerboundary[i];
+ loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i];
+ }
+
+ i = 0;
+ while (i >= 0) {
+ norm = 0;
+ for (k = 0; k < dimension; k++) {
+ av[k] = 0;
+ diff_av[k] = 0;
+ loop_flag_y[k] = loop_flag_x[k] - HALF_Y_SIZE * width[k];
+ }
+
+ int j = 0;
+ while (j >= 0) {
+ norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y);
+ for (k = 0; k < dimension; k++) {
+ if (sigma_square[k].get_value(loop_flag_y) > EPSILON || sigma_square[k].get_value(loop_flag_y) < -EPSILON)
+ av[k] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * ( (loop_flag_x[k] + 0.5 * width[k]) - x_av[k].get_value(loop_flag_y)) / sigma_square[k].get_value(loop_flag_y);
+
+ diff_av[k] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * (loop_flag_x[k] - loop_flag_y[k]);
+ }
+
+ // iterate over any dimensions
+ j = dimension - 1;
+ while (j >= 0) {
+ loop_flag_y[j] += width[j];
+ if (loop_flag_y[j] > loop_flag_x[j] + HALF_Y_SIZE * width[j] - width[j] + EPSILON) {
+ loop_flag_y[j] = loop_flag_x[j] - HALF_Y_SIZE * width[j];
+ j--;
+ }
+ else
+ break;
+ }
+ }
+
+ std::vector<double> grad_temp(dimension, 0);
+ for (k = 0; k < dimension; k++) {
+ diff_av[k] /= (norm > 0 ? norm : 1);
+ av[k] = cvm::boltzmann() * temperature * av[k] / (norm > 0 ? norm : 1);
+ grad_temp[k] = av[k] - krestr[k] * diff_av[k];
+ }
+ grad.set_value(loop_flag_x, grad_temp);
+ count.set_value(loop_flag_x, norm);
+
+ // iterate over any dimensions
+ i = dimension - 1;
+ while (i >= 0) {
+ loop_flag_x[i] += width[i];
+ if (loop_flag_x[i] > upperboundary[i] - width[i] + EPSILON) {
+ loop_flag_x[i] = lowerboundary[i];
+ i--;
+ }
+ else
+ break;
+ }
+ }
+ }
+
+
+ // calculate 1D pmf
+ void calc_1D_pmf()
+ {
+ std::vector<double> last_position(1, 0);
+ std::vector<double> position(1, 0);
+
+ double min = 0;
+ double dG = 0;
+ double i;
+
+ oneD_pmf.set_value(lowerboundary, 0);
+ last_position = lowerboundary;
+ for (i = lowerboundary[0] + width[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
+ position[0] = i + EPSILON;
+ if (restart == false || input_count.get_value(last_position) == 0) {
+ dG = oneD_pmf.get_value(last_position) + grad.get_value(last_position)[0] * width[0];
+ }
+ else {
+ dG = oneD_pmf.get_value(last_position) + ((grad.get_value(last_position)[0] * count.get_value(last_position) + input_grad.get_value(last_position)[0] * input_count.get_value(last_position)) / (count.get_value(last_position) + input_count.get_value(last_position))) * width[0];
+ }
+ if (dG < min)
+ min = dG;
+ oneD_pmf.set_value(position, dG);
+ last_position[0] = i + EPSILON;
+ }
+
+ for (i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
+ position[0] = i + EPSILON;
+ oneD_pmf.set_value(position, oneD_pmf.get_value(position) - min);
+ }
+ }
+
+ // write 1D pmf
+ void write_1D_pmf() {
+ std::string pmf_filename = output_filename + ".UI.pmf";
+
+ // only for colvars module!
+ if (written_1D) cvm::backup_file(pmf_filename.c_str());
+
+ std::ostream* ofile_pmf = cvm::proxy->output_stream(pmf_filename.c_str());
+
+ std::vector<double> position(1, 0);
+ for (double i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
+ *ofile_pmf << i << " ";
+ position[0] = i + EPSILON;
+ *ofile_pmf << oneD_pmf.get_value(position) << std::endl;
+ }
+ cvm::proxy->close_output_stream(pmf_filename.c_str());
+
+ written_1D = true;
+ }
+
+ // write heads of the output files
+ void writehead(std::ostream& os) const {
+ os << "# " << dimension << std::endl;
+ for (int i = 0; i < dimension; i++) {
+ os << "# " << lowerboundary[i] << " " << width[i] << " " << int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON) << " " << 0 << std::endl;
+ }
+ os << std::endl;
+ }
+
+ // write interal data, used for testing
+ void write_interal_data() {
+ std::string internal_filename = output_filename + ".UI.internal";
+
+ std::ostream* ofile_internal = cvm::proxy->output_stream(internal_filename.c_str());
+
+ std::vector<double> loop_flag(dimension, 0);
+ for (int i = 0; i < dimension; i++) {
+ loop_flag[i] = lowerboundary[i];
+ }
+
+ int n = 0;
+ while (n >= 0) {
+ for (int j = 0; j < dimension; j++) {
+ *ofile_internal << loop_flag[j] + 0.5 * width[j] << " ";
+ }
+
+ for (int k = 0; k < dimension; k++) {
+ *ofile_internal << grad.get_value(loop_flag)[k] << " ";
+ }
+
+ std::vector<double> ii(dimension,0);
+ for (double i = loop_flag[0] - 10; i < loop_flag[0] + 10 + EPSILON; i+= width[0]) {
+ for (double j = loop_flag[1] - 10; j< loop_flag[1] + 10 + EPSILON; j+=width[1]) {
+ ii[0] = i;
+ ii[1] = j;
+ *ofile_internal << i <<" "<<j<<" "<< distribution_x_y.get_value(loop_flag,ii)<< " ";
+ }
+ }
+ *ofile_internal << std::endl;
+
+ // iterate over any dimensions
+ n = dimension - 1;
+ while (n >= 0) {
+ loop_flag[n] += width[n];
+ if (loop_flag[n] > upperboundary[n] - width[n] + EPSILON) {
+ loop_flag[n] = lowerboundary[n];
+ n--;
+ }
+ else
+ break;
+ }
+ }
+ cvm::proxy->close_output_stream(internal_filename.c_str());
+ }
+
+ // write output files
+ void write_files() {
+ std::string grad_filename = output_filename + ".UI.grad";
+ std::string hist_filename = output_filename + ".UI.hist.grad";
+ std::string count_filename = output_filename + ".UI.count";
+
+ int i, j;
+//
+ // only for colvars module!
+ if (written) cvm::backup_file(grad_filename.c_str());
+ //if (written) cvm::backup_file(hist_filename.c_str());
+ if (written) cvm::backup_file(count_filename.c_str());
+
+ std::ostream* ofile = cvm::proxy->output_stream(grad_filename.c_str());
+ std::ostream* ofile_hist = cvm::proxy->output_stream(hist_filename.c_str(), std::ios::app);
+ std::ostream* ofile_count = cvm::proxy->output_stream(count_filename.c_str());
+
+ writehead(*ofile);
+ writehead(*ofile_hist);
+ writehead(*ofile_count);
+
+ if (dimension == 1) {
+ calc_1D_pmf();
+ write_1D_pmf();
+ }
+
+ std::vector<double> loop_flag(dimension, 0);
+ for (i = 0; i < dimension; i++) {
+ loop_flag[i] = lowerboundary[i];
+ }
+
+ i = 0;
+ while (i >= 0) {
+ for (j = 0; j < dimension; j++) {
+ *ofile << loop_flag[j] + 0.5 * width[j] << " ";
+ *ofile_hist << loop_flag[j] + 0.5 * width[j] << " ";
+ *ofile_count << loop_flag[j] + 0.5 * width[j] << " ";
+ }
+
+ if (restart == false) {
+ for (j = 0; j < dimension; j++) {
+ *ofile << grad.get_value(loop_flag)[j] << " ";
+ *ofile_hist << grad.get_value(loop_flag)[j] << " ";
+ }
+ *ofile << std::endl;
+ *ofile_hist << std::endl;
+ *ofile_count << count.get_value(loop_flag) << " " <<std::endl;
+ }
+ else {
+ double final_grad = 0;
+ for (j = 0; j < dimension; j++) {
+ int total_count_temp = (count.get_value(loop_flag) + input_count.get_value(loop_flag));
+ if (input_count.get_value(loop_flag) == 0)
+ final_grad = grad.get_value(loop_flag)[j];
+ else
+ final_grad = ((grad.get_value(loop_flag)[j] * count.get_value(loop_flag) + input_grad.get_value(loop_flag)[j] * input_count.get_value(loop_flag)) / total_count_temp);
+ *ofile << final_grad << " ";
+ *ofile_hist << final_grad << " ";
+ }
+ *ofile << std::endl;
+ *ofile_hist << std::endl;
+ *ofile_count << (count.get_value(loop_flag) + input_count.get_value(loop_flag)) << " " <<std::endl;
+ }
+
+ // iterate over any dimensions
+ i = dimension - 1;
+ while (i >= 0) {
+ loop_flag[i] += width[i];
+ if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) {
+ loop_flag[i] = lowerboundary[i];
+ i--;
+ *ofile << std::endl;
+ *ofile_hist << std::endl;
+ *ofile_count << std::endl;
+ }
+ else
+ break;
+ }
+ }
+ cvm::proxy->close_output_stream(grad_filename.c_str());
+ cvm::proxy->close_output_stream(hist_filename.c_str());
+ cvm::proxy->close_output_stream(count_filename.c_str());
+
+ written = true;
+ }
+
+ // read input files
+ void read_inputfiles(const std::vector<std::string> input_filename)
+ {
+ char sharp;
+ double nothing;
+ int dimension_temp;
+ int i, j, k, l, m;
+
+ std::vector<double> loop_bin_size(dimension, 0);
+ std::vector<double> position_temp(dimension, 0);
+ std::vector<double> grad_temp(dimension, 0);
+ int count_temp = 0;
+ for (i = 0; i < int(input_filename.size()); i++) {
+ int size = 1 , size_temp = 0;
+
+ std::string count_filename = input_filename[i] + ".UI.count";
+ std::string grad_filename = input_filename[i] + ".UI.grad";
+
+ std::ifstream count_file(count_filename.c_str(), std::ios::in);
+ std::ifstream grad_file(grad_filename.c_str(), std::ios::in);
+
+ count_file >> sharp >> dimension_temp;
+ grad_file >> sharp >> dimension_temp;
+
+ for (j = 0; j < dimension; j++) {
+ count_file >> sharp >> nothing >> nothing >> size_temp >> nothing;
+ grad_file >> sharp >> nothing >> nothing >> nothing >> nothing;
+ size *= size_temp;
+ }
+
+ for (j = 0; j < size; j++) {
+ do {
+ for (k = 0; k < dimension; k++) {
+ count_file >> position_temp[k];
+ grad_file >> nothing;
+ }
+
+ for (l = 0; l < dimension; l++) {
+ grad_file >> grad_temp[l];
+ }
+ count_file >> count_temp;
+ }
+ while (position_temp[i] < lowerboundary[i] - EPSILON || position_temp[i] > upperboundary[i] + EPSILON);
+
+ if (count_temp == 0) {
+ continue;
+ }
+
+ for (m = 0; m < dimension; m++) {
+ grad_temp[m] = (grad_temp[m] * count_temp + input_grad.get_value(position_temp)[m] * input_count.get_value(position_temp)) / (count_temp + input_count.get_value(position_temp));
+ }
+ input_grad.set_value(position_temp, grad_temp);
+ input_count.increase_value(position_temp, count_temp);
+ }
+
+ count_file.close();
+ grad_file.close();
+ }
+ }
+ };
+};
+
+#endif
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 9b4a922e3fbb2322455c17adca9194d2032ab6a1..d2a0f0a807d2727500d96e6a19a56eaa0af23f39 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -817,6 +817,18 @@ int cvm::atom_group::create_sorted_ids(void)
}
+int cvm::atom_group::overlap(const atom_group &g1, const atom_group &g2){
+ for (cvm::atom_const_iter ai1 = g1.begin(); ai1 != g1.end(); ai1++) {
+ for (cvm::atom_const_iter ai2 = g2.begin(); ai2 != g2.end(); ai2++) {
+ if (ai1->id == ai2->id) {
+ return (ai1->id + 1); // 1-based index to allow boolean usage
+ }
+ }
+ }
+ return 0;
+}
+
+
void cvm::atom_group::center_ref_pos()
{
ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 6113fb38a938b37a56d54825c51d43252d5b54db..71c587e23084e516c6b66358e97327c8d404822f 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -214,6 +214,12 @@ public:
{
return ag_features;
}
+ static void delete_features() {
+ for (size_t i=0; i < ag_features.size(); i++) {
+ delete ag_features[i];
+ }
+ ag_features.clear();
+ }
protected:
@@ -280,6 +286,10 @@ public:
/// Allocates and populates the sorted list of atom ids
int create_sorted_ids(void);
+ /// Detect whether two groups share atoms
+ /// If yes, returns 1-based number of a common atom; else, returns 0
+ static int overlap(const atom_group &g1, const atom_group &g2);
+
/// \brief When updating atomic coordinates, translate them to align with the
/// center of mass of the reference coordinates
bool b_center;
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 636727ca39b4f90b8b9686a64b51845857843c13..902ea4605d2f72aec3d79a546463ef2321d5180d 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -10,6 +10,7 @@
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias.h"
+#include "colvargrid.h"
colvarbias::colvarbias(char const *key)
@@ -31,12 +32,14 @@ int colvarbias::init(std::string const &conf)
{
colvarparse::init(conf);
+ size_t i = 0;
+
if (name.size() == 0) {
// first initialization
cvm::log("Initializing a new \""+bias_type+"\" instance.\n");
- rank = cvm::num_biases_type(bias_type);
+ rank = cvm::main()->num_biases_type(bias_type);
get_keyval(conf, "name", name, bias_type+cvm::to_str(rank));
{
@@ -62,7 +65,7 @@ int colvarbias::init(std::string const &conf)
INPUT_ERROR);
return INPUT_ERROR;
}
- for (size_t i = 0; i < colvar_names.size(); i++) {
+ for (i = 0; i < colvar_names.size(); i++) {
add_colvar(colvar_names[i]);
}
}
@@ -148,6 +151,13 @@ int colvarbias::clear()
}
+int colvarbias::clear_state_data()
+{
+ // no mutable content to delete for base class
+ return COLVARS_OK;
+}
+
+
int colvarbias::add_colvar(std::string const &cv_name)
{
if (colvar *cv = cvm::colvar_by_name(cv_name)) {
@@ -164,6 +174,8 @@ int colvarbias::add_colvar(std::string const &cv_name)
colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
colvar_forces.back().reset();
+ previous_colvar_forces.push_back(colvar_forces.back());
+
cv->biases.push_back(this); // add back-reference to this bias to colvar
if (is_enabled(f_cvb_apply_force)) {
@@ -204,7 +216,8 @@ int colvarbias::update()
void colvarbias::communicate_forces()
{
- for (size_t i = 0; i < num_variables(); i++) {
+ size_t i = 0;
+ for (i = 0; i < num_variables(); i++) {
if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+
variables(i)->name+"\".\n");
@@ -216,6 +229,9 @@ void colvarbias::communicate_forces()
// aware of this bias' time_step_factor
variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]);
}
+ for (i = 0; i < num_variables(); i++) {
+ previous_colvar_forces[i] = colvar_forces[i];
+ }
}
@@ -389,6 +405,248 @@ std::ostream & colvarbias::write_traj(std::ostream &os)
return os;
}
+
+
+colvarbias_ti::colvarbias_ti(char const *key)
+ : colvarbias(key)
+{
+ provide(f_cvb_calc_ti_samples);
+ ti_avg_forces = NULL;
+ ti_count = NULL;
+}
+
+
+colvarbias_ti::~colvarbias_ti()
+{
+ colvarbias_ti::clear_state_data();
+}
+
+
+int colvarbias_ti::clear_state_data()
+{
+ if (ti_avg_forces != NULL) {
+ delete ti_avg_forces;
+ ti_avg_forces = NULL;
+ }
+ if (ti_count != NULL) {
+ delete ti_count;
+ ti_count = NULL;
+ }
+ return COLVARS_OK;
+}
+
+
+int colvarbias_ti::init(std::string const &conf)
+{
+ int error_code = COLVARS_OK;
+
+ get_keyval_feature(this, conf, "writeTISamples",
+ f_cvb_write_ti_samples,
+ is_enabled(f_cvb_write_ti_samples));
+
+ get_keyval_feature(this, conf, "writeTIPMF",
+ f_cvb_write_ti_pmf,
+ is_enabled(f_cvb_write_ti_pmf));
+
+ if ((num_variables() > 1) && is_enabled(f_cvb_write_ti_pmf)) {
+ return cvm::error("Error: only 1-dimensional PMFs can be written "
+ "on the fly.\n"
+ "Consider using writeTISamples instead and "
+ "post-processing the sampled free-energy gradients.\n",
+ COLVARS_NOT_IMPLEMENTED);
+ } else {
+ error_code |= init_grids();
+ }
+
+ if (is_enabled(f_cvb_write_ti_pmf)) {
+ enable(f_cvb_write_ti_samples);
+ }
+
+ if (is_enabled(f_cvb_calc_ti_samples)) {
+ std::vector<std::string> const time_biases =
+ cvm::main()->time_dependent_biases();
+ if (time_biases.size() > 0) {
+ if ((time_biases.size() > 1) || (time_biases[0] != this->name)) {
+ for (size_t i = 0; i < num_variables(); i++) {
+ if (! variables(i)->is_enabled(f_cv_subtract_applied_force)) {
+ return cvm::error("Error: cannot collect TI samples while other "
+ "time-dependent biases are active and not all "
+ "variables have subtractAppliedForces on.\n",
+ INPUT_ERROR);
+ }
+ }
+ }
+ }
+ }
+
+ return error_code;
+}
+
+
+int colvarbias_ti::init_grids()
+{
+ if (is_enabled(f_cvb_calc_ti_samples)) {
+ if (ti_avg_forces == NULL) {
+ ti_bin.resize(num_variables());
+ ti_system_forces.resize(num_variables());
+ for (size_t icv = 0; icv < num_variables(); icv++) {
+ ti_system_forces[icv].type(variables(icv)->value());
+ ti_system_forces[icv].is_derivative();
+ ti_system_forces[icv].reset();
+ }
+ ti_avg_forces = new colvar_grid_gradient(colvars);
+ ti_count = new colvar_grid_count(colvars);
+ ti_avg_forces->samples = ti_count;
+ ti_count->has_parent_data = true;
+ }
+ }
+
+ return COLVARS_OK;
+}
+
+
+int colvarbias_ti::update()
+{
+ return update_system_forces(NULL);
+}
+
+
+int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
+ *subtract_forces)
+{
+ if (! is_enabled(f_cvb_calc_ti_samples)) {
+ return COLVARS_OK;
+ }
+
+ has_data = true;
+
+ if (cvm::debug()) {
+ cvm::log("Updating system forces for bias "+this->name+"\n");
+ }
+
+ // Collect total colvar forces from the previous step
+ size_t i;
+ if (cvm::step_relative() > 0) {
+ if (ti_avg_forces->index_ok(ti_bin)) {
+ for (i = 0; i < num_variables(); i++) {
+ if (variables(i)->is_enabled(f_cv_subtract_applied_force)) {
+ // this colvar is already subtracting all applied forces
+ ti_system_forces[i] = variables(i)->total_force();
+ } else {
+ ti_system_forces[i] = variables(i)->total_force() -
+ ((subtract_forces != NULL) ?
+ (*subtract_forces)[i] : previous_colvar_forces[i]);
+ }
+ }
+ ti_avg_forces->acc_value(ti_bin, ti_system_forces);
+ }
+ }
+
+ // Set the index to be used in the next iteration, when total forces come in
+ for (i = 0; i < num_variables(); i++) {
+ ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
+ }
+
+ return COLVARS_OK;
+}
+
+
+std::string const colvarbias_ti::get_state_params() const
+{
+ return std::string("");
+}
+
+
+int colvarbias_ti::set_state_params(std::string const &state_conf)
+{
+ return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias_ti::write_state_data(std::ostream &os)
+{
+ if (! is_enabled(f_cvb_calc_ti_samples)) {
+ return os;
+ }
+ os << "\nhistogram\n";
+ ti_count->write_raw(os);
+ os << "\nsystem_forces\n";
+ ti_avg_forces->write_raw(os);
+ return os;
+}
+
+
+std::istream & colvarbias_ti::read_state_data(std::istream &is)
+{
+ if (! is_enabled(f_cvb_calc_ti_samples)) {
+ return is;
+ }
+ if (cvm::debug()) {
+ cvm::log("Reading state data for the TI estimator.\n");
+ }
+ if (! read_state_data_key(is, "histogram")) {
+ return is;
+ }
+ if (! ti_count->read_raw(is)) {
+ return is;
+ }
+ if (! read_state_data_key(is, "system_forces")) {
+ return is;
+ }
+ if (! ti_avg_forces->read_raw(is)) {
+ return is;
+ }
+ if (cvm::debug()) {
+ cvm::log("Done reading state data for the TI estimator.\n");
+ }
+ return is;
+}
+
+
+int colvarbias_ti::write_output_files()
+{
+ if (!has_data) {
+ // nothing to write
+ return COLVARS_OK;
+ }
+
+ std::string const ti_output_prefix = cvm::output_prefix()+"."+this->name;
+
+ std::ostream *os = NULL;
+
+ if (is_enabled(f_cvb_write_ti_samples)) {
+ std::string const ti_count_file_name(ti_output_prefix+".ti.count");
+ os = cvm::proxy->output_stream(ti_count_file_name);
+ if (os) {
+ ti_count->write_multicol(*os);
+ cvm::proxy->close_output_stream(ti_count_file_name);
+ }
+
+ std::string const ti_grad_file_name(ti_output_prefix+".ti.grad");
+ os = cvm::proxy->output_stream(ti_grad_file_name);
+ if (os) {
+ ti_avg_forces->write_multicol(*os);
+ cvm::proxy->close_output_stream(ti_grad_file_name);
+ }
+ }
+
+ if (is_enabled(f_cvb_write_ti_pmf)) {
+ std::string const pmf_file_name(ti_output_prefix+".ti.pmf");
+ cvm::log("Writing TI PMF to file \""+pmf_file_name+"\".\n");
+ os = cvm::proxy->output_stream(pmf_file_name);
+ if (os) {
+ // get the FE gradient
+ ti_avg_forces->multiply_constant(-1.0);
+ ti_avg_forces->write_1D_integral(*os);
+ ti_avg_forces->multiply_constant(-1.0);
+ cvm::proxy->close_output_stream(pmf_file_name);
+ }
+ }
+
+ return COLVARS_OK;
+}
+
+
// Static members
std::vector<colvardeps::feature *> colvarbias::cvb_features;
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index a147cd3210486cc74d80f5e0c0e5f859e0ea1b70..083b9d73036da056a2c76d3cfb8cef600049c1f3 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -109,6 +109,9 @@ public:
/// \brief Delete everything
virtual int clear();
+ /// \brief Delete only the allocatable data (save memory)
+ virtual int clear_state_data();
+
/// Destructor
virtual ~colvarbias();
@@ -183,6 +186,12 @@ public:
{
return cvb_features;
}
+ static void delete_features() {
+ for (size_t i=0; i < cvb_features.size(); i++) {
+ delete cvb_features[i];
+ }
+ cvb_features.clear();
+ }
protected:
@@ -194,6 +203,9 @@ protected:
/// \brief Current forces from this bias to the variables
std::vector<colvarvalue> colvar_forces;
+ /// \brief Forces last applied by this bias to the variables
+ std::vector<colvarvalue> previous_colvar_forces;
+
/// \brief Current energy of this bias (colvar_forces should be obtained by deriving this)
cvm::real bias_energy;
@@ -209,4 +221,48 @@ protected:
};
+
+class colvar_grid_gradient;
+class colvar_grid_count;
+
+/// \brief Base class for unconstrained thermodynamic-integration FE estimator
+class colvarbias_ti : public virtual colvarbias {
+public:
+
+ colvarbias_ti(char const *key);
+ virtual ~colvarbias_ti();
+
+ virtual int clear_state_data();
+
+ virtual int init(std::string const &conf);
+ virtual int init_grids();
+ virtual int update();
+
+ /// Subtract applied forces (either last forces or argument) from the total
+ /// forces
+ virtual int update_system_forces(std::vector<colvarvalue> const
+ *subtract_forces);
+
+ virtual std::string const get_state_params() const;
+ virtual int set_state_params(std::string const &state_conf);
+ virtual std::ostream & write_state_data(std::ostream &os);
+ virtual std::istream & read_state_data(std::istream &is);
+ virtual int write_output_files();
+
+protected:
+
+ /// \brief Forces exerted from the system to the associated variables
+ std::vector<colvarvalue> ti_system_forces;
+
+ /// Averaged system forces
+ colvar_grid_gradient *ti_avg_forces;
+
+ /// Histogram of sampled data
+ colvar_grid_count *ti_count;
+
+ /// Because total forces may be from the last simulation step,
+ /// store the index of the variables then
+ std::vector<int> ti_bin;
+};
+
#endif
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index a96fc21d644e750d1a574596e4c3f0b99c87930f..5f5fe149784362e95807c2cfc7e6ef207ca90726 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -14,6 +14,8 @@
colvarbias_abf::colvarbias_abf(char const *key)
: colvarbias(key),
+ b_UI_estimator(false),
+ b_CZAR_estimator(false),
system_force(NULL),
gradients(NULL),
samples(NULL),
@@ -159,6 +161,7 @@ int colvarbias_abf::init(std::string const &conf)
// Data for eABF 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);
@@ -187,8 +190,38 @@ int colvarbias_abf::init(std::string const &conf)
read_gradients_samples();
}
- cvm::log("Finished ABF setup.\n");
+ // if extendedLangrangian is on, then call UI estimator
+ if (b_extended) {
+ get_keyval(conf, "UIestimator", b_UI_estimator, false);
+
+ if (b_UI_estimator) {
+ std::vector<double> UI_lowerboundary;
+ std::vector<double> UI_upperboundary;
+ std::vector<double> UI_width;
+ std::vector<double> UI_krestr;
+
+ bool UI_restart = (input_prefix.size() > 0);
+
+ for (size_t i = 0; i < colvars.size(); i++)
+ {
+ UI_lowerboundary.push_back(colvars[i]->lower_boundary);
+ UI_upperboundary.push_back(colvars[i]->upper_boundary);
+ UI_width.push_back(colvars[i]->width);
+ UI_krestr.push_back(colvars[i]->force_constant());
+ }
+ eabf_UI = UIestimator::UIestimator(UI_lowerboundary,
+ UI_upperboundary,
+ UI_width,
+ UI_krestr, // force constant in eABF
+ output_prefix, // the prefix of output files
+ cvm::restart_out_freq,
+ UI_restart, // whether restart from a .count and a .grad file
+ input_prefix, // the prefixes of input files
+ cvm::temperature());
+ }
+ }
+ cvm::log("Finished ABF setup.\n");
return COLVARS_OK;
}
@@ -332,7 +365,7 @@ int colvarbias_abf::update()
}
// update the output prefix; TODO: move later to setup_output() function
- if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
+ if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
// This is the only bias computing PMFs
output_prefix = cvm::output_prefix();
} else {
@@ -364,6 +397,20 @@ int colvarbias_abf::update()
cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
}
+ // 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++)
+ {
+ x[i] = colvars[i]->actual_value();
+ y[i] = colvars[i]->value();
+ }
+ eabf_UI.update_output_filename(output_prefix);
+ eabf_UI.update(cvm::step_absolute(), x, y);
+ }
+
return COLVARS_OK;
}
@@ -479,8 +526,8 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
cvm::proxy->close_output_stream(pmf_out_name);
}
- if (z_gradients) {
- // Write eABF-related quantities
+ if (b_CZAR_estimator) {
+ // Write eABF CZAR-related quantities
std::string z_samples_out_name = prefix + ".zcount";
@@ -588,7 +635,7 @@ void colvarbias_abf::read_gradients_samples()
is.close();
}
- if (z_gradients) {
+ if (b_CZAR_estimator) {
// Read eABF z-averaged data for CZAR
cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);
@@ -621,7 +668,7 @@ std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
os << "\ngradient\n";
gradients->write_raw(os, 8);
- if (z_gradients) {
+ if (b_CZAR_estimator) {
os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
os << "\nz_samples\n";
z_samples->write_raw(os, 8);
@@ -655,7 +702,7 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
return is;
}
- if (z_gradients) {
+ if (b_CZAR_estimator) {
if (! read_state_data_key(is, "z_samples")) {
return is;
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index 41a5475fa70087ca709e4ca2be9f0f118874a311..1defe72268c487e82a536d2f5e441e1067897272 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -17,6 +17,7 @@
#include "colvarbias.h"
#include "colvargrid.h"
+#include "colvar_UIestimator.h"
typedef cvm::real* gradient_t;
@@ -50,6 +51,12 @@ private:
/// Write CZAR output file for stratified eABF (.zgrad)
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?
bool cap_force;
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index 66806fc9fca01d57a596344ff2fd62a120272442..b0d154dfc9ebf591f49b15ba517f4bcf7ba10317 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -33,7 +33,7 @@
colvarbias_meta::colvarbias_meta(char const *key)
- : colvarbias(key)
+ : colvarbias(key), colvarbias_ti(key)
{
new_hills_begin = hills.end();
hills_traj_os = NULL;
@@ -44,6 +44,7 @@ colvarbias_meta::colvarbias_meta(char const *key)
int colvarbias_meta::init(std::string const &conf)
{
colvarbias::init(conf);
+ colvarbias_ti::init(conf);
enable(f_cvb_calc_pmf);
@@ -104,7 +105,7 @@ int colvarbias_meta::init(std::string const &conf)
get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
- "please use \"keepFreeEnergyFile\" instead.");
+ "please use \"keepFreeEnergyFiles\" instead.");
}
get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
@@ -230,15 +231,7 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
colvarbias_meta::~colvarbias_meta()
{
- if (hills_energy) {
- delete hills_energy;
- hills_energy = NULL;
- }
-
- if (hills_energy_gradients) {
- delete hills_energy_gradients;
- hills_energy_gradients = NULL;
- }
+ colvarbias_meta::clear_state_data();
if (replica_hills_os) {
cvm::proxy->close_output_stream(replica_hills_file);
@@ -250,13 +243,31 @@ colvarbias_meta::~colvarbias_meta()
hills_traj_os = NULL;
}
- if(target_dist) {
+ if (target_dist) {
delete target_dist;
target_dist = NULL;
}
}
+int colvarbias_meta::clear_state_data()
+{
+ if (hills_energy) {
+ delete hills_energy;
+ hills_energy = NULL;
+ }
+
+ if (hills_energy_gradients) {
+ delete hills_energy_gradients;
+ hills_energy_gradients = NULL;
+ }
+
+ hills.clear();
+ hills_off_grid.clear();
+
+ return COLVARS_OK;
+}
+
// **********************************************************************
// Hill management member functions
@@ -336,6 +347,9 @@ int colvarbias_meta::update()
// update base class
error_code |= colvarbias::update();
+ // update the TI estimator (if defined)
+ error_code |= colvarbias_ti::update();
+
// update grid definition, if needed
error_code |= update_grid_params();
// add new biasing energy/forces
@@ -1000,6 +1014,10 @@ void colvarbias_meta::update_replicas_registry()
(replicas.back())->hills_energy = new colvar_grid_scalar(colvars);
(replicas.back())->hills_energy_gradients = new colvar_grid_gradient(colvars);
}
+ if (is_enabled(f_cvb_calc_ti_samples)) {
+ (replicas.back())->enable(f_cvb_calc_ti_samples);
+ (replicas.back())->colvarbias_ti::init_grids();
+ }
}
}
} else {
@@ -1374,6 +1392,8 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
}
}
+ colvarbias_ti::read_state_data(is);
+
if (cvm::debug())
cvm::log("colvarbias_meta::read_restart() done\n");
@@ -1474,7 +1494,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
int colvarbias_meta::setup_output()
{
output_prefix = cvm::output_prefix();
- if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) {
+ if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) {
// if this is not the only free energy integrator, append
// this bias's name, to distinguish it from the output of the other
// biases producing a .pmf file
@@ -1631,6 +1651,7 @@ std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
}
}
+ colvarbias_ti::write_state_data(os);
return os;
}
@@ -1651,6 +1672,7 @@ int colvarbias_meta::write_state_to_replicas()
int colvarbias_meta::write_output_files()
{
+ colvarbias_ti::write_output_files();
if (dump_fes) {
write_pmf();
}
diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h
index 249f7342bc4b982b3747c56a11f9362b63914a41..78b2d35d41bb581d9b6da2297cf5886285415de5 100644
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@@ -19,7 +19,10 @@
#include "colvargrid.h"
/// Metadynamics bias (implementation of \link colvarbias \endlink)
-class colvarbias_meta : public colvarbias {
+class colvarbias_meta
+ : public virtual colvarbias,
+ public virtual colvarbias_ti
+{
public:
@@ -35,10 +38,13 @@ public:
Communication comm;
colvarbias_meta(char const *key);
+ virtual ~colvarbias_meta();
+
virtual int init(std::string const &conf);
virtual int init_well_tempered_params(std::string const &conf);
virtual int init_ebmeta_params(std::string const &conf);
- virtual ~colvarbias_meta();
+
+ virtual int clear_state_data();
virtual int update();
virtual int update_grid_params();
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index 70beca29fa79eb25b328c0afa1d3c8cefdf8226a..23534f56eb43c1cd0a2e257ed8581449026d0c8f 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -14,7 +14,7 @@
colvarbias_restraint::colvarbias_restraint(char const *key)
- : colvarbias(key)
+ : colvarbias(key), colvarbias_ti(key)
{
}
@@ -24,6 +24,8 @@ int colvarbias_restraint::init(std::string const &conf)
colvarbias::init(conf);
enable(f_cvb_apply_force);
+ colvarbias_ti::init(conf);
+
if (cvm::debug())
cvm::log("Initializing a new restraint bias.\n");
@@ -86,7 +88,7 @@ std::ostream & colvarbias_restraint::write_traj(std::ostream &os)
colvarbias_restraint_centers::colvarbias_restraint_centers(char const *key)
- : colvarbias(key), colvarbias_restraint(key)
+ : colvarbias(key), colvarbias_ti(key), colvarbias_restraint(key)
{
}
@@ -145,7 +147,7 @@ int colvarbias_restraint_centers::change_configuration(std::string const &conf)
colvarbias_restraint_k::colvarbias_restraint_k(char const *key)
- : colvarbias(key), colvarbias_restraint(key)
+ : colvarbias(key), colvarbias_ti(key), colvarbias_restraint(key)
{
force_k = -1.0;
}
@@ -237,6 +239,7 @@ int colvarbias_restraint_moving::set_state_params(std::string const &conf)
colvarbias_restraint_centers_moving::colvarbias_restraint_centers_moving(char const *key)
: colvarbias(key),
+ colvarbias_ti(key),
colvarbias_restraint(key),
colvarbias_restraint_centers(key),
colvarbias_restraint_moving(key)
@@ -284,14 +287,17 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf)
target_centers[i],
0.5);
}
+
+ get_keyval(conf, "outputAccumulatedWork", b_output_acc_work,
+ b_output_acc_work); // TODO this conflicts with stages
+
} else {
target_centers.clear();
- return COLVARS_OK;
}
+ // Output restraint centers even when they do not change; some NAMD REUS
+ // scripts expect this behavior
get_keyval(conf, "outputCenters", b_output_centers, b_output_centers);
- get_keyval(conf, "outputAccumulatedWork", b_output_acc_work,
- b_output_acc_work); // TODO this conflicts with stages
return COLVARS_OK;
}
@@ -475,6 +481,7 @@ std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key)
: colvarbias(key),
+ colvarbias_ti(key),
colvarbias_restraint(key),
colvarbias_restraint_k(key),
colvarbias_restraint_moving(key)
@@ -712,6 +719,7 @@ std::ostream & colvarbias_restraint::write_state(std::ostream &os)
colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
: colvarbias(key),
+ colvarbias_ti(key),
colvarbias_restraint(key),
colvarbias_restraint_centers(key),
colvarbias_restraint_moving(key),
@@ -743,17 +751,22 @@ int colvarbias_restraint_harmonic::init(std::string const &conf)
int colvarbias_restraint_harmonic::update()
{
+ int error_code = COLVARS_OK;
+
+ // update the TI estimator (if defined)
+ error_code |= colvarbias_ti::update();
+
// update parameters (centers or force constant)
- colvarbias_restraint_centers_moving::update();
- colvarbias_restraint_k_moving::update();
+ error_code |= colvarbias_restraint_centers_moving::update();
+ error_code |= colvarbias_restraint_k_moving::update();
// update restraint energy and forces
- colvarbias_restraint::update();
+ error_code |= colvarbias_restraint::update();
// update accumulated work using the current forces
- colvarbias_restraint_centers_moving::update_acc_work();
+ error_code |= colvarbias_restraint_centers_moving::update_acc_work();
- return COLVARS_OK;
+ return error_code;
}
@@ -798,6 +811,18 @@ int colvarbias_restraint_harmonic::set_state_params(std::string const &conf)
}
+std::ostream & colvarbias_restraint_harmonic::write_state_data(std::ostream &os)
+{
+ return colvarbias_ti::write_state_data(os);
+}
+
+
+std::istream & colvarbias_restraint_harmonic::read_state_data(std::istream &is)
+{
+ return colvarbias_ti::read_state_data(is);
+}
+
+
std::ostream & colvarbias_restraint_harmonic::write_traj_label(std::ostream &os)
{
colvarbias_restraint::write_traj_label(os);
@@ -845,6 +870,7 @@ cvm::real colvarbias_restraint_harmonic::energy_difference(std::string const &co
colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char const *key)
: colvarbias(key),
+ colvarbias_ti(key),
colvarbias_restraint(key),
colvarbias_restraint_k(key),
colvarbias_restraint_moving(key),
@@ -967,11 +993,15 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
int colvarbias_restraint_harmonic_walls::update()
{
- colvarbias_restraint_k_moving::update();
+ int error_code = COLVARS_OK;
- colvarbias_restraint::update();
+ error_code |= colvarbias_ti::update();
- return COLVARS_OK;
+ error_code |= colvarbias_restraint_k_moving::update();
+
+ error_code |= colvarbias_restraint::update();
+
+ return error_code;
}
@@ -1065,6 +1095,18 @@ int colvarbias_restraint_harmonic_walls::set_state_params(std::string const &con
}
+std::ostream & colvarbias_restraint_harmonic_walls::write_state_data(std::ostream &os)
+{
+ return colvarbias_ti::write_state_data(os);
+}
+
+
+std::istream & colvarbias_restraint_harmonic_walls::read_state_data(std::istream &is)
+{
+ return colvarbias_ti::read_state_data(is);
+}
+
+
std::ostream & colvarbias_restraint_harmonic_walls::write_traj_label(std::ostream &os)
{
colvarbias_restraint::write_traj_label(os);
@@ -1084,6 +1126,7 @@ std::ostream & colvarbias_restraint_harmonic_walls::write_traj(std::ostream &os)
colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
: colvarbias(key),
+ colvarbias_ti(key),
colvarbias_restraint(key),
colvarbias_restraint_centers(key),
colvarbias_restraint_moving(key),
@@ -1120,17 +1163,22 @@ int colvarbias_restraint_linear::init(std::string const &conf)
int colvarbias_restraint_linear::update()
{
+ int error_code = COLVARS_OK;
+
+ // update the TI estimator (if defined)
+ error_code |= colvarbias_ti::update();
+
// update parameters (centers or force constant)
- colvarbias_restraint_centers_moving::update();
- colvarbias_restraint_k_moving::update();
+ error_code |= colvarbias_restraint_centers_moving::update();
+ error_code |= colvarbias_restraint_k_moving::update();
// update restraint energy and forces
- colvarbias_restraint::update();
+ error_code |= colvarbias_restraint::update();
// update accumulated work using the current forces
- colvarbias_restraint_centers_moving::update_acc_work();
+ error_code |= colvarbias_restraint_centers_moving::update_acc_work();
- return COLVARS_OK;
+ return error_code;
}
@@ -1196,6 +1244,18 @@ int colvarbias_restraint_linear::set_state_params(std::string const &conf)
}
+std::ostream & colvarbias_restraint_linear::write_state_data(std::ostream &os)
+{
+ return colvarbias_ti::write_state_data(os);
+}
+
+
+std::istream & colvarbias_restraint_linear::read_state_data(std::istream &is)
+{
+ return colvarbias_ti::read_state_data(is);
+}
+
+
std::ostream & colvarbias_restraint_linear::write_traj_label(std::ostream &os)
{
colvarbias_restraint::write_traj_label(os);
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index 8c3a1537fc881cbca3dc8ef66fcf3caeb4d49c90..b10649cab112f8a10bfedaff0c10823cc7a43f04 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -16,7 +16,8 @@
/// see derived classes for specific types
/// (implementation of \link colvarbias \endlink)
class colvarbias_restraint
- : public virtual colvarbias
+ : public virtual colvarbias,
+ public virtual colvarbias_ti
{
public:
@@ -95,7 +96,7 @@ protected:
/// Options to change the restraint configuration over time (shared between centers and k moving)
class colvarbias_restraint_moving
- : public virtual colvarparse {
+ : public virtual colvarparse, public virtual colvardeps {
public:
colvarbias_restraint_moving(char const *key);
@@ -226,6 +227,8 @@ public:
virtual int update();
virtual std::string const get_state_params() const;
virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_state_data(std::ostream &os);
+ virtual std::istream & read_state_data(std::istream &os);
virtual std::ostream & write_traj_label(std::ostream &os);
virtual std::ostream & write_traj(std::ostream &os);
virtual int change_configuration(std::string const &conf);
@@ -252,6 +255,8 @@ public:
virtual void communicate_forces();
virtual std::string const get_state_params() const;
virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_state_data(std::ostream &os);
+ virtual std::istream & read_state_data(std::istream &os);
virtual std::ostream & write_traj_label(std::ostream &os);
virtual std::ostream & write_traj(std::ostream &os);
@@ -292,6 +297,8 @@ public:
virtual std::string const get_state_params() const;
virtual int set_state_params(std::string const &conf);
+ virtual std::ostream & write_state_data(std::ostream &os);
+ virtual std::istream & read_state_data(std::istream &os);
virtual std::ostream & write_traj_label(std::ostream &os);
virtual std::ostream & write_traj(std::ostream &os);
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 3c1ec2495c9c5af29c5549b9c9ddaa355d2f1502..b94d798be9f0a3cd49abc27bb8cdaf41f88e8c87 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -140,7 +140,12 @@ public:
{
return cvc_features;
}
-
+ static void delete_features() {
+ for (size_t i=0; i < cvc_features.size(); i++) {
+ delete cvc_features[i];
+ }
+ cvc_features.clear();
+ }
/// \brief Obtain data needed for the calculation for the backend
virtual void read_data();
diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp
index 369d489e279c04d5051d8fef83ca0d16e5cf67ed..edfea96795b77d7dbd13fd174594d2b9ff007fc5 100644
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@@ -87,6 +87,12 @@ colvar::coordnum::coordnum(std::string const &conf)
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
+ if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
+ cvm::error("Error: group1 and group2 share a common atom (number: " +
+ cvm::to_str(atom_number) + ")\n");
+ return;
+ }
+
if (group1->b_dummy) {
cvm::error("Error: only group2 is allowed to be a dummy atom\n");
return;
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 18d154515a3cfee7ade375083f70ab932da50067..ce8055843f93d389117f38926e81d97c66bc7954 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -1066,8 +1066,9 @@ void colvar::rmsd::calc_force_invgrads()
void colvar::rmsd::calc_Jacobian_derivative()
{
// divergence of the rotated coordinates (including only derivatives of the rotation matrix)
- cvm::real divergence = 0.0;
+ cvm::real rotation_term = 0.0;
+ // The rotation term only applies is coordinates are rotated
if (atoms->b_rotate) {
// gradient of the rotation matrix
@@ -1104,7 +1105,7 @@ void colvar::rmsd::calc_Jacobian_derivative()
for (size_t alpha = 0; alpha < 3; alpha++) {
for (size_t beta = 0; beta < 3; beta++) {
- divergence += grad_rot_mat[beta][alpha][alpha] * y[beta];
+ rotation_term += grad_rot_mat[beta][alpha][alpha] * y[beta];
// Note: equation was derived for inverse rotation (see colvars paper)
// so here the matrix is transposed
// (eq would give divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];)
@@ -1112,7 +1113,13 @@ void colvar::rmsd::calc_Jacobian_derivative()
}
}
}
- jd.real_value = x.real_value > 0.0 ? (3.0 * atoms->size() - 4.0 - divergence) / x.real_value : 0.0;
+
+ // The translation term only applies is coordinates are centered
+ cvm::real translation_term = atoms->b_center ? 3.0 : 0.0;
+
+ jd.real_value = x.real_value > 0.0 ?
+ (3.0 * atoms->size() - 1.0 - translation_term - rotation_term) / x.real_value :
+ 0.0;
}
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 8f241a6255953f2ae0d4815a050daa12c077f866..ac906e7be755fc382d1cb383730522da38860092 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -413,15 +413,27 @@ void colvardeps::init_cvb_requires() {
init_feature(f_cvb_apply_force, "apply force", f_type_user);
f_req_children(f_cvb_apply_force, f_cv_gradient);
- init_feature(f_cvb_get_total_force, "obtain total force");
+ 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_history_dependent, "history-dependent", f_type_static);
+ init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
+
init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
f_req_children(f_cvb_scalar_variables, f_cv_scalar);
init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
+
+ init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic);
+ f_req_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
+ f_req_children(f_cvb_calc_ti_samples, f_cv_grid);
+
+ init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user);
+ f_req_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
+
+ init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user);
+ f_req_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
}
// Initialize feature_states for each instance
@@ -431,6 +443,9 @@ void colvardeps::init_cvb_requires() {
// Most features are available, so we set them so
// and list exceptions below
}
+
+ // only compute TI samples when deriving from colvarbias_ti
+ feature_states[f_cvb_calc_ti_samples].available = false;
}
@@ -504,9 +519,6 @@ void colvardeps::init_cv_requires() {
init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
- // There is no well-defined way to implement f_cv_subtract_applied_force
- // in the case of extended-Lagrangian colvars
- f_req_exclude(f_cv_subtract_applied_force, f_cv_extended_Lagrangian);
init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
f_req_self(f_cv_lower_boundary, f_cv_scalar);
@@ -514,7 +526,7 @@ void colvardeps::init_cv_requires() {
init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
f_req_self(f_cv_upper_boundary, f_cv_scalar);
- init_feature(f_cv_grid, "grid", f_type_user);
+ init_feature(f_cv_grid, "grid", f_type_dynamic);
f_req_self(f_cv_grid, f_cv_lower_boundary);
f_req_self(f_cv_grid, f_cv_upper_boundary);
@@ -693,7 +705,6 @@ void colvardeps::print_state() {
}
-
void colvardeps::add_child(colvardeps *child) {
children.push_back(child);
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index dfb10d00e421f7635fe4bce5ec88b5f45cfc39eb..bd892fbca8cef746bb4dfad818b1e4577507a735 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -180,8 +180,6 @@ public:
protected:
-
-
/// Parse a keyword and enable a feature accordingly
bool get_keyval_feature(colvarparse *cvp,
std::string const &conf, char const *key,
@@ -229,10 +227,18 @@ public:
f_cvb_get_total_force,
/// \brief depends on simulation history
f_cvb_history_dependent,
+ /// \brief depends on time
+ f_cvb_time_dependent,
/// \brief requires scalar colvars
f_cvb_scalar_variables,
/// \brief whether this bias will compute a PMF
f_cvb_calc_pmf,
+ /// \brief whether this bias will compute TI samples
+ f_cvb_calc_ti_samples,
+ /// \brief whether this bias will write TI samples
+ f_cvb_write_ti_samples,
+ /// \brief whether this bias should write the TI PMF
+ f_cvb_write_ti_pmf,
f_cvb_ntot
};
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 6f06cb1066fffbc4ad8982a643264f496264657a..a01104dba819c2745a9d63d434b8e6832c9d85b7 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1403,6 +1403,15 @@ public:
/// Constructor from a vector of colvars
colvar_grid_gradient(std::vector<colvar *> &colvars);
+ /// \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++) {
+ data[address(ix) + imult] += values[imult].real_value;
+ }
+ if (samples)
+ samples->incr_count(ix);
+ }
+
/// \brief Accumulate the gradient
inline void acc_grad(std::vector<int> const &ix, cvm::real const *grads) {
for (size_t imult = 0; imult < mult; imult++) {
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 780dc28afaebd4c6db1788c0495e73f8b7904d5b..f293aac63f2ce433410898098127ab856cb89b84 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -22,7 +22,7 @@
#include "colvarbias_restraint.h"
#include "colvarscript.h"
#include "colvaratoms.h"
-
+#include "colvarcomp.h"
colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
@@ -274,9 +274,9 @@ int colvarmodule::parse_global_params(std::string const &conf)
parse->get_keyval(conf, "colvarsRestartFrequency",
restart_out_freq, restart_out_freq);
- // if this is true when initializing, it means
- // we are continuing after a reset(): default to true
- parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append);
+ // Deprecate append flag
+ parse->get_keyval(conf, "colvarsTrajAppend",
+ cv_traj_append, cv_traj_append, colvarparse::parse_silent);
parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);
@@ -409,22 +409,12 @@ int colvarmodule::parse_biases(std::string const &conf)
cvm::decrease_depth();
}
- size_t i;
-
- size_t n_hist_dep_biases = 0;
- std::vector<std::string> hist_dep_biases_names;
- for (i = 0; i < biases.size(); i++) {
- if (biases[i]->is_enabled(colvardeps::f_cvb_apply_force) &&
- biases[i]->is_enabled(colvardeps::f_cvb_history_dependent)) {
- n_hist_dep_biases++;
- hist_dep_biases_names.push_back(biases[i]->name);
- }
- }
- if (n_hist_dep_biases > 1) {
- cvm::log("WARNING: there are "+cvm::to_str(n_hist_dep_biases)+
- " history-dependent biases with non-zero force parameters:\n"+
- cvm::to_str(hist_dep_biases_names)+"\n"+
- "Please make sure that their forces do not counteract each other.\n");
+ std::vector<std::string> const time_biases = time_dependent_biases();
+ if (time_biases.size() > 1) {
+ cvm::log("WARNING: there are "+cvm::to_str(time_biases.size())+
+ " time-dependent biases with non-zero force parameters:\n"+
+ cvm::to_str(time_biases)+"\n"+
+ "Please ensure that their forces do not counteract each other.\n");
}
if (biases.size() || use_scripted_forces) {
@@ -441,7 +431,7 @@ int colvarmodule::parse_biases(std::string const &conf)
}
-int colvarmodule::num_biases_feature(int feature_id)
+int colvarmodule::num_biases_feature(int feature_id) const
{
colvarmodule *cv = cvm::main();
size_t n = 0;
@@ -456,7 +446,7 @@ int colvarmodule::num_biases_feature(int feature_id)
}
-int colvarmodule::num_biases_type(std::string const &type)
+int colvarmodule::num_biases_type(std::string const &type) const
{
colvarmodule *cv = cvm::main();
size_t n = 0;
@@ -471,6 +461,22 @@ int colvarmodule::num_biases_type(std::string const &type)
}
+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++) {
+ 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) ||
+ biases[i]->is_enabled(colvardeps::f_cvb_time_dependent))) {
+ biases_names.push_back(biases[i]->name);
+ }
+ }
+ return biases_names;
+}
+
+
int colvarmodule::catch_input_errors(int result)
{
if (result != COLVARS_OK || get_error()) {
@@ -673,8 +679,15 @@ int colvarmodule::calc()
}
// write restart files, if needed
- if (restart_out_freq && restart_out_name.size()) {
- error_code |= write_restart_files();
+ if (restart_out_freq && (cvm::step_relative() > 0) &&
+ ((cvm::step_absolute() % restart_out_freq) == 0) ) {
+ if (restart_out_name.size()) {
+ // Write restart file, if different from main output
+ error_code |= write_restart_file(restart_out_name);
+ } else {
+ error_code |= write_restart_file(output_prefix()+".colvars.state");
+ }
+ write_output_files();
}
return error_code;
@@ -916,21 +929,16 @@ int colvarmodule::calc_scripted_forces()
}
-int colvarmodule::write_restart_files()
+int colvarmodule::write_restart_file(std::string const &out_name)
{
- if ( (cvm::step_relative() > 0) &&
- ((cvm::step_absolute() % restart_out_freq) == 0) ) {
- cvm::log("Writing the state file \""+
- restart_out_name+"\".\n");
- proxy->backup_file(restart_out_name);
- std::ostream *restart_out_os = proxy->output_stream(restart_out_name);
- if (!restart_out_os) return cvm::get_error();
- if (!write_restart(*restart_out_os)) {
- return cvm::error("Error: in writing restart file.\n", FILE_ERROR);
- }
- proxy->close_output_stream(restart_out_name);
+ cvm::log("Saving collective variables state to \""+out_name+"\".\n");
+ proxy->backup_file(out_name);
+ std::ostream *restart_out_os = proxy->output_stream(out_name);
+ if (!restart_out_os) return cvm::get_error();
+ if (!write_restart(*restart_out_os)) {
+ return cvm::error("Error: in writing restart file.\n", FILE_ERROR);
}
-
+ proxy->close_output_stream(out_name);
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@@ -1011,7 +1019,15 @@ colvarmodule::~colvarmodule()
{
if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) ||
(proxy->smp_thread_id() == 0)) {
+
+ // Delete contents of static arrays
+ colvarbias::delete_features();
+ colvar::delete_features();
+ colvar::cvc::delete_features();
+ atom_group::delete_features();
+
reset();
+
delete parse;
parse = NULL;
proxy = NULL;
@@ -1261,7 +1277,7 @@ continue the previous simulation.\n\n");
to:\n\
\""+ proxy->input_prefix()+".colvars.state\"\n");
output_prefix() = output_prefix()+".tmp";
- write_output_files();
+ write_restart_file(output_prefix()+".colvars.state");
cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
}
@@ -1277,24 +1293,13 @@ int colvarmodule::backup_file(char const *filename)
int colvarmodule::write_output_files()
{
- // if this is a simulation run (i.e. not a postprocessing), output data
- // must be written to be able to restart the simulation
- std::string const out_name =
- (output_prefix().size() ?
- std::string(output_prefix()+".colvars.state") :
- std::string("colvars.state"));
- cvm::log("Saving collective variables state to \""+out_name+"\".\n");
-
- std::ostream * os = proxy->output_stream(out_name);
- os->setf(std::ios::scientific, std::ios::floatfield);
- this->write_restart(*os);
- proxy->close_output_stream(out_name);
+ int error_code = COLVARS_OK;
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
- (*cvi)->write_output_files();
+ error_code |= (*cvi)->write_output_files();
}
cvm::decrease_depth();
@@ -1302,8 +1307,8 @@ int colvarmodule::write_output_files()
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
- (*bi)->write_output_files();
- (*bi)->write_state_to_replicas();
+ error_code |= (*bi)->write_output_files();
+ error_code |= (*bi)->write_state_to_replicas();
}
cvm::decrease_depth();
@@ -1403,15 +1408,12 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
cvi != colvars.end();
cvi++) {
(*cvi)->write_restart(os);
- error_code |= (*cvi)->write_output_files();
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_state(os);
- error_code |= (*bi)->write_state_to_replicas();
- error_code |= (*bi)->write_output_files();
}
cvm::decrease_depth();
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 0f6efd14c4ea1b7574304c7bc16f067b735b5076..3f13ae88cf07cb438db71f5716f8bc465090c0ae 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -293,10 +293,13 @@ private:
public:
/// Return how many biases have this feature enabled
- static int num_biases_feature(int feature_id);
+ int num_biases_feature(int feature_id) const;
/// Return how many biases are defined with this type
- static int num_biases_type(std::string const &type);
+ int num_biases_type(std::string const &type) const;
+
+ /// Return the names of time-dependent biases with forces enabled
+ std::vector<std::string> const time_dependent_biases() const;
private:
/// Useful wrapper to interrupt parsing if any error occurs
@@ -334,9 +337,9 @@ public:
/// Write all trajectory files
int write_traj_files();
- /// Write all restart files
- int write_restart_files();
- /// Write all FINAL output files
+ /// Write a state file useful to resume the simulation
+ int write_restart_file(std::string const &out_name);
+ /// Write all other output files
int write_output_files();
/// Backup a file before writing it
static int backup_file(char const *filename);
@@ -580,7 +583,7 @@ public:
/// from static functions in the colvarmodule class
static colvarproxy *proxy;
- /// \brief Accessor for the above
+ /// \brief Access the one instance of the Colvars module
static colvarmodule *main();
};
diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp
index fa24091d5233da18a1af36f0dffceca3f89a26fb..8160144c6bab1e11d36f4d1b7712e75bce90e6c8 100644
--- a/lib/colvars/colvarproxy.cpp
+++ b/lib/colvars/colvarproxy.cpp
@@ -10,6 +10,10 @@
#include <sstream>
#include <string.h>
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarscript.h"
@@ -40,6 +44,12 @@ bool colvarproxy_system::total_forces_enabled() const
}
+bool colvarproxy_system::total_forces_same_step() const
+{
+ return false;
+}
+
+
cvm::real colvarproxy_system::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
@@ -204,7 +214,13 @@ void colvarproxy_atom_groups::clear_atom_group(int index)
colvarproxy_smp::colvarproxy_smp()
{
- b_smp_active = true;
+ b_smp_active = true; // May be disabled by user option
+ omp_lock_state = NULL;
+#if defined(_OPENMP)
+ if (smp_thread_id() == 0) {
+ omp_init_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state));
+ }
+#endif
}
@@ -213,60 +229,143 @@ colvarproxy_smp::~colvarproxy_smp() {}
int colvarproxy_smp::smp_enabled()
{
+#if defined(_OPENMP)
+ if (b_smp_active) {
+ return COLVARS_OK;
+ }
+ return COLVARS_ERROR;
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
int colvarproxy_smp::smp_colvars_loop()
{
+#if defined(_OPENMP)
+ colvarmodule *cv = cvm::main();
+ colvarproxy *proxy = cv->proxy;
+#pragma omp parallel for
+ for (size_t i = 0; i < cv->variables_active_smp()->size(); i++) {
+ colvar *x = (*(cv->variables_active_smp()))[i];
+ int x_item = (*(cv->variables_active_smp_items()))[i];
+ if (cvm::debug()) {
+ cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+
+ cvm::to_str(proxy->smp_num_threads())+
+ "]: calc_colvars_items_smp(), i = "+cvm::to_str(i)+", cv = "+
+ x->name+", cvc = "+cvm::to_str(x_item)+"\n");
+ }
+ x->calc_cvcs(x_item, 1);
+ }
+ return cvm::get_error();
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
int colvarproxy_smp::smp_biases_loop()
{
+#if defined(_OPENMP)
+ colvarmodule *cv = cvm::main();
+#pragma omp parallel
+ {
+#pragma omp for
+ for (size_t i = 0; i < cv->biases_active()->size(); i++) {
+ colvarbias *b = (*(cv->biases_active()))[i];
+ if (cvm::debug()) {
+ cvm::log("Calculating bias \""+b->name+"\" on thread "+
+ cvm::to_str(smp_thread_id())+"\n");
+ }
+ b->update();
+ }
+ }
+ return cvm::get_error();
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
int colvarproxy_smp::smp_biases_script_loop()
{
+#if defined(_OPENMP)
+ colvarmodule *cv = cvm::main();
+#pragma omp parallel
+ {
+#pragma omp single nowait
+ {
+ cv->calc_scripted_forces();
+ }
+#pragma omp for
+ for (size_t i = 0; i < cv->biases_active()->size(); i++) {
+ colvarbias *b = (*(cv->biases_active()))[i];
+ if (cvm::debug()) {
+ cvm::log("Calculating bias \""+b->name+"\" on thread "+
+ cvm::to_str(smp_thread_id())+"\n");
+ }
+ b->update();
+ }
+ }
+ return cvm::get_error();
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
+
+
int colvarproxy_smp::smp_thread_id()
{
+#if defined(_OPENMP)
+ return omp_get_thread_num();
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
int colvarproxy_smp::smp_num_threads()
{
+#if defined(_OPENMP)
+ return omp_get_max_threads();
+#else
return COLVARS_NOT_IMPLEMENTED;
+#endif
}
int colvarproxy_smp::smp_lock()
{
+#if defined(_OPENMP)
+ omp_set_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state));
+#endif
return COLVARS_OK;
}
int colvarproxy_smp::smp_trylock()
{
+#if defined(_OPENMP)
+ return omp_test_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state)) ?
+ COLVARS_OK : COLVARS_ERROR;
+#else
return COLVARS_OK;
+#endif
}
int colvarproxy_smp::smp_unlock()
{
+#if defined(_OPENMP)
+ omp_unset_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state));
+#endif
return COLVARS_OK;
}
-
colvarproxy_replicas::colvarproxy_replicas() {}
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 95d13cd7e05c15ba209333f4a6550e11d29b5977..e51ddfbe3bfc16a9d9a77264b86912a2cd8aa34f 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -80,6 +80,9 @@ public:
/// Are total forces being used?
virtual bool total_forces_enabled() const;
+
+ /// Are total forces from the current step available?
+ virtual bool total_forces_same_step() const;
};
@@ -372,6 +375,11 @@ public:
/// Release the lock
virtual int smp_unlock();
+
+protected:
+
+ /// Lock state for OpenMP
+ void *omp_lock_state;
};
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index 312c0fd1a0c86e832ccbac15605867ddd101c824..9dfeb17898295ee80f33832738eac202fe660041 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2017-08-06"
+#define COLVARS_VERSION "2017-10-11"
// 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 89302a16a2ead1c29bbc501ec6d8e2aa9d043cb8..9570acd8327252a99e6394f9bd29d6f6a0bd2f4e 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -11,7 +11,10 @@
#include <stdlib.h>
#include <string.h>
+#define COLVARSCRIPT_CPP
#include "colvarscript.h"
+#undef COLVARSCRIPT_CPP
+
#include "colvarproxy.h"
#include "colvardeps.h"
@@ -21,6 +24,11 @@ colvarscript::colvarscript(colvarproxy *p)
colvars(p->colvars),
proxy_error(0)
{
+ comm_help.resize(colvarscript::cv_n_commands);
+ comm_fns.resize(colvarscript::cv_n_commands);
+#define COLVARSCRIPT_INIT_FN
+#include "colvarscript.h"
+#undef COLVARSCRIPT_INIT_FN
}
@@ -66,8 +74,7 @@ int colvarscript::run(int objc, unsigned char *const objv[])
}
if (objc < 2) {
- result = help_string();
- return COLVARS_OK;
+ return exec_command(cv_help, NULL, objc, objv);
}
std::string const cmd(obj_to_str(objv[1]));
@@ -167,17 +174,7 @@ int colvarscript::run(int objc, unsigned char *const objv[])
/// Parse config from string
if (cmd == "config") {
- if (objc < 3) {
- result = "Missing arguments\n" + help_string();
- return COLVARSCRIPT_ERROR;
- }
- std::string const conf(obj_to_str(objv[2]));
- if (colvars->read_config_string(conf) == COLVARS_OK) {
- return COLVARS_OK;
- } else {
- result = "Error parsing configuration string";
- return COLVARSCRIPT_ERROR;
- }
+ return exec_command(cv_config, NULL, objc, objv);
}
/// Load an input state file
@@ -204,6 +201,8 @@ int colvarscript::run(int objc, unsigned char *const objv[])
proxy->output_prefix() = obj_to_str(objv[2]);
int error = 0;
error |= colvars->setup_output();
+ error |= colvars->write_restart_file(colvars->output_prefix()+
+ ".colvars.state");
error |= colvars->write_output_files();
return error ? COLVARSCRIPT_ERROR : COLVARS_OK;
}
@@ -255,6 +254,10 @@ int colvarscript::run(int objc, unsigned char *const objv[])
}
}
+ if (cmd == "help") {
+ return exec_command(cv_help, NULL, objc, objv);
+ }
+
result = "Syntax error\n" + help_string();
return COLVARSCRIPT_ERROR;
}
@@ -295,7 +298,9 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
- colvars->write_traj_label(*(colvars->cv_traj_os));
+ if (colvars->cv_traj_os != NULL) {
+ colvars->write_traj_label(*(colvars->cv_traj_os));
+ }
return COLVARS_OK;
}
@@ -374,7 +379,6 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
int colvarscript::proc_bias(colvarbias *b, int objc, unsigned char *const objv[]) {
- std::string const key(obj_to_str(objv[0]));
std::string const subcmd(obj_to_str(objv[2]));
if (subcmd == "energy") {
@@ -425,7 +429,9 @@ int colvarscript::proc_bias(colvarbias *b, int objc, unsigned char *const objv[]
// the bias destructor takes care of the cleanup at cvm level
delete b;
// TODO this could be done by the destructors
- colvars->write_traj_label(*(colvars->cv_traj_os));
+ if (colvars->cv_traj_os != NULL) {
+ colvars->write_traj_label(*(colvars->cv_traj_os));
+ }
return COLVARS_OK;
}
@@ -528,7 +534,7 @@ int colvarscript::proc_features(colvardeps *obj,
}
-std::string colvarscript::help_string()
+std::string colvarscript::help_string() const
{
std::string buf;
buf = "Usage: cv <subcommand> [args...]\n\
@@ -538,7 +544,7 @@ Managing the Colvars module:\n\
config <string> -- read configuration from the given string\n\
reset -- delete all internal configuration\n\
delete -- delete this Colvars module instance\n\
- version -- return version of colvars code\n\
+ version -- return version of Colvars code\n\
\n\
Input and output:\n\
list -- return a list of all variables\n\
@@ -564,6 +570,8 @@ Accessing collective variables:\n\
colvar <name> type -- return the type of colvar <name>\n\
colvar <name> delete -- delete colvar <name>\n\
colvar <name> addforce <F> -- apply given force on colvar <name>\n\
+ colvar <name> getappliedforce -- return applied force of colvar <name>\n\
+ colvar <name> gettotalforce -- return total force of colvar <name>\n\
colvar <name> getconfig -- return config string of colvar <name>\n\
colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\
colvar <name> get <f> -- get the value of the colvar feature <f>\n\
diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h
index 94d451809cdd9feb853dbf5462342a52140d907a..39cd08934059f25e0fe598dc0a915ec79036631a 100644
--- a/lib/colvars/colvarscript.h
+++ b/lib/colvars/colvarscript.h
@@ -8,21 +8,27 @@
// Colvars repository at GitHub.
#ifndef COLVARSCRIPT_H
-#define COLVARSCRIPT_H
+//#define COLVARSCRIPT_H // Delay definition until later
#include <string>
+#include <vector>
+#include <map>
+
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias.h"
#include "colvarproxy.h"
+
// Only these error values are part of the scripting interface
#define COLVARSCRIPT_ERROR -1
#define COLVARSCRIPT_OK 0
+
class colvarscript {
private:
+
colvarproxy *proxy;
colvarmodule *colvars;
@@ -35,16 +41,93 @@ public:
colvarscript(colvarproxy * p);
inline ~colvarscript() {}
- /// If an error is caught by the proxy through fatal_error(), this is set to COLVARSCRIPT_ERROR
+ /// If an error is caught by the proxy through fatal_error(), this is set to
+ /// COLVARSCRIPT_ERROR
int proxy_error;
- /// If an error is returned by one of the methods, it should set this to the error message
+ /// If an error is returned by one of the methods, it should set this to the
+ /// error message
std::string result;
/// Run script command with given positional arguments (objects)
int run(int objc, unsigned char *const objv[]);
+ /// Set the return value of the script command to the given string
+ inline void set_str_result(std::string const &s)
+ {
+ result = s;
+ }
+
+ /// Build and return a short help
+ std::string help_string(void) const;
+
+ /// Use scripting language to get the string representation of an object
+ inline char const *obj_to_str(unsigned char *const obj)
+ {
+ return cvm::proxy->script_obj_to_str(obj);
+ }
+
+ enum command {
+ cv_help,
+ cv_version,
+ cv_config,
+ cv_configfile,
+ cv_reset,
+ cv_delete,
+ cv_list,
+ cv_list_biases,
+ cv_load,
+ cv_save,
+ cv_update,
+ cv_addenergy,
+ cv_getenergy,
+ cv_printframe,
+ cv_printframelabels,
+ cv_frame,
+ cv_colvar,
+ cv_colvar_value,
+ cv_colvar_update,
+ cv_colvar_type,
+ cv_colvar_delete,
+ cv_colvar_addforce,
+ cv_colvar_getappliedforce,
+ cv_colvar_gettotalforce,
+ cv_colvar_cvcflags,
+ cv_colvar_getconfig,
+ cv_colvar_get,
+ cv_colvar_set,
+ cv_bias,
+ cv_bias_energy,
+ cv_bias_update,
+ cv_bias_delete,
+ cv_bias_getconfig,
+ cv_bias_get,
+ cv_bias_set,
+ cv_n_commands
+ };
+
+ /// Execute a script command
+ inline int exec_command(command c,
+ void *pobj,
+ int objc, unsigned char * const *objv)
+ {
+ return (*(comm_fns[c]))(pobj, objc, objv);
+ }
+
+ /// Get help for a command (TODO reformat for each language?)
+ inline std::string command_help(colvarscript::command c) const
+ {
+ return comm_help[c];
+ }
+
+ /// Clear all object results
+ inline void clear_results()
+ {
+ result.clear();
+ }
+
private:
+
/// Run subcommands on colvar
int proc_colvar(colvar *cv, int argc, unsigned char *const argv[]);
@@ -55,17 +138,146 @@ private:
int proc_features(colvardeps *obj,
int argc, unsigned char *const argv[]);
- /// Build and return a short help
- std::string help_string(void);
+ /// Internal identifiers of command strings
+ std::map<std::string, command> comm_str_map;
-public:
+ /// Help strings for each command
+ std::vector<std::string> comm_help;
- inline char const *obj_to_str(unsigned char *const obj)
- {
- return cvm::proxy->script_obj_to_str(obj);
- }
+ /// Number of arguments for each command
+ std::vector<size_t> comm_n_args;
+
+ /// Arguments for each command
+ std::vector< std::vector<std::string> > comm_args;
+
+ /// Implementations of each command
+ std::vector<int (*)(void *, int, unsigned char * const *)> comm_fns;
};
+/// Get a pointer to the main colvarscript object
+inline static colvarscript *colvarscript_obj()
+{
+ return cvm::main()->proxy->script;
+}
+
+/// Get a pointer to the colvar object pointed to by pobj
+inline static colvar *colvar_obj(void *pobj)
+{
+ return reinterpret_cast<colvar *>(pobj);
+}
+
+/// Get a pointer to the colvarbias object pointed to by pobj
+inline static colvarbias *colvarbias_obj(void *pobj)
+{
+ return reinterpret_cast<colvarbias *>(pobj);
+}
+
+
+#define CVSCRIPT_COMM_FNAME(COMM) cvscript_ ## COMM
+
+#define CVSCRIPT_COMM_PROTO(COMM) \
+ int CVSCRIPT_COMM_FNAME(COMM)(void *, int, unsigned char *const *);
+
+#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
+ CVSCRIPT_COMM_PROTO(COMM)
+
+#undef COLVARSCRIPT_H
+#endif // #ifndef COLVARSCRIPT_H
+
+
+#ifdef COLVARSCRIPT_CPP
+#define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
+ int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \
+ int objc, unsigned char *const objv[]) \
+ { \
+ colvarscript *script = colvarscript_obj(); \
+ script->clear_results(); \
+ if (objc < 2+N_ARGS_MIN) /* "cv" and "COMM" are 1st and 2nd */ { \
+ script->set_str_result("Missing arguments\n" + \
+ script->command_help(colvarscript::COMM)); \
+ return COLVARSCRIPT_ERROR; \
+ } \
+ if (objc > 2+N_ARGS_MAX) { \
+ script->set_str_result("Too many arguments\n" + \
+ script->command_help(colvarscript::COMM)); \
+ return COLVARSCRIPT_ERROR; \
+ } \
+ FN_BODY; \
+ }
+#undef CVSCRIPT
+#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
+ CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)
+#endif // #ifdef COLVARSCRIPT_CPP
+
+
+#ifdef COLVARSCRIPT_INIT_FN
+#define CVSCRIPT_COMM_INIT(COMM,HELP,ARGS) { \
+ comm_str_map[#COMM] = COMM; \
+ comm_help[COMM] = HELP; \
+ comm_fns[COMM] = &(CVSCRIPT_COMM_FNAME(COMM)); \
+ }
+#undef CVSCRIPT
+#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
+ CVSCRIPT_COMM_INIT(COMM,HELP,ARGS)
#endif
+
+
+#if !defined(COLVARSCRIPT_H) || defined(COLVARSCRIPT_INIT_FN)
+#define COLVARSCRIPT_H
+
+#ifndef COLVARSCRIPT_INIT_FN
+#ifdef __cplusplus
+extern "C" {
+#endif
+#endif
+
+ // Add optional arguments for command-specific help?
+ CVSCRIPT(cv_help,
+ "Print the help message",
+ 0, 0,
+ {},
+ script->set_str_result(script->help_string());
+ return COLVARS_OK;
+ )
+
+ CVSCRIPT(cv_config,
+ "Read configuration from the given string",
+ 1, 1,
+ { "conf (str) - Configuration string" },
+ std::string const conf(script->obj_to_str(objv[2]));
+ if (cvm::main()->read_config_string(conf) == COLVARS_OK) {
+ return COLVARS_OK;
+ }
+ script->set_str_result("Error parsing configuration string");
+ return COLVARSCRIPT_ERROR;
+ )
+
+ CVSCRIPT(cv_addenergy,
+ "Add an energy to the MD engine",
+ 1, 1,
+ { "E (float) - Amount of energy to add" },
+ cvm::main()->total_bias_energy +=
+ strtod(script->obj_to_str(objv[2]), NULL);
+ return COLVARS_OK;
+ )
+
+ CVSCRIPT(cv_getenergy,
+ "Get the current Colvars energy",
+ 1, 1,
+ { "E (float) - Store the energy in this variable" },
+ double *energy = reinterpret_cast<double *>(objv[2]);
+ *energy = cvm::main()->total_bias_energy;
+ return COLVARS_OK;
+ )
+
+#ifndef COLVARSCRIPT_INIT_FN
+#ifdef __cplusplus
+} // extern "C"
+#endif
+#endif
+
+#undef CVSCRIPT
+
+#endif // #ifndef COLVARSCRIPT_H
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index fe3160eb4b3238ef036c2ea56d6b3954d73cb121..97257d18ad66d1da99b0b74586190e735687850b 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -705,7 +705,7 @@ public:
{
std::stringstream stream(s);
size_t i = 0;
- while ((stream >> data[i]) && (i < data.size())) {
+ while ((i < data.size()) && (stream >> data[i])) {
i++;
}
if (i < data.size()) {
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index 17dff305673456b1312ef2be36393f6c8e54281e..ad2366f3d4963162e86f2ba10053aedf1b671455 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -120,12 +120,6 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
if (restart_output_prefix_str.rfind(".*") != std::string::npos)
restart_output_prefix_str.erase(restart_output_prefix_str.rfind(".*"),2);
-#if defined(_OPENMP)
- if (smp_thread_id() == 0) {
- omp_init_lock(&smp_lock_state);
- }
-#endif
-
// initialize multi-replica support, if available
if (replica_enabled()) {
MPI_Comm_rank(inter_comm, &inter_me);
@@ -239,6 +233,9 @@ void colvarproxy_lammps::serialize_status(std::string &rst)
std::ostringstream os;
colvars->write_restart(os);
rst = os.str();
+
+ // TODO separate this as its own function?
+ colvars->write_output_files();
}
// set status from string
@@ -331,89 +328,6 @@ int colvarproxy_lammps::backup_file(char const *filename)
}
-#if defined(_OPENMP)
-
-
-// SMP support
-
-int colvarproxy_lammps::smp_enabled()
-{
- if (b_smp_active) {
- return COLVARS_OK;
- }
- return COLVARS_ERROR;
-}
-
-
-int colvarproxy_lammps::smp_colvars_loop()
-{
- colvarmodule *cv = this->colvars;
- colvarproxy_lammps *proxy = (colvarproxy_lammps *) cv->proxy;
-#pragma omp parallel for
- for (size_t i = 0; i < cv->variables_active_smp()->size(); i++) {
- colvar *x = (*(cv->variables_active_smp()))[i];
- int x_item = (*(cv->variables_active_smp_items()))[i];
- if (cvm::debug()) {
- cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
- "]: calc_colvars_items_smp(), i = "+cvm::to_str(i)+", cv = "+
- x->name+", cvc = "+cvm::to_str(x_item)+"\n");
- }
- x->calc_cvcs(x_item, 1);
- }
- return cvm::get_error();
-}
-
-
-int colvarproxy_lammps::smp_biases_loop()
-{
- colvarmodule *cv = this->colvars;
-#pragma omp parallel for
- for (size_t i = 0; i < cv->biases_active()->size(); i++) {
- colvarbias *b = (*(cv->biases_active()))[i];
- if (cvm::debug()) {
- cvm::log("Calculating bias \""+b->name+"\" on thread "+
- cvm::to_str(smp_thread_id())+"\n");
- }
- b->update();
- }
- return cvm::get_error();
-}
-
-
-int colvarproxy_lammps::smp_thread_id()
-{
- return omp_get_thread_num();
-}
-
-
-int colvarproxy_lammps::smp_num_threads()
-{
- return omp_get_max_threads();
-}
-
-
-int colvarproxy_lammps::smp_lock()
-{
- omp_set_lock(&smp_lock_state);
- return COLVARS_OK;
-}
-
-
-int colvarproxy_lammps::smp_trylock()
-{
- return omp_test_lock(&smp_lock_state) ? COLVARS_OK : COLVARS_ERROR;
-}
-
-
-int colvarproxy_lammps::smp_unlock()
-{
- omp_unset_lock(&smp_lock_state);
- return COLVARS_OK;
-}
-
-#endif
-
-
// multi-replica support
void colvarproxy_lammps::replica_comm_barrier() {
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index 6cdf0edfe8b97272aab39841d8229e8f08d37331..849af410f216894309d98ba158a4b3e2d1012d9f 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -25,10 +25,6 @@
#include <vector>
#include <iostream>
-#if defined(_OPENMP)
-#include <omp.h>
-#endif
-
/* struct for packed data communication of coordinates and forces. */
struct commdata {
int tag,type;
@@ -91,7 +87,8 @@ class colvarproxy_lammps : public colvarproxy {
// methods for lammps to move data or trigger actions in the proxy
public:
void set_temperature(double t) { t_target = t; };
- bool total_forces_enabled() const { return total_force_requested; };
+ bool total_forces_enabled() const { return total_force_requested; };
+ bool total_forces_same_step() const { return true; };
bool want_exit() const { return do_exit; };
// perform colvars computation. returns biasing energy
@@ -140,21 +137,6 @@ class colvarproxy_lammps : public colvarproxy {
// implementation of optional methods from base class
public:
-#if defined(_OPENMP)
- // SMP support
- int smp_enabled();
- int smp_colvars_loop();
- int smp_biases_loop();
- int smp_thread_id();
- int smp_num_threads();
-protected:
- omp_lock_t smp_lock_state;
-public:
- int smp_lock();
- int smp_trylock();
- int smp_unlock();
-#endif
-
// Multi-replica support
// Indicate if multi-replica support is available and active
virtual bool replica_enabled() { return (inter_comm != MPI_COMM_NULL); }
diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h
index 0eb6f2d95ac6bbb4c98fe5c75fc2308a7758342c..e89e00ba039e3c03b32557d28f857216780d6519 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-07-19"
+#define COLVARPROXY_VERSION "2017-10-11"
// 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/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp
index 59e6c46b76113bd1594dcf2acb983ee75212b15b..5897b9c5f0191ae2bb77d3ca1870fa595d987a0c 100644
--- a/src/USER-COLVARS/fix_colvars.cpp
+++ b/src/USER-COLVARS/fix_colvars.cpp
@@ -913,6 +913,7 @@ void FixColvars::write_restart(FILE *fp)
if (me == 0) {
std::string rest_text("");
proxy->serialize_status(rest_text);
+ // TODO call write_output_files()
const char *cvm_state = rest_text.c_str();
int len = strlen(cvm_state) + 1; // need to include terminating NULL byte.
fwrite(&len,sizeof(int),1,fp);