forked from lijiext/lammps
Merge pull request #693 from giacomofiorin/colvars-update
Update Colvars to version 2017-10-11
This commit is contained in:
Steve Plimpton
GitHub
commit
f8f13d929f
Binary file not shown.
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -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++) {
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
||||
};
|
||||
|
|
|
|||
|
|
@ -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() {}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
};
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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\
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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()) {
|
||||
|
|
|
|||
|
|
@ -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() {
|
||||
|
|
|
|||
|
|
@ -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); }
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue