diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 5e90636e1..fcd7d57f0 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1,1934 +1,1943 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvarscript.h"
#include <algorithm>
colvar::colvar(std::string const &conf)
: colvarparse(conf)
{
size_t i, j;
cvm::log("Initializing a new collective variable.\n");
get_keyval(conf, "name", this->name,
(std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));
if (cvm::colvar_by_name(this->name) != NULL) {
cvm::error("Error: this colvar cannot have the same name, \""+this->name+
"\", as another colvar.\n",
INPUT_ERROR);
return;
}
// all tasks disabled by default
for (i = 0; i < task_ntot; i++) {
tasks[i] = false;
}
kinetic_energy = 0.0;
potential_energy = 0.0;
// read the configuration and set up corresponding instances, for
// each type of component implemented
#define initialize_components(def_desc,def_config_key,def_class_name) \
{ \
size_t def_count = 0; \
std::string def_conf = ""; \
size_t pos = 0; \
while ( this->key_lookup(conf, \
def_config_key, \
def_conf, \
pos) ) { \
if (!def_conf.size()) continue; \
cvm::log("Initializing " \
"a new \""+std::string(def_config_key)+"\" component"+ \
(cvm::debug() ? ", with configuration:\n"+def_conf \
: ".\n")); \
cvm::increase_depth(); \
cvc *cvcp = new colvar::def_class_name(def_conf); \
if (cvcp != NULL) { \
cvcs.push_back(cvcp); \
cvcp->check_keywords(def_conf, def_config_key); \
if (cvm::get_error()) { \
cvm::error("Error: in setting up component \"" \
def_config_key"\".\n"); \
return; \
} \
cvm::decrease_depth(); \
} else { \
cvm::error("Error: in allocating component \"" \
def_config_key"\".\n", \
MEMORY_ERROR); \
return; \
} \
if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) { \
if ( (cvcp->function_type != std::string("distance_z")) && \
(cvcp->function_type != std::string("dihedral")) && \
(cvcp->function_type != std::string("spin_angle")) ) { \
cvm::error("Error: invalid use of period and/or " \
"wrapAround in a \""+ \
std::string(def_config_key)+ \
"\" component.\n"+ \
"Period: "+cvm::to_str(cvcp->period) + \
" wrapAround: "+cvm::to_str(cvcp->wrap_center), \
INPUT_ERROR); \
return; \
} \
} \
if ( ! cvcs.back()->name.size()) \
cvcs.back()->name = std::string(def_config_key)+ \
(cvm::to_str(++def_count)); \
if (cvm::debug()) \
cvm::log("Done initializing a \""+ \
std::string(def_config_key)+ \
"\" component"+ \
(cvm::debug() ? \
", named \""+cvcs.back()->name+"\"" \
: "")+".\n"); \
def_conf = ""; \
if (cvm::debug()) \
cvm::log("Parsed "+cvm::to_str(cvcs.size())+ \
" components at this time.\n"); \
} \
}
initialize_components("distance", "distance", distance);
initialize_components("distance vector", "distanceVec", distance_vec);
initialize_components("Cartesian coordinates", "cartesian", cartesian);
initialize_components("distance vector "
"direction", "distanceDir", distance_dir);
initialize_components("distance projection "
"on an axis", "distanceZ", distance_z);
initialize_components("distance projection "
"on a plane", "distanceXY", distance_xy);
initialize_components("average distance weighted by inverse power",
"distanceInv", distance_inv);
initialize_components("N1xN2-long vector of pairwise distances",
"distancePairs", distance_pairs);
initialize_components("coordination "
"number", "coordNum", coordnum);
initialize_components("self-coordination "
"number", "selfCoordNum", selfcoordnum);
initialize_components("angle", "angle", angle);
initialize_components("dihedral", "dihedral", dihedral);
initialize_components("hydrogen bond", "hBond", h_bond);
// initialize_components ("alpha helix", "alphaDihedrals", alpha_dihedrals);
initialize_components("alpha helix", "alpha", alpha_angles);
initialize_components("dihedral principal "
"component", "dihedralPC", dihedPC);
initialize_components("orientation", "orientation", orientation);
initialize_components("orientation "
"angle", "orientationAngle",orientation_angle);
initialize_components("orientation "
"projection", "orientationProj",orientation_proj);
initialize_components("tilt", "tilt", tilt);
initialize_components("spin angle", "spinAngle", spin_angle);
initialize_components("RMSD", "rmsd", rmsd);
// initialize_components ("logarithm of MSD", "logmsd", logmsd);
initialize_components("radius of "
"gyration", "gyration", gyration);
initialize_components("moment of "
"inertia", "inertia", inertia);
initialize_components("moment of inertia around an axis",
"inertiaZ", inertia_z);
initialize_components("eigenvector", "eigenvector", eigenvector);
if (!cvcs.size()) {
cvm::error("Error: no valid components were provided "
"for this collective variable.\n",
INPUT_ERROR);
return;
}
cvm::log("All components initialized.\n");
// Setup colvar as scripted function of components
if (get_keyval(conf, "scriptedFunction", scripted_function,
"", colvarparse::parse_silent)) {
enable(task_scripted);
cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
std::string type_str;
get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
x.type(colvarvalue::type_notset);
int t;
for (t = 0; t < colvarvalue::type_all; t++) {
if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
x.type(colvarvalue::Type(t));
break;
}
}
if (x.type() == colvarvalue::type_notset) {
cvm::error("Could not parse scripted colvar type.");
return;
}
x_reported.type(x.type());
cvm::log(std::string("Expecting colvar value of type ")
+ colvarvalue::type_desc(x.type()));
if (x.type() == colvarvalue::type_vector) {
int size;
if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
cvm::error("Error: no size specified for vector scripted function.");
return;
}
x.vector1d_value.resize(size);
}
// Sort array of cvcs based on values of componentExp
std::vector<cvc *> temp_vec;
for (i = 1; i <= cvcs.size(); i++) {
for (j = 0; j < cvcs.size(); j++) {
if (cvcs[j]->sup_np == int(i)) {
temp_vec.push_back(cvcs[j]);
break;
}
}
}
if (temp_vec.size() != cvcs.size()) {
cvm::error("Could not find order numbers for all components "
"in componentExp values.");
return;
}
cvcs = temp_vec;
// Build ordered list of component values that will be
// passed to the script
for (j = 0; j < cvcs.size(); j++) {
sorted_cvc_values.push_back(&(cvcs[j]->value()));
}
b_homogeneous = false;
// Scripted functions are deemed non-periodic
b_periodic = false;
period = 0.0;
b_inverse_gradients = false;
b_Jacobian_force = false;
}
if (!tasks[task_scripted]) {
colvarvalue const &cvc_value = (cvcs[0])->value();
if (cvm::debug())
cvm::log ("This collective variable is a "+
colvarvalue::type_desc(cvc_value.type())+
((cvc_value.size() > 1) ? " with "+
cvm::to_str(cvc_value.size())+" individual components.\n" :
".\n"));
x.type(cvc_value);
x_reported.type(cvc_value);
}
// If using scripted biases, any colvar may receive bias forces
// and will need its gradient
if (cvm::scripted_forces()) {
enable(task_gradients);
}
// check for linear combinations
b_linear = !tasks[task_scripted];
for (i = 0; i < cvcs.size(); i++) {
if ((cvcs[i])->b_debug_gradients)
enable(task_gradients);
if ((cvcs[i])->sup_np != 1) {
if (cvm::debug() && b_linear)
cvm::log("Warning: You are using a non-linear polynomial "
"combination to define this collective variable, "
"some biasing methods may be unavailable.\n");
b_linear = false;
if ((cvcs[i])->sup_np < 0) {
cvm::log("Warning: you chose a negative exponent in the combination; "
"if you apply forces, the simulation may become unstable "
"when the component \""+
(cvcs[i])->function_type+"\" approaches zero.\n");
}
}
}
// Colvar is homogeneous iff:
// - it is not scripted
// - it is linear
// - all cvcs have coefficient 1 or -1
// i.e. sum or difference of cvcs
b_homogeneous = !tasks[task_scripted] && b_linear;
for (i = 0; i < cvcs.size(); i++) {
if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {
b_homogeneous = false;
}
}
// Colvar is deemed periodic iff:
// - it is homogeneous
// - all cvcs are periodic
// - all cvcs have the same period
b_periodic = cvcs[0]->b_periodic && b_homogeneous;
period = cvcs[0]->period;
for (i = 1; i < cvcs.size(); i++) {
if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
b_periodic = false;
period = 0.0;
}
}
// these will be set to false if any of the cvcs has them false
b_inverse_gradients = !tasks[task_scripted];
b_Jacobian_force = !tasks[task_scripted];
// check the available features of each cvc
for (i = 0; i < cvcs.size(); i++) {
if ((cvcs[i])->b_periodic && !b_periodic) {
cvm::log("Warning: although this component is periodic, the colvar will "
"not be treated as periodic, either because the exponent is not "
"1, or because multiple components are present. Make sure that "
"you know what you are doing!");
}
if (! (cvcs[i])->b_inverse_gradients)
b_inverse_gradients = false;
if (! (cvcs[i])->b_Jacobian_derivative)
b_Jacobian_force = false;
// components may have different types only for scripted functions
if (!tasks[task_scripted] && (colvarvalue::check_types(cvcs[i]->value(),
cvcs[0]->value())) ) {
cvm::error("ERROR: you are definining this collective variable "
"by using components of different types. "
"You must use the same type in order to "
" sum them together.\n", INPUT_ERROR);
return;
}
}
// at this point, the colvar's type is defined
f.type(value());
fb.type(value());
get_keyval(conf, "width", width, 1.0);
if (width <= 0.0) {
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
}
lower_boundary.type(value());
lower_wall.type(value());
upper_boundary.type(value());
upper_wall.type(value());
if (value().type() == colvarvalue::type_scalar) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
enable(task_lower_boundary);
}
get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0);
if (lower_wall_k > 0.0) {
get_keyval(conf, "lowerWall", lower_wall, lower_boundary);
enable(task_lower_wall);
}
if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
enable(task_upper_boundary);
}
get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0);
if (upper_wall_k > 0.0) {
get_keyval(conf, "upperWall", upper_wall, upper_boundary);
enable(task_upper_wall);
}
}
if (tasks[task_lower_boundary]) {
get_keyval(conf, "hardLowerBoundary", hard_lower_boundary, false);
}
if (tasks[task_upper_boundary]) {
get_keyval(conf, "hardUpperBoundary", hard_upper_boundary, false);
}
// consistency checks for boundaries and walls
if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) {
if (lower_boundary >= upper_boundary) {
cvm::error("Error: the upper boundary, "+
cvm::to_str(upper_boundary)+
", is not higher than the lower boundary, "+
cvm::to_str(lower_boundary)+".\n",
INPUT_ERROR);
}
}
if (tasks[task_lower_wall] && tasks[task_upper_wall]) {
if (lower_wall >= upper_wall) {
cvm::error("Error: the upper wall, "+
cvm::to_str(upper_wall)+
", is not higher than the lower wall, "+
cvm::to_str(lower_wall)+".\n",
INPUT_ERROR);
}
if (dist2(lower_wall, upper_wall) < 1.0E-12) {
cvm::log("Lower wall and upper wall are equal "
"in the periodic domain of the colvar: disabling walls.\n");
disable(task_lower_wall);
disable(task_upper_wall);
}
}
get_keyval(conf, "expandBoundaries", expand_boundaries, false);
if (expand_boundaries && periodic_boundaries()) {
cvm::error("Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.\n",
INPUT_ERROR);
}
if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
cvm::error("Error: inconsistent configuration "
"(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).\n",
INPUT_ERROR);
}
{
bool b_extended_lagrangian;
get_keyval(conf, "extendedLagrangian", b_extended_lagrangian, false);
if (b_extended_lagrangian) {
cvm::real temp, tolerance, period;
cvm::log("Enabling the extended Lagrangian term for colvar \""+
this->name+"\".\n");
enable(task_extended_lagrangian);
xr.type(value());
vr.type(value());
fr.type(value());
const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
if (temp <= 0.0) {
if (found)
cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
else
cvm::error("Error: a positive temperature must be provided, either "
"by enabling a thermostat, or through \"extendedTemp\".\n",
INPUT_ERROR);
}
get_keyval(conf, "extendedFluctuation", tolerance);
if (tolerance <= 0.0) {
cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
}
ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
get_keyval(conf, "extendedTimeConstant", period, 200.0);
if (period <= 0.0) {
cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
}
ext_mass = (cvm::boltzmann() * temp * period * period)
/ (4.0 * PI * PI * tolerance * tolerance);
cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");
{
bool b_output_energy;
get_keyval(conf, "outputEnergy", b_output_energy, false);
if (b_output_energy) {
enable(task_output_energy);
}
}
get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
if (ext_gamma < 0.0) {
cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
}
if (ext_gamma != 0.0) {
enable(task_langevin);
ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1
ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());
}
}
}
{
bool b_output_value;
get_keyval(conf, "outputValue", b_output_value, true);
if (b_output_value) {
enable(task_output_value);
}
}
{
bool b_output_velocity;
get_keyval(conf, "outputVelocity", b_output_velocity, false);
if (b_output_velocity) {
enable(task_output_velocity);
}
}
{
bool b_output_system_force;
get_keyval(conf, "outputSystemForce", b_output_system_force, false);
if (b_output_system_force) {
enable(task_output_system_force);
}
}
{
bool b_output_applied_force;
get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);
if (b_output_applied_force) {
enable(task_output_applied_force);
}
}
if (cvm::b_analysis)
parse_analysis(conf);
if (cvm::debug())
cvm::log("Done initializing collective variable \""+this->name+"\".\n");
}
void colvar::build_atom_list(void)
{
// If atomic gradients are requested, build full list of atom ids from all cvcs
std::list<int> temp_id_list;
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
temp_id_list.push_back(cvcs[i]->atom_groups[j]->at(k).id);
}
}
}
temp_id_list.sort();
temp_id_list.unique();
// atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
unsigned int id_i = 0;
std::list<int>::iterator li;
for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
atom_ids[id_i] = *li;
id_i++;
}
temp_id_list.clear();
atomic_gradients.resize(atom_ids.size());
if (atom_ids.size()) {
if (cvm::debug())
cvm::log("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
} else {
cvm::log("Warning: colvar components communicated no atom IDs.\n");
}
}
int colvar::parse_analysis(std::string const &conf)
{
// if (cvm::debug())
// cvm::log ("Parsing analysis flags for collective variable \""+
// this->name+"\".\n");
runave_length = 0;
bool b_runave = false;
if (get_keyval(conf, "runAve", b_runave) && b_runave) {
enable(task_runave);
get_keyval(conf, "runAveLength", runave_length, 1000);
get_keyval(conf, "runAveStride", runave_stride, 1);
if ((cvm::restart_out_freq % runave_stride) != 0) {
cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
}
std::string runave_outfile;
get_keyval(conf, "runAveOutputFile", runave_outfile,
std::string(cvm::output_prefix+"."+
this->name+".runave.traj"));
size_t const this_cv_width = x.output_width(cvm::cv_width);
runave_os.open(runave_outfile.c_str());
runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
<< " "
<< cvm::wrap_string("running average", this_cv_width)
<< " "
<< cvm::wrap_string("running stddev", this_cv_width)
<< "\n";
}
acf_length = 0;
bool b_acf = false;
if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {
enable(task_corrfunc);
std::string acf_colvar_name;
get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
if (acf_colvar_name == this->name) {
cvm::log("Calculating auto-correlation function.\n");
} else {
cvm::log("Calculating correlation function with \""+
this->name+"\".\n");
}
std::string acf_type_str;
get_keyval(conf, "corrFuncType", acf_type_str, to_lower_cppstr(std::string("velocity")));
if (acf_type_str == to_lower_cppstr(std::string("coordinate"))) {
acf_type = acf_coor;
} else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
acf_type = acf_vel;
enable(task_fdiff_velocity);
if (acf_colvar_name.size())
(cvm::colvar_by_name(acf_colvar_name))->enable(task_fdiff_velocity);
} else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
acf_type = acf_p2coor;
} else {
cvm::log("Unknown type of correlation function, \""+
acf_type_str+"\".\n");
cvm::set_error_bits(INPUT_ERROR);
}
get_keyval(conf, "corrFuncOffset", acf_offset, 0);
get_keyval(conf, "corrFuncLength", acf_length, 1000);
get_keyval(conf, "corrFuncStride", acf_stride, 1);
if ((cvm::restart_out_freq % acf_stride) != 0) {
cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
}
get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
get_keyval(conf, "corrFuncOutputFile", acf_outfile,
std::string(cvm::output_prefix+"."+this->name+
".corrfunc.dat"));
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvar::enable(colvar::task const &t)
{
switch (t) {
case task_output_system_force:
enable(task_system_force);
break;
case task_report_Jacobian_force:
enable(task_Jacobian_force);
enable(task_system_force);
if (cvm::debug())
cvm::log("Adding the Jacobian force to the system force, "
"rather than applying its opposite silently.\n");
break;
case task_Jacobian_force:
// checks below do not apply to extended-system colvars
if ( !tasks[task_extended_lagrangian] ) {
enable(task_gradients);
if (!b_Jacobian_force) {
cvm::error("Error: colvar \""+this->name+
"\" does not have Jacobian forces implemented.\n",
INPUT_ERROR);
}
if (!b_linear) {
cvm::error("Error: colvar \""+this->name+
"\" must be defined as a linear combination "
"to calculate the Jacobian force.\n",
INPUT_ERROR);
}
if (cvm::debug())
cvm::log("Enabling calculation of the Jacobian force "
"on this colvar.\n");
}
fj.type(value());
break;
case task_system_force:
if (!tasks[task_extended_lagrangian]) {
if (!b_inverse_gradients) {
cvm::error("Error: one or more of the components of "
"colvar \""+this->name+
"\" does not support system force calculation.\n",
INPUT_ERROR);
}
cvm::request_system_force();
}
ft.type(value());
ft_reported.type(value());
break;
case task_output_applied_force:
case task_lower_wall:
case task_upper_wall:
// all of the above require gradients
enable(task_gradients);
break;
case task_fdiff_velocity:
x_old.type(value());
v_fdiff.type(value());
v_reported.type(value());
break;
case task_output_velocity:
enable(task_fdiff_velocity);
break;
case task_grid:
if (value().type() != colvarvalue::type_scalar) {
cvm::error("Cannot calculate a grid for collective variable, \""+
this->name+"\", because its value is not a scalar number.\n",
INPUT_ERROR);
}
break;
case task_extended_lagrangian:
enable(task_gradients);
v_reported.type(value());
break;
case task_lower_boundary:
case task_upper_boundary:
if (value().type() != colvarvalue::type_scalar) {
cvm::error("Error: this colvar is not a scalar value "
"and cannot produce a grid.\n",
INPUT_ERROR);
}
break;
case task_output_value:
case task_runave:
case task_corrfunc:
case task_ntot:
case task_langevin:
case task_output_energy:
case task_scripted:
case task_gradients:
break;
case task_collect_gradients:
if (value().type() != colvarvalue::type_scalar) {
cvm::error("Collecting atomic gradients for non-scalar collective variable \""+
this->name+"\" is not yet implemented.\n",
INPUT_ERROR);
}
enable(task_gradients);
if (atom_ids.size() == 0) {
build_atom_list();
}
break;
}
tasks[t] = true;
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void colvar::disable(colvar::task const &t)
{
// check dependencies
switch (t) {
case task_gradients:
disable(task_upper_wall);
disable(task_lower_wall);
disable(task_output_applied_force);
disable(task_system_force);
disable(task_Jacobian_force);
break;
case task_system_force:
disable(task_output_system_force);
break;
case task_Jacobian_force:
disable(task_report_Jacobian_force);
break;
case task_fdiff_velocity:
disable(task_output_velocity);
break;
case task_lower_boundary:
case task_upper_boundary:
disable(task_grid);
break;
case task_extended_lagrangian:
case task_report_Jacobian_force:
case task_output_value:
case task_output_velocity:
case task_output_applied_force:
case task_output_system_force:
case task_runave:
case task_corrfunc:
case task_grid:
case task_lower_wall:
case task_upper_wall:
case task_ntot:
case task_langevin:
case task_output_energy:
case task_collect_gradients:
case task_scripted:
break;
}
tasks[t] = false;
}
void colvar::setup() {
// loop over all components to reset masses of all groups
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
atoms.read_positions();
atoms.reset_mass(name,i,ig);
}
}
}
colvar::~colvar()
{
for (size_t i = 0; i < cvcs.size(); i++) {
delete cvcs[i];
}
// remove reference to this colvar from the CVM
for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
cvi != cvm::colvars.end();
++cvi) {
if ( *cvi == this) {
cvm::colvars.erase(cvi);
break;
}
}
}
// ******************** CALC FUNCTIONS ********************
void colvar::calc()
{
size_t i, ig;
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\".\n");
// prepare atom groups for calculation
if (cvm::debug())
cvm::log("Collecting data from atom groups.\n");
+ // Update the enabled/disabled status of cvcs if necessary
+ if (cvc_flags.size()) {
+ bool any = false;
+ for (i = 0; i < cvcs.size(); i++) {
+ cvcs[i]->b_enabled = cvc_flags[i];
+ any = any || cvc_flags[i];
+ }
+ if (!any) {
+ cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
+ return;
+ }
+ cvc_flags.resize(0);
+ }
+
for (i = 0; i < cvcs.size(); i++) {
+ if (!cvcs[i]->b_enabled) continue;
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
atoms.reset_atoms_data();
atoms.read_positions();
if (atoms.b_center || atoms.b_rotate) {
atoms.calc_apply_roto_translation();
}
// each atom group will take care of its own ref_pos_group, if defined
}
}
//// Don't try to get atom velocities, as no back-end currently implements it
// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
// for (i = 0; i < cvcs.size(); i++) {
// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
// cvcs[i]->atom_groups[ig]->read_velocities();
// }
// }
// }
if (tasks[task_system_force]) {
for (i = 0; i < cvcs.size(); i++) {
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvcs[i]->atom_groups[ig]->read_system_forces();
}
}
}
// calculate the value of the colvar
if (cvm::debug())
cvm::log("Calculating colvar components.\n");
x.reset();
// First, update component values
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
(cvcs[i])->calc_value();
cvm::decrease_depth();
if (cvm::debug())
cvm::log("Colvar component no. "+cvm::to_str(i+1)+
" within colvar \""+this->name+"\" has value "+
cvm::to_str((cvcs[i])->value(),
cvm::cv_width, cvm::cv_prec)+".\n");
}
// Then combine them appropriately
if (tasks[task_scripted]) {
// cvcs combined by user script
int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Scripted colvars are not implemented.");
return;
}
if (res != COLVARS_OK) {
cvm::error("Error running scripted colvar");
return;
}
} else if (x.type() == colvarvalue::type_scalar) {
// polynomial combination allowed
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
x += (cvcs[i])->sup_coeff *
( ((cvcs[i])->sup_np != 1) ?
std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
(cvcs[i])->value().real_value );
}
} else {
// only linear combination allowed
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
}
}
if (cvm::debug())
cvm::log("Colvar \""+this->name+"\" has value "+
cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
if (tasks[task_gradients]) {
if (cvm::debug())
cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
for (i = 0; i < cvcs.size(); i++) {
// calculate the gradients of each component
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
(cvcs[i])->calc_gradients();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
cvcs[i]->atom_groups[ig]->calc_fit_gradients();
}
cvm::decrease_depth();
}
if (cvm::debug())
cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
if (tasks[task_collect_gradients]) {
if (tasks[task_scripted]) {
cvm::error("Collecting atomic gradients is not implemented for "
"scripted colvars.");
return;
}
// Collect the atomic gradients inside colvar object
for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
atomic_gradients[a].reset();
}
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
// Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
// If necessary, apply inverse rotation to get atomic
// gradient in the laboratory frame
if (cvcs[i]->atom_groups[j]->b_rotate) {
cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse();
for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
atomic_gradients[a] += coeff *
rot_inv.rotate(cvcs[i]->atom_groups[j]->at(k).grad);
}
} else {
for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad;
}
}
}
}
}
}
if (tasks[task_system_force] && !tasks[task_extended_lagrangian]) {
// If extended Lagrangian is enabled, system force calculation is trivial
// and done together with integration of the extended coordinate.
if (tasks[task_scripted]) {
// TODO see if this could reasonably be done in a generic way
// from generic inverse gradients
cvm::error("System force is not implemented for "
"scripted colvars.");
return;
}
if (cvm::debug())
cvm::log("Calculating system force of colvar \""+this->name+"\".\n");
ft.reset();
// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
// Disabled check to allow for explicit system force calculation
// even with extended Lagrangian
if (cvm::step_relative() > 0) {
// get from the cvcs the system forces from the PREVIOUS step
for (i = 0; i < cvcs.size(); i++) {
(cvcs[i])->calc_force_invgrads();
// linear combination is assumed
cvm::increase_depth();
ft += (cvcs[i])->system_force() / ((cvcs[i])->sup_coeff * cvm::real(cvcs.size()));
cvm::decrease_depth();
}
}
if (tasks[task_report_Jacobian_force]) {
// add the Jacobian force to the system force, and don't apply any silent
// correction internally: biases such as colvarbias_abf will handle it
ft += fj;
}
if (cvm::debug())
cvm::log("Done calculating system force of colvar \""+this->name+"\".\n");
}
if (tasks[task_fdiff_velocity]) {
// calculate the velocity by finite differences
if (cvm::step_relative() == 0)
x_old = x;
else {
v_fdiff = fdiff_velocity(x_old, x);
v_reported = v_fdiff;
}
}
if (tasks[task_extended_lagrangian]) {
// initialize the restraint center in the first step to the value
// just calculated from the cvcs
// TODO: put it in the restart information
if (cvm::step_relative() == 0) {
xr = x;
vr = 0.0; // (already 0; added for clarity)
}
// report the restraint center as "value"
x_reported = xr;
v_reported = vr;
// the "system force" with the extended Lagrangian is just the
// harmonic term acting on the extended coordinate
// Note: this is the force for current timestep
ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
} else {
x_reported = x;
ft_reported = ft;
}
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\".\n");
}
cvm::real colvar::update()
{
if (cvm::debug())
cvm::log("Updating colvar \""+this->name+"\".\n");
// set to zero the applied force
f.type(value());
f.reset();
// add the biases' force, which at this point should already have
// been summed over each bias using this colvar
f += fb;
if (tasks[task_lower_wall] || tasks[task_upper_wall]) {
// wall force
colvarvalue fw(value());
fw.reset();
if (cvm::debug())
cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
// if the two walls are applied concurrently, decide which is the
// closer one (on a periodic colvar, both walls may be applicable
// at the same time)
if ( (!tasks[task_upper_wall]) ||
(this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
cvm::real const grad = this->dist2_lgrad(x, lower_wall);
if (grad < 0.0) {
fw = -0.5 * lower_wall_k * grad;
if (cvm::debug())
cvm::log("Applying a lower wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
f += fw;
}
} else {
cvm::real const grad = this->dist2_lgrad(x, upper_wall);
if (grad > 0.0) {
fw = -0.5 * upper_wall_k * grad;
if (cvm::debug())
cvm::log("Applying an upper wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
f += fw;
}
}
}
if (tasks[task_Jacobian_force]) {
size_t i;
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
(cvcs[i])->calc_Jacobian_derivative();
cvm::decrease_depth();
}
size_t ncvc = 0;
fj.reset();
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
// linear combination is assumed
fj += 1.0 / cvm::real((cvcs[i])->sup_coeff) *
(cvcs[i])->Jacobian_derivative();
ncvc++;
}
fj *= (1.0/cvm::real(ncvc)) * cvm::boltzmann() * cvm::temperature();
// the instantaneous Jacobian force was not included in the reported system force;
// instead, it is subtracted from the applied force (silent Jacobian correction)
if (! tasks[task_report_Jacobian_force])
f -= fj;
}
if (tasks[task_extended_lagrangian]) {
cvm::real dt = cvm::dt();
cvm::real f_ext;
// the total force is applied to the fictitious mass, while the
// atoms only feel the harmonic force
// fr: extended coordinate force (without harmonic spring), for output in trajectory
// f_ext: total force on extended coordinate (including harmonic spring)
// f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
// Because of leapfrog, kinetic energy at time i is approximate
kinetic_energy = 0.5 * ext_mass * vr * vr;
potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
// leap to v_(i+1/2)
if (tasks[task_langevin]) {
vr -= dt * ext_gamma * vr.real_value;
vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
}
vr += (0.5 * dt) * f_ext / ext_mass;
xr += dt * vr;
xr.apply_constraints();
if (this->b_periodic) this->wrap(xr);
}
if (tasks[task_fdiff_velocity]) {
// set it for the next step
x_old = x;
}
if (cvm::debug())
cvm::log("Done updating colvar \""+this->name+"\".\n");
return (potential_energy + kinetic_energy);
}
void colvar::communicate_forces()
{
size_t i;
if (cvm::debug())
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
if (tasks[task_scripted]) {
std::vector<cvm::matrix2d<cvm::real> > func_grads;
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
cvcs[i]->value().size()));
}
int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
if (res != COLVARS_OK) {
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Colvar gradient scripts are not implemented.");
} else {
cvm::error("Error running colvar gradient script");
}
return;
}
+ int grad_index = 0; // index in the scripted gradients, to account for some components being disabled
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
// cvc force is colvar force times colvar/cvc Jacobian
// (vector-matrix product)
- (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[i],
+ (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
cvcs[i]->value().type()));
cvm::decrease_depth();
}
} else if (x.type() == colvarvalue::type_scalar) {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
cvm::real((cvcs[i])->sup_np) *
(std::pow((cvcs[i])->value().real_value,
(cvcs[i])->sup_np-1)) );
cvm::decrease_depth();
}
} else {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->b_enabled) continue;
cvm::increase_depth();
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
cvm::decrease_depth();
}
}
if (cvm::debug())
cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
}
int colvar::set_cvc_flags(std::vector<bool> const &flags) {
- size_t i;
if (flags.size() != cvcs.size()) {
cvm::error("ERROR: Wrong number of CVC flags provided.");
return COLVARS_ERROR;
}
- bool e = false;
- for (i = 0; i < cvcs.size(); i++) {
- cvcs[i]->b_enabled = flags[i];
- e = e || flags[i];
- }
- if (!e) {
- cvm::error("ERROR: All CVCs are disabled for this colvar.");
- return COLVARS_ERROR;
- }
+ // We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
+ // so we store the flags that will be enforced at the next call to calc()
+ cvc_flags = flags;
return COLVARS_OK;
}
// ******************** METRIC FUNCTIONS ********************
// Use the metrics defined by \link cvc \endlink objects
bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
{
if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) {
cvm::log("Error: requesting to histogram the "
"collective variable \""+this->name+"\", but a "
"pair of lower and upper boundaries must be "
"defined.\n");
cvm::set_error_bits(INPUT_ERROR);
}
if (period > 0.0) {
if ( ((std::sqrt(this->dist2(lb, ub))) / this->width)
< 1.0E-10 ) {
return true;
}
}
return false;
}
bool colvar::periodic_boundaries() const
{
if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) {
cvm::error("Error: requesting to histogram the "
"collective variable \""+this->name+"\", but a "
"pair of lower and upper boundaries must be "
"defined.\n");
}
return periodic_boundaries(lower_boundary, upper_boundary);
}
cvm::real colvar::dist2(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (b_homogeneous) {
return (cvcs[0])->dist2(x1, x2);
} else {
return x1.dist2(x2);
}
}
colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (b_homogeneous) {
return (cvcs[0])->dist2_lgrad(x1, x2);
} else {
return x1.dist2_grad(x2);
}
}
colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (b_homogeneous) {
return (cvcs[0])->dist2_rgrad(x1, x2);
} else {
return x2.dist2_grad(x1);
}
}
void colvar::wrap(colvarvalue &x) const
{
if (b_homogeneous) {
(cvcs[0])->wrap(x);
}
return;
}
// ******************** INPUT FUNCTIONS ********************
std::istream & colvar::read_restart(std::istream &is)
{
size_t const start_pos = is.tellg();
std::string conf;
if ( !(is >> colvarparse::read_block("colvar", conf)) ) {
// this is not a colvar block
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
{
std::string check_name = "";
if ( (get_keyval(conf, "name", check_name,
std::string(""), colvarparse::parse_silent)) &&
(check_name != name) ) {
cvm::error("Error: the state file does not match the "
"configuration file, at colvar \""+name+"\".\n");
}
if (check_name.size() == 0) {
cvm::error("Error: Collective variable in the "
"restart file without any identifier.\n");
}
}
if ( !(get_keyval(conf, "x", x,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"the value of the colvar \""+
name+"\" .\n");
} else {
cvm::log("Restarting collective variable \""+name+"\" from value: "+
cvm::to_str(x)+"\n");
}
if (tasks[colvar::task_extended_lagrangian]) {
if ( !(get_keyval(conf, "extended_x", xr,
colvarvalue(x.type()), colvarparse::parse_silent)) &&
!(get_keyval(conf, "extended_v", vr,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"\"extended_x\" or \"extended_v\" for the colvar \""+
name+"\", but you requested \"extendedLagrangian\".\n");
}
}
if (tasks[task_extended_lagrangian]) {
x_reported = xr;
} else {
x_reported = x;
}
if (tasks[task_output_velocity]) {
if ( !(get_keyval(conf, "v", v_fdiff,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"the velocity for the colvar \""+
name+"\", but you requested \"outputVelocity\".\n");
}
if (tasks[task_extended_lagrangian]) {
v_reported = vr;
} else {
v_reported = v_fdiff;
}
}
return is;
}
std::istream & colvar::read_traj(std::istream &is)
{
size_t const start_pos = is.tellg();
if (tasks[task_output_value]) {
if (!(is >> x)) {
cvm::log("Error: in reading the value of colvar \""+
this->name+"\" from trajectory.\n");
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
if (tasks[task_extended_lagrangian]) {
is >> xr;
x_reported = xr;
} else {
x_reported = x;
}
}
if (tasks[task_output_velocity]) {
is >> v_fdiff;
if (tasks[task_extended_lagrangian]) {
is >> vr;
v_reported = vr;
} else {
v_reported = v_fdiff;
}
}
if (tasks[task_output_system_force]) {
is >> ft;
ft_reported = ft;
}
if (tasks[task_output_applied_force]) {
is >> f;
}
return is;
}
// ******************** OUTPUT FUNCTIONS ********************
std::ostream & colvar::write_restart(std::ostream &os) {
os << "colvar {\n"
<< " name " << name << "\n"
<< " x "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< x << "\n";
if (tasks[task_output_velocity]) {
os << " v "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< v_reported << "\n";
}
if (tasks[task_extended_lagrangian]) {
os << " extended_x "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< xr << "\n"
<< " extended_v "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< vr << "\n";
}
os << "}\n\n";
return os;
}
std::ostream & colvar::write_traj_label(std::ostream & os)
{
size_t const this_cv_width = x.output_width(cvm::cv_width);
os << " ";
if (tasks[task_output_value]) {
os << " "
<< cvm::wrap_string(this->name, this_cv_width);
if (tasks[task_extended_lagrangian]) {
// extended DOF
os << " r_"
<< cvm::wrap_string(this->name, this_cv_width-2);
}
}
if (tasks[task_output_velocity]) {
os << " v_"
<< cvm::wrap_string(this->name, this_cv_width-2);
if (tasks[task_extended_lagrangian]) {
// extended DOF
os << " vr_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
}
if (tasks[task_output_energy]) {
os << " Ep_"
<< cvm::wrap_string(this->name, this_cv_width-3)
<< " Ek_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
if (tasks[task_output_system_force]) {
os << " fs_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
if (tasks[task_output_applied_force]) {
os << " fa_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
return os;
}
std::ostream & colvar::write_traj(std::ostream &os)
{
os << " ";
if (tasks[task_output_value]) {
if (tasks[task_extended_lagrangian]) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< x;
}
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< x_reported;
}
if (tasks[task_output_velocity]) {
if (tasks[task_extended_lagrangian]) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< v_fdiff;
}
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< v_reported;
}
if (tasks[task_output_energy]) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< potential_energy
<< " "
<< kinetic_energy;
}
if (tasks[task_output_system_force]) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< ft_reported;
}
if (tasks[task_output_applied_force]) {
if (tasks[task_extended_lagrangian]) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< fr;
} else {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< f;
}
}
return os;
}
int colvar::write_output_files()
{
if (cvm::b_analysis) {
if (acf.size()) {
cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
cvm::ofstream acf_os(acf_outfile.c_str());
if (! acf_os.good()) {
cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR);
}
write_acf(acf_os);
acf_os.close();
}
if (runave_os.good()) {
runave_os.close();
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
// ******************** ANALYSIS FUNCTIONS ********************
void colvar::analyze()
{
if (tasks[task_runave]) {
calc_runave();
}
if (tasks[task_corrfunc]) {
calc_acf();
}
}
inline void history_add_value(size_t const &history_length,
std::list<colvarvalue> &history,
colvarvalue const &new_value)
{
history.push_front(new_value);
if (history.size() > history_length)
history.pop_back();
}
inline void history_incr(std::list< std::list<colvarvalue> > &history,
std::list< std::list<colvarvalue> >::iterator &history_p)
{
if ((++history_p) == history.end())
history_p = history.begin();
}
int colvar::calc_acf()
{
// using here an acf_stride-long list of vectors for either
// coordinates(acf_x_history) or velocities (acf_v_history); each vector can
// contain up to acf_length values, which are contiguous in memory
// representation but separated by acf_stride in the time series;
// the pointer to each vector is changed at every step
if (acf_x_history.empty() && acf_v_history.empty()) {
// first-step operations
colvar *cfcv = (acf_colvar_name.size() ?
cvm::colvar_by_name(acf_colvar_name) :
this);
if (colvarvalue::check_types(cfcv->value(), value())) {
cvm::error("Error: correlation function between \""+cfcv->name+
"\" and \""+this->name+"\" cannot be calculated, "
"because their value types are different.\n",
INPUT_ERROR);
}
acf_nframes = 0;
cvm::log("Colvar \""+this->name+"\": initializing ACF calculation.\n");
if (acf.size() < acf_length+1)
acf.resize(acf_length+1, 0.0);
size_t i;
switch (acf_type) {
case acf_vel:
// allocate space for the velocities history
for (i = 0; i < acf_stride; i++) {
acf_v_history.push_back(std::list<colvarvalue>());
}
acf_v_history_p = acf_v_history.begin();
break;
case acf_coor:
case acf_p2coor:
// allocate space for the coordinates history
for (i = 0; i < acf_stride; i++) {
acf_x_history.push_back(std::list<colvarvalue>());
}
acf_x_history_p = acf_x_history.begin();
break;
default:
break;
}
} else {
colvar *cfcv = (acf_colvar_name.size() ?
cvm::colvar_by_name(acf_colvar_name) :
this);
switch (acf_type) {
case acf_vel:
if (tasks[task_fdiff_velocity]) {
// calc() should do this already, but this only happens in a
// simulation; better do it again in case a trajectory is
// being read
v_reported = v_fdiff = fdiff_velocity(x_old, cfcv->value());
}
calc_vel_acf((*acf_v_history_p), cfcv->velocity());
// store this value in the history
history_add_value(acf_length+acf_offset, *acf_v_history_p, cfcv->velocity());
// if stride is larger than one, cycle among different histories
history_incr(acf_v_history, acf_v_history_p);
break;
case acf_coor:
calc_coor_acf((*acf_x_history_p), cfcv->value());
history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
history_incr(acf_x_history, acf_x_history_p);
break;
case acf_p2coor:
calc_p2coor_acf((*acf_x_history_p), cfcv->value());
history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
history_incr(acf_x_history, acf_x_history_p);
break;
default:
break;
}
}
if (tasks[task_fdiff_velocity]) {
// set it for the next step
x_old = x;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
colvarvalue const &v)
{
// loop over stored velocities and add to the ACF, but only the
// length is sufficient to hold an entire row of ACF values
if (v_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator vs_i = v_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++vs_i;
// current vel with itself
*(acf_i) += v.norm2();
++acf_i;
// inner products of previous velocities with current (acf_i and
// vs_i are updated)
colvarvalue::inner_opt(v, vs_i, v_list.end(), acf_i);
acf_nframes++;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
colvarvalue const &x)
{
// same as above but for coordinates
if (x_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator xs_i = x_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++xs_i;
*(acf_i++) += x.norm2();
colvarvalue::inner_opt(x, xs_i, x_list.end(), acf_i);
acf_nframes++;
}
}
void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
colvarvalue const &x)
{
// same as above but with second order Legendre polynomial instead
// of just the scalar product
if (x_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator xs_i = x_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++xs_i;
// value of P2(0) = 1
*(acf_i++) += 1.0;
colvarvalue::p2leg_opt(x, xs_i, x_list.end(), acf_i);
acf_nframes++;
}
}
void colvar::write_acf(std::ostream &os)
{
if (!acf_nframes)
cvm::log("Warning: ACF was not calculated (insufficient frames).\n");
os.setf(std::ios::scientific, std::ios::floatfield);
os << "# Autocorrelation function for collective variable \""
<< this->name << "\"\n";
// one frame is used for normalization, the statistical sample is
// hence decreased
os << "# nframes = " << (acf_normalize ?
acf_nframes - 1 :
acf_nframes) << "\n";
cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
std::vector<cvm::real>::iterator acf_i;
size_t it = acf_offset;
for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
os << std::setw(cvm::it_width) << acf_stride * (it++) << " "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< ( acf_normalize ?
(*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
(*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
}
}
void colvar::calc_runave()
{
if (x_history.empty()) {
runave.type(value().type());
runave.reset();
// first-step operations
if (cvm::debug())
cvm::log("Colvar \""+this->name+
"\": initializing running average calculation.\n");
acf_nframes = 0;
x_history.push_back(std::list<colvarvalue>());
x_history_p = x_history.begin();
} else {
if ( (cvm::step_relative() % runave_stride) == 0) {
if ((*x_history_p).size() >= runave_length-1) {
runave = x;
std::list<colvarvalue>::iterator xs_i;
for (xs_i = (*x_history_p).begin();
xs_i != (*x_history_p).end(); ++xs_i) {
runave += (*xs_i);
}
runave *= 1.0 / cvm::real(runave_length);
runave.apply_constraints();
runave_variance = 0.0;
runave_variance += this->dist2(x, runave);
for (xs_i = (*x_history_p).begin();
xs_i != (*x_history_p).end(); ++xs_i) {
runave_variance += this->dist2(x, (*xs_i));
}
runave_variance *= 1.0 / cvm::real(runave_length-1);
runave_os << std::setw(cvm::it_width) << cvm::step_relative()
<< " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< runave << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< std::sqrt(runave_variance) << "\n";
}
history_add_value(runave_length, *x_history_p, x);
}
}
}
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 79b00c026..8a3e40518 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -1,580 +1,583 @@
/// -*- c++ -*-
#ifndef COLVAR_H
#define COLVAR_H
#include <iostream>
#include <iomanip>
#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
/// \brief A collective variable (main class); to be defined, it needs
/// at least one object of a derived class of colvar::cvc; it
/// calculates and returns a \link colvarvalue \endlink object
///
/// This class parses the configuration, defines the behaviour and
/// stores the value (\link colvar::x \endlink) and all related data
/// of a collective variable. How the value is calculated is defined
/// in \link colvar::cvc \endlink and its derived classes. The
/// \link colvar \endlink object contains pointers to multiple \link
/// colvar::cvc \endlink derived objects, which can be combined
/// together into one collective variable. This makes possible to
/// implement new collective variables at runtime based on the
/// existing ones. Currently, this possibility is limited to a
/// polynomial, using the coefficients cvc::sup_coeff and the
/// exponents cvc::sup_np. In case of non-scalar variables,
/// only exponents equal to 1 are accepted.
///
/// Please note that most of its members are \link colvarvalue
/// \endlink objects, i.e. they can handle different data types
/// together, and must all be set to the same type of colvar::value()
/// before using them together in assignments or other operations; this is usually done
/// automatically in the constructor. If you add a new member of
/// \link colvarvalue \endlink type, you should also add its
/// initialization line in the \link colvar \endlink constructor.
class colvar : public colvarparse {
public:
/// Name
std::string name;
/// \brief Current value (previously set by calc() or by read_traj())
colvarvalue const & value() const;
/// \brief Current actual value (not extended DOF)
colvarvalue const & actual_value() const;
/// \brief Current velocity (previously set by calc() or by read_traj())
colvarvalue const & velocity() const;
/// \brief Current system force (previously obtained from calc() or
/// read_traj()). Note: this is calculated using the atomic forces
/// from the last simulation step.
///
/// Total atom forces are read from the MD program, the total force
/// acting on the collective variable is calculated summing those
/// from all colvar components, the bias and walls forces are
/// subtracted.
colvarvalue const & system_force() const;
/// \brief
/// \brief Typical fluctuation amplitude for this collective
/// variable (e.g. local width of a free energy basin)
///
/// In metadynamics calculations, \link colvarbias_meta \endlink,
/// this value is used to calculate the width of a gaussian. In ABF
/// calculations, \link colvarbias_abf \endlink, it is used to
/// calculate the grid spacing in the direction of this collective
/// variable.
cvm::real width;
/// \brief True if this \link colvar \endlink is a linear
/// combination of \link cvc \endlink elements
bool b_linear;
/// \brief True if this \link colvar \endlink is a linear
/// combination of cvcs with coefficients 1 or -1
bool b_homogeneous;
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating inverse gradients
bool b_inverse_gradients;
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating Jacobian forces
bool b_Jacobian_force;
/// \brief Options controlling the behaviour of the colvar during
/// the simulation, which are set from outside the cvcs
enum task {
/// \brief Gradients are calculated and temporarily stored, so
/// that external forces can be applied
task_gradients,
/// \brief Collect atomic gradient data from all cvcs into vector
/// atomic_gradients
task_collect_gradients,
/// \brief Calculate the velocity with finite differences
task_fdiff_velocity,
/// \brief The system force is calculated, projecting the atomic
/// forces on the inverse gradients
task_system_force,
/// \brief The variable has a harmonic restraint around a moving
/// center with fictitious mass; bias forces will be applied to
/// the center
task_extended_lagrangian,
/// \brief The extended system coordinate undergoes Langevin
/// dynamics
task_langevin,
/// \brief Output the potential and kinetic energies
/// (for extended Lagrangian colvars only)
task_output_energy,
/// \brief Compute analytically the "force" arising from the
/// geometric entropy component (for example, that from the angular
/// states orthogonal to a distance vector)
task_Jacobian_force,
/// \brief Report the Jacobian force as part of the system force
/// if disabled, apply a correction internally to cancel it
task_report_Jacobian_force,
/// \brief Output the value to the trajectory file (on by default)
task_output_value,
/// \brief Output the velocity to the trajectory file
task_output_velocity,
/// \brief Output the applied force to the trajectory file
task_output_applied_force,
/// \brief Output the system force to the trajectory file
task_output_system_force,
/// \brief A lower boundary is defined
task_lower_boundary,
/// \brief An upper boundary is defined
task_upper_boundary,
/// \brief Provide a discretization of the values of the colvar to
/// be used by the biases or in analysis (needs lower and upper
/// boundary)
task_grid,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes below the lower wall
task_lower_wall,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes above the upper wall
task_upper_wall,
/// \brief Compute running average
task_runave,
/// \brief Compute time correlation function
task_corrfunc,
/// \brief Value and gradients computed by user script
task_scripted,
/// \brief Number of possible tasks
task_ntot
};
/// Tasks performed by this colvar
bool tasks[task_ntot];
/// List of biases that depend on this colvar
std::vector<colvarbias *> biases;
protected:
/*
extended:
H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
+ \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
+ \sum_{t'<t} W * exp(-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
+ \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M
normal:
H = H_{0} + \sum_{t'<t} W * exp(-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
+ \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M
output:
H = H_{0} (only output S(x), no forces)
Here:
S(x(t)) = x
s(t) = xr
DS = Ds = delta
*/
/// Value of the colvar
colvarvalue x;
// TODO: implement functionality to treat these
// /// Vector of individual values from CVCs
// colvarvalue x_cvc;
// /// Jacobian matrix of individual values from CVCs
// colvarvalue dx_cvc;
/// Cached reported value (x may be manipulated)
colvarvalue x_reported;
/// Finite-difference velocity
colvarvalue v_fdiff;
inline colvarvalue fdiff_velocity(colvarvalue const &xold,
colvarvalue const &xnew)
{
// using the gradient of the square distance to calculate the
// velocity (non-scalar variables automatically taken into
// account)
cvm::real dt = cvm::dt();
return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) *
0.5 * dist2_lgrad(xnew, xold) );
}
/// Cached reported velocity
colvarvalue v_reported;
// Options for task_extended_lagrangian
/// Restraint center
colvarvalue xr;
/// Velocity of the restraint center
colvarvalue vr;
/// Mass of the restraint center
cvm::real ext_mass;
/// Restraint force constant
cvm::real ext_force_k;
/// Friction coefficient for Langevin extended dynamics
cvm::real ext_gamma;
/// Amplitude of Gaussian white noise for Langevin extended dynamics
cvm::real ext_sigma;
/// \brief Harmonic restraint force
colvarvalue fr;
/// \brief Jacobian force, when task_Jacobian_force is enabled
colvarvalue fj;
/// Cached reported system force
colvarvalue ft_reported;
public:
/// \brief Bias force; reset_bias_force() should be called before
/// the biases are updated
colvarvalue fb;
/// \brief Total \em applied force; fr (if task_extended_lagrangian
/// is defined), fb (if biases are applied) and the walls' forces
/// (if defined) contribute to it
colvarvalue f;
/// \brief Total force, as derived from the atomic trajectory;
/// should equal the total system force plus \link f \endlink
colvarvalue ft;
/// Period, if it is a constant
cvm::real period;
/// \brief Same as above, but also takes into account components
/// with a variable period, such as distanceZ
bool b_periodic;
/// \brief Expand the boundaries of multiples of width, to keep the
/// value always within range
bool expand_boundaries;
/// \brief Location of the lower boundary
colvarvalue lower_boundary;
/// \brief Location of the lower wall
colvarvalue lower_wall;
/// \brief Force constant for the lower boundary potential (|x-xb|^2)
cvm::real lower_wall_k;
/// \brief Whether this colvar has a hard lower boundary
bool hard_lower_boundary;
/// \brief Location of the upper boundary
colvarvalue upper_boundary;
/// \brief Location of the upper wall
colvarvalue upper_wall;
/// \brief Force constant for the upper boundary potential (|x-xb|^2)
cvm::real upper_wall_k;
/// \brief Whether this colvar has a hard upper boundary
bool hard_upper_boundary;
/// \brief Is the interval defined by the two boundaries periodic?
bool periodic_boundaries() const;
/// \brief Is the interval defined by the two boundaries periodic?
bool periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const;
/// Constructor
colvar(std::string const &conf);
/// Enable the specified task
int enable(colvar::task const &t);
/// Disable the specified task
void disable(colvar::task const &t);
/// Get ready for a run and possibly re-initialize internal data
void setup();
/// Destructor
~colvar();
/// \brief Calculate the colvar value and all the other requested
/// quantities
void calc();
/// \brief Calculate just the value, and store it in the argument
void calc_value(colvarvalue &xn);
/// Set the total biasing force to zero
void reset_bias_force();
/// Add to the total force from biases
void add_bias_force(colvarvalue const &force);
/// \brief Collect all forces on this colvar, integrate internal
/// equations of motion of internal degrees of freedom; see also
/// colvar::communicate_forces()
/// return colvar energy if extended Lagrandian active
cvm::real update();
/// \brief Communicate forces (previously calculated in
/// colvar::update()) to the external degrees of freedom
void communicate_forces();
/// \brief Enables and disables individual CVCs based on flags
int set_cvc_flags(std::vector<bool> const &flags);
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval
///
/// Handles correctly symmetries and periodic boundary conditions
void wrap(colvarvalue &x) const;
/// Read the analysis tasks
int parse_analysis(std::string const &conf);
/// Perform analysis tasks
void analyze();
/// Read the value from a collective variable trajectory file
std::istream & read_traj(std::istream &is);
/// Output formatted values to the trajectory file
std::ostream & write_traj(std::ostream &os);
/// Write a label to the trajectory file (comment line)
std::ostream & write_traj_label(std::ostream &os);
/// Read the collective variable from a restart file
std::istream & read_restart(std::istream &is);
/// Write the collective variable to a restart file
std::ostream & write_restart(std::ostream &os);
/// Write output files (if defined, e.g. in analysis mode)
int write_output_files();
protected:
/// Previous value (to calculate velocities during analysis)
colvarvalue x_old;
/// Time series of values and velocities used in correlation
/// functions
std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
/// Time series of values and velocities used in correlation
/// functions (pointers)x
std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;
/// Time series of values and velocities used in running averages
std::list< std::list<colvarvalue> > x_history;
/// Time series of values and velocities used in correlation
/// functions (pointers)x
std::list< std::list<colvarvalue> >::iterator x_history_p;
/// \brief Collective variable with which the correlation is
/// calculated (default: itself)
std::string acf_colvar_name;
/// Length of autocorrelation function (ACF)
size_t acf_length;
/// After how many steps the ACF starts
size_t acf_offset;
/// How many timesteps separate two ACF values
size_t acf_stride;
/// Number of frames for each ACF point
size_t acf_nframes;
/// Normalize the ACF to a maximum value of 1?
bool acf_normalize;
/// ACF values
std::vector<cvm::real> acf;
/// Name of the file to write the ACF
std::string acf_outfile;
/// Type of autocorrelation function (ACF)
enum acf_type_e {
/// Unset type
acf_notset,
/// Velocity ACF, scalar product between v(0) and v(t)
acf_vel,
/// Coordinate ACF, scalar product between x(0) and x(t)
acf_coor,
/// \brief Coordinate ACF, second order Legendre polynomial
/// between x(0) and x(t) (does not work with scalar numbers)
acf_p2coor
};
/// Type of autocorrelation function (ACF)
acf_type_e acf_type;
/// \brief Velocity ACF, scalar product between v(0) and v(t)
int calc_vel_acf(std::list<colvarvalue> &v_history,
colvarvalue const &v);
/// \brief Coordinate ACF, scalar product between x(0) and x(t)
/// (does not work with scalar numbers)
void calc_coor_acf(std::list<colvarvalue> &x_history,
colvarvalue const &x);
/// \brief Coordinate ACF, second order Legendre polynomial between
/// x(0) and x(t) (does not work with scalar numbers)
void calc_p2coor_acf(std::list<colvarvalue> &x_history,
colvarvalue const &x);
/// Calculate the auto-correlation function (ACF)
int calc_acf();
/// Save the ACF to a file
void write_acf(std::ostream &os);
/// Length of running average series
size_t runave_length;
/// Timesteps to skip between two values in the running average series
size_t runave_stride;
/// Name of the file to write the running average
cvm::ofstream runave_os;
/// Current value of the running average
colvarvalue runave;
/// Current value of the square deviation from the running average
cvm::real runave_variance;
/// Calculate the running average and its standard deviation
void calc_runave();
/// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
cvm::real kinetic_energy;
cvm::real potential_energy;
public:
// collective variable component base class
class cvc;
// currently available collective variable components
// scalar colvar components
class distance;
class distance_z;
class distance_xy;
class distance_inv;
class distance_pairs;
class angle;
class dihedral;
class coordnum;
class selfcoordnum;
class h_bond;
class rmsd;
class orientation_angle;
class orientation_proj;
class tilt;
class spin_angle;
class gyration;
class inertia;
class inertia_z;
class eigenvector;
class alpha_dihedrals;
class alpha_angles;
class dihedPC;
// non-scalar components
class distance_vec;
class distance_dir;
class cartesian;
class orientation;
protected:
/// \brief Array of \link cvc \endlink objects
std::vector<cvc *> cvcs;
+ /// \brief Flags to enable or disable cvcs at next colvar evaluation
+ std::vector<bool> cvc_flags;
+
/// \brief Initialize the sorted list of atom IDs for atoms involved
/// in all cvcs (called when enabling task_collect_gradients)
void build_atom_list(void);
private:
/// Name of scripted function to be used
std::string scripted_function;
/// Current cvc values in the order requested by script
/// when using scriptedFunction
std::vector<const colvarvalue *> sorted_cvc_values;
public:
/// \brief Sorted array of (zero-based) IDs for all atoms involved
std::vector<int> atom_ids;
/// \brief Array of atomic gradients collected from all cvcs
/// with appropriate components, rotations etc.
/// For scalar variables only!
std::vector<cvm::rvector> atomic_gradients;
inline size_t n_components() const {
return cvcs.size();
}
};
inline colvarvalue const & colvar::value() const
{
return x_reported;
}
inline colvarvalue const & colvar::actual_value() const
{
return x;
}
inline colvarvalue const & colvar::velocity() const
{
return v_reported;
}
inline colvarvalue const & colvar::system_force() const
{
return ft_reported;
}
inline void colvar::add_bias_force(colvarvalue const &force)
{
fb += force;
}
inline void colvar::reset_bias_force() {
fb.type(value());
fb.reset();
}
#endif
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 6cc667b51..e2e660b98 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -1,244 +1,172 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvargrid.h"
colvar_grid_count::colvar_grid_count()
: colvar_grid<size_t>()
{
mult = 1;
}
colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
size_t const &def_count)
: colvar_grid<size_t>(nx_i, def_count, 1)
{}
colvar_grid_count::colvar_grid_count(std::vector<colvar *> &colvars,
size_t const &def_count)
: colvar_grid<size_t>(colvars, def_count, 1)
{}
-std::istream & colvar_grid_count::read_restart(std::istream &is)
-{
- size_t const start_pos = is.tellg();
- std::string key, conf;
- if ((is >> key) && (key == std::string("grid_parameters"))) {
- is.seekg(start_pos, std::ios::beg);
- is >> colvarparse::read_block("grid_parameters", conf);
- parse_params(conf);
- } else {
- cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
- is.seekg(start_pos, std::ios::beg);
- }
- read_raw(is);
- return is;
-}
-
-std::ostream & colvar_grid_count::write_restart(std::ostream &os)
-{
- write_params(os);
- write_raw(os);
- return os;
-}
-
-
-
colvar_grid_scalar::colvar_grid_scalar()
: colvar_grid<cvm::real>(), samples(NULL), grad(NULL)
{}
colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
: colvar_grid<cvm::real>(g), samples(NULL), grad(NULL)
{
grad = new cvm::real[nd];
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
: colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
{
grad = new cvm::real[nd];
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
: colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
{
grad = new cvm::real[nd];
}
colvar_grid_scalar::~colvar_grid_scalar()
{
if (grad) {
delete [] grad;
grad = NULL;
}
}
-std::istream & colvar_grid_scalar::read_restart(std::istream &is)
-{
- size_t const start_pos = is.tellg();
- std::string key, conf;
- if ((is >> key) && (key == std::string("grid_parameters"))) {
- is.seekg(start_pos, std::ios::beg);
- is >> colvarparse::read_block("grid_parameters", conf);
- parse_params(conf);
- } else {
- cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
- is.seekg(start_pos, std::ios::beg);
- }
- read_raw(is);
- return is;
-}
-
-std::ostream & colvar_grid_scalar::write_restart(std::ostream &os)
-{
- write_params(os);
- write_raw(os);
- return os;
-}
-
-
cvm::real colvar_grid_scalar::maximum_value() const
{
cvm::real max = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] > max) max = data[i];
}
return max;
}
cvm::real colvar_grid_scalar::minimum_value() const
{
cvm::real min = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] < min) min = data[i];
}
return min;
}
cvm::real colvar_grid_scalar::integral() const
{
cvm::real sum = 0.0;
for (size_t i = 0; i < nt; i++) {
sum += data[i];
}
cvm::real bin_volume = 1.0;
for (size_t id = 0; id < widths.size(); id++) {
bin_volume *= widths[id];
}
return bin_volume * sum;
}
cvm::real colvar_grid_scalar::entropy() const
{
cvm::real sum = 0.0;
for (size_t i = 0; i < nt; i++) {
sum += -1.0 * data[i] * std::log(data[i]);
}
cvm::real bin_volume = 1.0;
for (size_t id = 0; id < widths.size(); id++) {
bin_volume *= widths[id];
}
return bin_volume * sum;
}
colvar_grid_gradient::colvar_grid_gradient()
: colvar_grid<cvm::real>(), samples(NULL)
{}
colvar_grid_gradient::colvar_grid_gradient(std::vector<int> const &nx_i)
: colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()), samples(NULL)
{}
colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars)
: colvar_grid<cvm::real>(colvars, 0.0, colvars.size()), samples(NULL)
{}
-std::istream & colvar_grid_gradient::read_restart(std::istream &is)
-{
- size_t const start_pos = is.tellg();
- std::string key, conf;
- if ((is >> key) && (key == std::string("grid_parameters"))) {
- is.seekg(start_pos, std::ios::beg);
- is >> colvarparse::read_block("grid_parameters", conf);
- parse_params(conf);
- } else {
- cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
- is.seekg(start_pos, std::ios::beg);
- }
- read_raw(is);
- return is;
-}
-
-std::ostream & colvar_grid_gradient::write_restart(std::ostream &os)
-{
- write_params(os);
- write_raw(os);
- return os;
-}
-
void colvar_grid_gradient::write_1D_integral(std::ostream &os)
{
cvm::real bin, min, integral;
std::vector<cvm::real> int_vals;
os << "# xi A(xi)\n";
if ( cv.size() != 1 ) {
cvm::fatal_error("Cannot write integral for multi-dimensional gradient grids.");
}
integral = 0.0;
int_vals.push_back( 0.0 );
bin = 0.0;
min = 0.0;
// correction for periodic colvars, so that the PMF is periodic
cvm::real corr;
if ( periodic[0] ) {
corr = average();
} else {
corr = 0.0;
}
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix), bin += 1.0 ) {
if (samples) {
size_t const samples_here = samples->value(ix);
if (samples_here)
integral += (value(ix) / cvm::real(samples_here) - corr) * cv[0]->width;
} else {
integral += (value(ix) - corr) * cv[0]->width;
}
if ( integral < min ) min = integral;
int_vals.push_back( integral );
}
bin = 0.0;
for ( int i = 0; i < nx[0]; i++, bin += 1.0 ) {
os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
<< std::setw(cvm::cv_width)
<< std::setprecision(cvm::cv_prec)
<< int_vals[i] - min << "\n";
}
os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
<< std::setw(cvm::cv_width)
<< std::setprecision(cvm::cv_prec)
<< int_vals[nx[0]] - min << "\n";
return;
}
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 29d21feb3..ac192ed0e 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1,1384 +1,1392 @@
/// -*- c++ -*-
#ifndef COLVARGRID_H
#define COLVARGRID_H
#include <iostream>
#include <iomanip>
#include <cmath>
#include "colvar.h"
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
/// \brief Grid of values of a function of several collective
/// variables \param T The data type
///
/// Only scalar colvars supported so far: vector colvars are treated as arrays
template <class T> class colvar_grid : public colvarparse {
protected:
/// Number of dimensions
size_t nd;
/// Number of points along each dimension
std::vector<int> nx;
/// Cumulative number of points along each dimension
std::vector<int> nxc;
/// \brief Multiplicity of each datum (allow the binning of
/// non-scalar types such as atomic gradients)
size_t mult;
/// Total number of grid points
size_t nt;
/// Low-level array of values
std::vector<T> data;
/// Newly read data (used for count grids, when adding several grids read from disk)
std::vector<size_t> new_data;
/// Colvars collected in this grid
std::vector<colvar *> cv;
/// Do we request actual value (for extended-system colvars)?
std::vector<bool> actual_value;
/// Get the low-level index corresponding to an index
inline size_t address(std::vector<int> const &ix) const
{
size_t addr = 0;
for (size_t i = 0; i < nd; i++) {
addr += ix[i]*nxc[i];
if (cvm::debug()) {
if (ix[i] >= nx[i]) {
cvm::error("Error: exceeding bounds in colvar_grid.\n", BUG_ERROR);
return 0;
}
}
}
return addr;
}
public:
/// Lower boundaries of the colvars in this grid
std::vector<colvarvalue> lower_boundaries;
/// Upper boundaries of the colvars in this grid
std::vector<colvarvalue> upper_boundaries;
/// Whether some colvars are periodic
std::vector<bool> periodic;
/// Whether some colvars have hard lower boundaries
std::vector<bool> hard_lower_boundaries;
/// Whether some colvars have hard upper boundaries
std::vector<bool> hard_upper_boundaries;
/// Widths of the colvars in this grid
std::vector<cvm::real> widths;
/// True if this is a count grid related to another grid of data
bool has_parent_data;
/// Whether this grid has been filled with data or is still empty
bool has_data;
/// Return the number of colvars
inline size_t number_of_colvars() const
{
return nd;
}
/// Return the number of points in the i-th direction, if provided, or
/// the total number
inline size_t number_of_points(int const icv = -1) const
{
if (icv < 0) {
return nt;
} else {
return nx[icv];
}
}
/// Get the sizes in each direction
inline std::vector<int> const & sizes() const
{
return nx;
}
/// Set the sizes in each direction
inline void set_sizes(std::vector<int> const &new_sizes)
{
nx = new_sizes;
}
/// Return the multiplicity of the type used
inline size_t multiplicity() const
{
return mult;
}
/// \brief Allocate data
int setup(std::vector<int> const &nx_i,
T const &t = T(),
size_t const &mult_i = 1)
{
if (cvm::debug()) {
cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+
", dimensions = "+cvm::to_str(nx_i)+".\n");
}
mult = mult_i;
data.clear();
nx = nx_i;
nd = nx.size();
nxc.resize(nd);
// setup dimensions
nt = mult;
for (int i = nd-1; i >= 0; i--) {
if (nx[i] <= 0) {
cvm::error("Error: providing an invalid number of grid points, "+
cvm::to_str(nx[i])+".\n", BUG_ERROR);
return COLVARS_ERROR;
}
nxc[i] = nt;
nt *= nx[i];
}
if (cvm::debug()) {
cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n");
}
data.reserve(nt);
data.assign(nt, t);
return COLVARS_OK;
}
/// \brief Allocate data (allow initialization also after construction)
int setup()
{
return setup(this->nx, T(), this->mult);
}
/// \brief Reset data (in case the grid is being reused)
void reset(T const &t = T())
{
data.assign(nt, t);
}
/// Default constructor
colvar_grid() : has_data(false)
{
save_delimiters = false;
nd = nt = 0;
mult = 1;
this->setup();
}
/// Destructor
virtual ~colvar_grid()
{}
/// \brief "Almost copy-constructor": only copies configuration
/// parameters from another grid, but doesn't reallocate stuff;
/// setup() must be called after that;
colvar_grid(colvar_grid<T> const &g) : nd(g.nd),
nx(g.nx),
mult(g.mult),
data(),
cv(g.cv),
actual_value(g.actual_value),
lower_boundaries(g.lower_boundaries),
upper_boundaries(g.upper_boundaries),
periodic(g.periodic),
hard_lower_boundaries(g.hard_lower_boundaries),
hard_upper_boundaries(g.hard_upper_boundaries),
widths(g.widths),
has_data(false)
{
save_delimiters = false;
}
/// \brief Constructor from explicit grid sizes \param nx_i Number
/// of grid points along each dimension \param t Initial value for
/// the function at each point (optional) \param mult_i Multiplicity
/// of each value
colvar_grid(std::vector<int> const &nx_i,
T const &t = T(),
size_t mult_i = 1)
: has_data(false)
{
save_delimiters = false;
this->setup(nx_i, t, mult_i);
}
/// \brief Constructor from a vector of colvars
colvar_grid(std::vector<colvar *> const &colvars,
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
: has_data(false)
{
save_delimiters = false;
this->init_from_colvars(colvars, t, mult_i, margin);
}
int init_from_colvars(std::vector<colvar *> const &colvars,
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
{
if (cvm::debug()) {
cvm::log("Reading grid configuration from collective variables.\n");
}
cv = colvars;
nd = colvars.size();
mult = mult_i;
size_t i;
for (i = 0; i < cv.size(); i++) {
if (!cv[i]->tasks[colvar::task_lower_boundary] ||
!cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::error("Tried to initialize a grid on a "
"variable with undefined boundaries.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
}
if (cvm::debug()) {
cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+
" collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
}
for (i = 0; i < cv.size(); i++) {
if (cv[i]->value().type() != colvarvalue::type_scalar) {
cvm::error("Colvar grids can only be automatically "
"constructed for scalar variables. "
"ABF and histogram can not be used; metadynamics "
"can be used with useGrids disabled.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
if (cv[i]->width <= 0.0) {
cvm::error("Tried to initialize a grid on a "
"variable with negative or zero width.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
widths.push_back(cv[i]->width);
hard_lower_boundaries.push_back(cv[i]->hard_lower_boundary);
hard_upper_boundaries.push_back(cv[i]->hard_upper_boundary);
periodic.push_back(cv[i]->periodic_boundaries());
// By default, get reported colvar value (for extended Lagrangian colvars)
actual_value.push_back(false);
// except if a colvar is specified twice in a row
// then the first instance is the actual value
if (i > 0 && cv[i-1] == cv[i]) {
actual_value[i-1] = true;
}
if (margin) {
if (periodic[i]) {
// Shift the grid by half the bin width (values at edges instead of center of bins)
lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back(cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
} else {
// Make this grid larger by one bin width
lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back(cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
}
} else {
lower_boundaries.push_back(cv[i]->lower_boundary);
upper_boundaries.push_back(cv[i]->upper_boundary);
}
}
this->init_from_boundaries();
return this->setup();
}
int init_from_boundaries(T const &t = T(),
size_t const &mult_i = 1)
{
if (cvm::debug()) {
cvm::log("Configuring grid dimensions from colvars boundaries.\n");
}
// these will have to be updated
nx.clear();
nxc.clear();
nt = 0;
for (size_t i = 0; i < lower_boundaries.size(); i++) {
cvm::real nbins = ( upper_boundaries[i].real_value -
lower_boundaries[i].real_value ) / widths[i];
int nbins_round = (int)(nbins+0.5);
if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
cvm::log("Warning: grid interval("+
cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
") is not commensurate to its bin width("+
cvm::to_str(widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
upper_boundaries[i].real_value = lower_boundaries[i].real_value +
(nbins_round * widths[i]);
}
if (cvm::debug())
cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+
" for the colvar no. "+cvm::to_str(i+1)+".\n");
nx.push_back(nbins_round);
}
return COLVARS_OK;
}
/// Wrap an index vector around periodic boundary conditions
/// also checks validity of non-periodic indices
inline void wrap(std::vector<int> & ix) const
{
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) {
ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
} else {
if (ix[i] < 0 || ix[i] >= nx[i]) {
cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: "
+ cvm::to_str(ix), BUG_ERROR);
return;
}
}
}
}
/// \brief Report the bin corresponding to the current value of variable i
inline int current_bin_scalar(int const i) const
{
return value_to_bin_scalar(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
}
/// \brief Use the lower boundary and the width to report which bin
/// the provided value is in
inline int value_to_bin_scalar(colvarvalue const &value, const int i) const
{
return (int) std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
}
/// \brief Same as the standard version, but uses another grid definition
inline int value_to_bin_scalar(colvarvalue const &value,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return (int) std::floor( (value.real_value - new_offset.real_value) / new_width );
}
/// \brief Use the two boundaries and the width to report the
/// central value corresponding to a bin index
inline colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const
{
return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
}
/// \brief Same as the standard version, but uses different parameters
inline colvarvalue bin_to_value_scalar(int const &i_bin,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return new_offset.real_value + new_width * (0.5 + i_bin);
}
/// Set the value at the point with index ix
inline void set_value(std::vector<int> const &ix,
T const &t,
size_t const &imult = 0)
{
data[this->address(ix)+imult] = t;
has_data = true;
}
/// \brief Get the change from this to other_grid
/// and store the result in this.
/// this_grid := other_grid - this_grid
/// Grids must have the same dimensions.
void delta_grid(colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error("Error: trying to subtract two grids with "
"different multiplicity.\n");
return;
}
if (other_grid.data.size() != this->data.size()) {
cvm::error("Error: trying to subtract two grids with "
"different size.\n");
return;
}
for (size_t i = 0; i < data.size(); i++) {
data[i] = other_grid.data[i] - data[i];
}
has_data = true;
}
/// \brief Copy data from another grid of the same type, AND
/// identical definition (boundaries, widths)
/// Added for shared ABF.
void copy_grid(colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error("Error: trying to copy two grids with "
"different multiplicity.\n");
return;
}
if (other_grid.data.size() != this->data.size()) {
cvm::error("Error: trying to copy two grids with "
"different size.\n");
return;
}
for (size_t i = 0; i < data.size(); i++) {
data[i] = other_grid.data[i];
}
has_data = true;
}
/// \brief Extract the grid data as they are represented in memory.
/// Put the results in "out_data".
void raw_data_out(T* out_data) const
{
for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
}
/// \brief Input the data as they are represented in memory.
void raw_data_in(const T* in_data)
{
for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
has_data = true;
}
/// \brief Size of the data as they are represented in memory.
size_t raw_data_num() const { return data.size(); }
/// \brief Get the binned value indexed by ix, or the first of them
/// if the multiplicity is larger than 1
inline T const & value(std::vector<int> const &ix,
size_t const &imult = 0) const
{
return data[this->address(ix) + imult];
}
/// \brief Add a constant to all elements (fast loop)
inline void add_constant(T const &t)
{
for (size_t i = 0; i < nt; i++)
data[i] += t;
has_data = true;
}
/// \brief Multiply all elements by a scalar constant (fast loop)
inline void multiply_constant(cvm::real const &a)
{
for (size_t i = 0; i < nt; i++)
data[i] *= a;
}
/// \brief Get the bin indices corresponding to the provided values of
/// the colvars
inline std::vector<int> const get_colvars_index(std::vector<colvarvalue> const &values) const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = value_to_bin_scalar(values[i], i);
}
return index;
}
/// \brief Get the bin indices corresponding to the current values
/// of the colvars
inline std::vector<int> const get_colvars_index() const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = current_bin_scalar(i);
}
return index;
}
/// \brief Get the minimal distance (in number of bins) from the
/// boundaries; a negative number is returned if the given point is
/// off-grid
inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
bool skip_hard_boundaries = false)
{
cvm::real minimum = 1.0E+16;
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) continue;
cvm::real dl = std::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
cvm::real du = std::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
if (values[i].real_value < lower_boundaries[i])
dl *= -1.0;
if (values[i].real_value > upper_boundaries[i])
du *= -1.0;
if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
minimum = dl;
if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
minimum = du;
}
return minimum;
}
/// \brief Add data from another grid of the same type
///
/// Note: this function maps other_grid inside this one regardless
/// of whether it fits or not.
void map_grid(colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error("Error: trying to merge two grids with values of "
"different multiplicity.\n");
return;
}
std::vector<colvarvalue> const &gb = this->lower_boundaries;
std::vector<cvm::real> const &gw = this->widths;
std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
std::vector<cvm::real> const &ogw = other_grid.widths;
std::vector<int> ix = this->new_index();
std::vector<int> oix = other_grid.new_index();
if (cvm::debug())
cvm::log("Remapping grid...\n");
for ( ; this->index_ok(ix); this->incr(ix)) {
for (size_t i = 0; i < nd; i++) {
oix[i] =
value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
ogb[i],
ogw[i]);
}
if (! other_grid.index_ok(oix)) {
continue;
}
for (size_t im = 0; im < mult; im++) {
this->set_value(ix, other_grid.value(oix, im), im);
}
}
has_data = true;
if (cvm::debug())
cvm::log("Remapping done.\n");
}
/// \brief Add data from another grid of the same type, AND
/// identical definition (boundaries, widths)
void add_grid(colvar_grid<T> const &other_grid,
cvm::real scale_factor = 1.0)
{
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error("Error: trying to sum togetehr two grids with values of "
"different multiplicity.\n");
return;
}
if (scale_factor != 1.0)
for (size_t i = 0; i < data.size(); i++) {
data[i] += scale_factor * other_grid.data[i];
}
else
// skip multiplication if possible
for (size_t i = 0; i < data.size(); i++) {
data[i] += other_grid.data[i];
}
has_data = true;
}
/// \brief Return the value suitable for output purposes (so that it
/// may be rescaled or manipulated without changing it permanently)
virtual inline T value_output(std::vector<int> const &ix,
size_t const &imult = 0)
{
return value(ix, imult);
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (the two may be different,
/// e.g. when using colvar_grid_count)
virtual inline void value_input(std::vector<int> const &ix,
T const &t,
size_t const &imult = 0,
bool add = false)
{
if ( add )
data[address(ix) + imult] += t;
else
data[address(ix) + imult] = t;
has_data = true;
}
// /// Get the pointer to the binned value indexed by ix
// inline T const *value_p (std::vector<int> const &ix)
// {
// return &(data[address (ix)]);
// }
/// \brief Get the index corresponding to the "first" bin, to be
/// used as the initial value for an index in looping
inline std::vector<int> const new_index() const
{
return std::vector<int> (nd, 0);
}
/// \brief Check that the index is within range in each of the
/// dimensions
inline bool index_ok(std::vector<int> const &ix) const
{
for (size_t i = 0; i < nd; i++) {
if ( (ix[i] < 0) || (ix[i] >= int(nx[i])) )
return false;
}
return true;
}
/// \brief Increment the index, in a way that will make it loop over
/// the whole nd-dimensional array
inline void incr(std::vector<int> &ix) const
{
for (int i = ix.size()-1; i >= 0; i--) {
ix[i]++;
if (ix[i] >= nx[i]) {
if (i > 0) {
ix[i] = 0;
continue;
} else {
// this is the last iteration, a non-valid index is being
// set for the outer index, which will be caught by
// index_ok()
ix[0] = nx[0];
return;
}
} else {
return;
}
}
}
/// \brief Write the grid parameters (number of colvars, boundaries, width and number of points)
std::ostream & write_params(std::ostream &os)
{
size_t i;
os << "grid_parameters {\n n_colvars " << nd << "\n";
os << " lower_boundaries ";
for (i = 0; i < nd; i++)
os << " " << lower_boundaries[i];
os << "\n";
os << " upper_boundaries ";
for (i = 0; i < nd; i++)
os << " " << upper_boundaries[i];
os << "\n";
os << " widths ";
for (i = 0; i < nd; i++)
os << " " << widths[i];
os << "\n";
os << " sizes ";
for (i = 0; i < nd; i++)
os << " " << nx[i];
os << "\n";
os << "}\n";
return os;
}
/// Read a grid definition from a config string
int parse_params(std::string const &conf)
{
if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");
std::vector<int> old_nx = nx;
std::vector<colvarvalue> old_lb = lower_boundaries;
{
size_t nd_in = 0;
colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
if (nd_in != nd) {
cvm::error("Error: trying to read data for a grid "
"that contains a different number of colvars ("+
cvm::to_str(nd_in)+") than the grid defined "
"in the configuration file("+cvm::to_str(nd)+
").\n");
return COLVARS_ERROR;
}
}
colvarparse::get_keyval(conf, "lower_boundaries",
lower_boundaries, lower_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval(conf, "upper_boundaries",
upper_boundaries, upper_boundaries, colvarparse::parse_silent);
// support also camel case
colvarparse::get_keyval(conf, "lowerBoundaries",
lower_boundaries, lower_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval(conf, "upperBoundaries",
upper_boundaries, upper_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);
if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
if (! actual_value.size()) actual_value.assign(nd, false);
if (! periodic.size()) periodic.assign(nd, false);
if (! widths.size()) widths.assign(nd, 1.0);
bool new_params = false;
if (old_nx.size()) {
for (size_t i = 0; i < nd; i++) {
if ( (old_nx[i] != nx[i]) ||
(std::sqrt(cv[i]->dist2(old_lb[i],
lower_boundaries[i])) > 1.0E-10) ) {
new_params = true;
}
}
} else {
new_params = true;
}
// reallocate the array in case the grid params have just changed
if (new_params) {
init_from_boundaries();
// data.resize(0); // no longer needed: setup calls clear()
return this->setup(nx, T(), mult);
}
return COLVARS_OK;
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with the current setting of the
/// colvars
void check_consistency()
{
for (size_t i = 0; i < nd; i++) {
if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
lower_boundaries[i])) > 1.0E-10) ||
(std::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
upper_boundaries[i])) > 1.0E-10) ||
(std::sqrt(cv[i]->dist2(cv[i]->width,
widths[i])) > 1.0E-10) ) {
cvm::error("Error: restart information for a grid is "
"inconsistent with that of its colvars.\n");
return;
}
}
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with that of another grid
void check_consistency(colvar_grid<T> const &other_grid)
{
for (size_t i = 0; i < nd; i++) {
// we skip dist2(), because periodicities and the like should
// matter: boundaries should be EXACTLY the same (otherwise,
// map_grid() should be used)
if ( (std::fabs(other_grid.lower_boundaries[i] -
lower_boundaries[i]) > 1.0E-10) ||
(std::fabs(other_grid.upper_boundaries[i] -
upper_boundaries[i]) > 1.0E-10) ||
(std::fabs(other_grid.widths[i] -
widths[i]) > 1.0E-10) ||
(data.size() != other_grid.data.size()) ) {
cvm::error("Error: inconsistency between "
"two grids that are supposed to be equal, "
"aside from the data stored.\n");
return;
}
}
}
+ /// \brief Read grid entry in restart file
+ std::istream & read_restart(std::istream &is)
+ {
+ size_t const start_pos = is.tellg();
+ std::string key, conf;
+ if ((is >> key) && (key == std::string("grid_parameters"))) {
+ is.seekg(start_pos, std::ios::beg);
+ is >> colvarparse::read_block("grid_parameters", conf);
+ parse_params(conf);
+ } else {
+ cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
+ is.seekg(start_pos, std::ios::beg);
+ }
+ read_raw(is);
+ return is;
+ }
+
+ /// \brief Write grid entry in restart file
+ std::ostream & write_restart(std::ostream &os)
+ {
+ write_params(os);
+ write_raw(os);
+ return os;
+ }
+
+
/// \brief Write the grid data without labels, as they are
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_raw(std::ostream &os,
size_t const buf_size = 3)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
std::vector<int> ix = new_index();
size_t count = 0;
for ( ; index_ok(ix); incr(ix)) {
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< value_output(ix, imult);
if (((++count) % buf_size) == 0)
os << "\n";
}
}
// write a final newline only if buffer is not empty
if ((count % buf_size) != 0)
os << "\n";
return os;
}
/// \brief Read data written by colvar_grid::write_raw()
std::istream & read_raw(std::istream &is)
{
size_t const start_pos = is.tellg();
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
for (size_t imult = 0; imult < mult; imult++) {
T new_value;
if (is >> new_value) {
value_input(ix, new_value, imult);
} else {
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
cvm::error("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
return is;
}
}
}
has_data = true;
return is;
}
/// \brief Write the grid in a format which is both human readable
/// and suitable for visualization e.g. with gnuplot
void write_multicol(std::ostream &os)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
os << std::setw(2) << "# " << nd << "\n";
for (size_t i = 0; i < nd; i++) {
os << "# "
<< std::setw(10) << lower_boundaries[i]
<< std::setw(10) << widths[i]
<< std::setw(10) << nx[i] << " "
<< periodic[i] << "\n";
}
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
if (ix.back() == 0) {
// if the last index is 0, add a new line to mark the new record
os << "\n";
}
for (size_t i = 0; i < nd; i++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< bin_to_value_scalar(ix[i], i);
}
os << " ";
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< value_output(ix, imult);
}
os << "\n";
}
}
/// \brief Read a grid written by colvar_grid::write_multicol()
/// Adding data if add is true, replacing if false
std::istream & read_multicol(std::istream &is, bool add = false)
{
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
std::string hash;
cvm::real lower, width, x;
size_t n, periodic;
bool remap;
std::vector<T> new_value;
std::vector<int> nx_read;
std::vector<int> bin;
if ( cv.size() != nd ) {
cvm::error("Cannot read grid file: missing reference to colvars.");
return is;
}
if ( !(is >> hash) || (hash != "#") ) {
cvm::error("Error reading grid at position "+
cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
return is;
}
is >> n;
if ( n != nd ) {
cvm::error("Error reading grid: wrong number of collective variables.\n");
return is;
}
nx_read.resize(n);
bin.resize(n);
new_value.resize(mult);
if (this->has_parent_data && add) {
new_data.resize(data.size());
}
remap = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> hash) || (hash != "#") ) {
cvm::error("Error reading grid at position "+
cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
return is;
}
is >> lower >> width >> nx_read[i] >> periodic;
if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
(std::fabs(width - widths[i] ) > 1.0e-10) ||
(nx_read[i] != nx[i]) ) {
cvm::log("Warning: reading from different grid definition (colvar "
+ cvm::to_str(i+1) + "); remapping data on new grid.\n");
remap = true;
}
}
if ( remap ) {
// re-grid data
while (is.good()) {
bool end_of_file = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> x) ) end_of_file = true;
bin[i] = value_to_bin_scalar(x, i);
}
if (end_of_file) break;
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
}
if ( index_ok(bin) ) {
for (size_t imult = 0; imult < mult; imult++) {
value_input(bin, new_value[imult], imult, add);
}
}
}
} else {
// do not re-grid the data but assume the same grid is used
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
for (size_t i = 0; i < nd; i++ ) {
is >> x;
}
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
value_input(ix, new_value[imult], imult, add);
}
}
}
has_data = true;
return is;
}
/// \brief Write the grid data without labels, as they are
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_opendx(std::ostream &os)
{
// write the header
os << "object 1 class gridpositions counts";
int icv;
for (icv = 0; icv < number_of_colvars(); icv++) {
os << " " << number_of_points(icv);
}
os << "\n";
os << "origin";
for (icv = 0; icv < number_of_colvars(); icv++) {
os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]);
}
os << "\n";
for (icv = 0; icv < number_of_colvars(); icv++) {
os << "delta";
for (size_t icv2 = 0; icv2 < number_of_colvars(); icv2++) {
if (icv == icv2) os << " " << widths[icv];
else os << " " << 0.0;
}
os << "\n";
}
os << "object 2 class gridconnections counts";
for (icv = 0; icv < number_of_colvars(); icv++) {
os << " " << number_of_points(icv);
}
os << "\n";
os << "object 3 class array type double rank 0 items "
<< number_of_points() << " data follows\n";
write_raw(os);
os << "object \"collective variables scalar field\" class field\n";
return os;
}
};
/// \brief Colvar_grid derived class to hold counters in discrete
/// n-dim colvar space
class colvar_grid_count : public colvar_grid<size_t>
{
public:
/// Default constructor
colvar_grid_count();
/// Destructor
virtual inline ~colvar_grid_count()
{}
/// Constructor
colvar_grid_count(std::vector<int> const &nx_i,
size_t const &def_count = 0);
/// Constructor from a vector of colvars
colvar_grid_count(std::vector<colvar *> &colvars,
size_t const &def_count = 0);
/// Increment the counter at given position
inline void incr_count(std::vector<int> const &ix)
{
++(data[this->address(ix)]);
}
/// \brief Get the binned count indexed by ix from the newly read data
inline size_t const & new_count(std::vector<int> const &ix,
size_t const &imult = 0)
{
return new_data[address(ix) + imult];
}
- /// \brief Read the grid from a restart
- std::istream & read_restart(std::istream &is);
-
- /// \brief Write the grid to a restart
- std::ostream & write_restart(std::ostream &os);
-
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input(std::vector<int> const &ix,
size_t const &t,
size_t const &imult = 0,
bool add = false)
{
if (add) {
data[address(ix)] += t;
if (this->has_parent_data) {
// save newly read data for inputting parent grid
new_data[address(ix)] = t;
}
} else {
data[address(ix)] = t;
}
has_data = true;
}
};
/// Class for accumulating a scalar function on a grid
class colvar_grid_scalar : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the associated sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_scalar();
/// Copy constructor (needed because of the grad pointer)
colvar_grid_scalar(colvar_grid_scalar const &g);
/// Destructor
~colvar_grid_scalar();
/// Constructor from specific sizes arrays
colvar_grid_scalar(std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_scalar(std::vector<colvar *> &colvars,
bool margin = 0);
/// Accumulate the value
inline void acc_value(std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0)
{
// only legal value of imult here is 0
data[address(ix)] += new_value;
if (samples)
samples->incr_count(ix);
has_data = true;
}
/// Return the gradient of the scalar field from finite differences
inline const cvm::real * gradient_finite_diff( const std::vector<int> &ix0 )
{
cvm::real A0, A1;
std::vector<int> ix;
if (nd != 2) {
cvm::error("Finite differences available in dimension 2 only.");
return grad;
}
for (unsigned int n = 0; n < nd; n++) {
ix = ix0;
A0 = data[address(ix)];
ix[n]++; wrap(ix);
A1 = data[address(ix)];
ix[1-n]++; wrap(ix);
A1 += data[address(ix)];
ix[n]--; wrap(ix);
A0 += data[address(ix)];
grad[n] = 0.5 * (A1 - A0) / widths[n];
}
return grad;
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual cvm::real value_output(std::vector<int> const &ix,
size_t const &imult = 0)
{
if (imult > 0) {
cvm::error("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
return 0.;
}
if (samples) {
return (samples->value(ix) > 0) ?
(data[address(ix)] / cvm::real(samples->value(ix))) :
0.0;
} else {
return data[address(ix)];
}
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual void value_input(std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (imult > 0) {
cvm::error("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
return;
}
if (add) {
if (samples)
data[address(ix)] += new_value * samples->new_count(ix);
else
data[address(ix)] += new_value;
} else {
if (samples)
data[address(ix)] = new_value * samples->value(ix);
else
data[address(ix)] = new_value;
}
has_data = true;
}
- /// \brief Read the grid from a restart
- std::istream & read_restart(std::istream &is);
-
- /// \brief Write the grid to a restart
- std::ostream & write_restart(std::ostream &os);
-
/// \brief Return the highest value
cvm::real maximum_value() const;
/// \brief Return the lowest value
cvm::real minimum_value() const;
/// \brief Calculates the integral of the map (uses widths if they are defined)
cvm::real integral() const;
/// \brief Assuming that the map is a normalized probability density,
/// calculates the entropy (uses widths if they are defined)
cvm::real entropy() const;
private:
// gradient
cvm::real * grad;
};
/// Class for accumulating the gradient of a scalar function on a grid
class colvar_grid_gradient : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_gradient();
/// Destructor
virtual inline ~colvar_grid_gradient()
{}
/// Constructor from specific sizes arrays
colvar_grid_gradient(std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_gradient(std::vector<colvar *> &colvars);
/// \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++) {
data[address(ix) + imult] += grads[imult];
}
if (samples)
samples->incr_count(ix);
}
/// \brief Accumulate the gradient based on the force (i.e. sums the
/// opposite of the force)
inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
for (size_t imult = 0; imult < mult; imult++) {
data[address(ix) + imult] -= forces[imult];
}
if (samples)
samples->incr_count(ix);
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual inline cvm::real value_output(std::vector<int> const &ix,
size_t const &imult = 0)
{
if (samples)
return (samples->value(ix) > 0) ?
(data[address(ix) + imult] / cvm::real(samples->value(ix))) :
0.0;
else
return data[address(ix) + imult];
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input(std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (add) {
if (samples)
data[address(ix) + imult] += new_value * samples->new_count(ix);
else
data[address(ix) + imult] += new_value;
} else {
if (samples)
data[address(ix) + imult] = new_value * samples->value(ix);
else
data[address(ix) + imult] = new_value;
}
has_data = true;
}
- /// \brief Read the grid from a restart
- std::istream & read_restart(std::istream &is);
-
- /// \brief Write the grid to a restart
- std::ostream & write_restart(std::ostream &os);
-
/// Compute and return average value for a 1D gradient grid
inline cvm::real average()
{
size_t n = 0;
if (nd != 1 || nx[0] == 0) {
return 0.0;
}
cvm::real sum = 0.0;
std::vector<int> ix = new_index();
if (samples) {
for ( ; index_ok(ix); incr(ix)) {
if ( (n = samples->value(ix)) )
sum += value(ix) / n;
}
} else {
for ( ; index_ok(ix); incr(ix)) {
sum += value(ix);
}
}
return (sum / cvm::real(nx[0]));
}
/// \brief If the grid is 1-dimensional, integrate it and write the
/// integral to a file
void write_1D_integral(std::ostream &os);
};
#endif
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index e1d25b1d9..8fe92cea1 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -1,625 +1,625 @@
/// -*- c++ -*-
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2015-07-21"
+#define COLVARS_VERSION "2015-08-12"
#endif
#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
#endif
/// \file colvarmodule.h
/// \brief Collective variables main module
///
/// This file declares the main class for defining and manipulating
/// collective variables: there should be only one instance of this
/// class, because several variables are made static (i.e. they are
/// shared between all object instances) to be accessed from other
/// objects.
// Internal method return codes
#define COLVARS_NOT_IMPLEMENTED -2
#define COLVARS_ERROR -1
#define COLVARS_OK 0
// On error, values of the colvars module error register
#define GENERAL_ERROR 1
#define FILE_ERROR (1<<1)
#define MEMORY_ERROR (1<<2)
#define BUG_ERROR (1<<3) // Inconsistent state indicating bug
#define INPUT_ERROR (1<<4) // out of bounds or inconsistent input
#define DELETE_COLVARS (1<<5) // Instruct the caller to delete cvm
#define FATAL_ERROR (1<<6) // Should be set, or not, together with other bits
#include <iostream>
#include <iomanip>
#include <string>
#include <cstring>
#include <sstream>
#include <fstream>
#include <cmath>
#include <vector>
#include <list>
#ifdef NAMD_VERSION
// use Lustre-friendly wrapper to POSIX write()
#include "fstream_namd.h"
#endif
class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
class colvarscript;
/// \brief Collective variables module (main class)
///
/// Class to control the collective variables calculation. An object
/// (usually one) of this class is spawned from the MD program,
/// containing all i/o routines and general interface.
///
/// At initialization, the colvarmodule object creates a proxy object
/// to provide a transparent interface between the MD program and the
/// child objects
class colvarmodule {
private:
/// Impossible to initialize the main object without arguments
colvarmodule();
public:
friend class colvarproxy;
// TODO colvarscript should be unaware of colvarmodule's internals
friend class colvarscript;
/// Defining an abstract real number allows to switch precision
typedef double real;
/// Residue identifier
typedef int residue_id;
class rvector;
template <class T> class vector1d;
template <class T> class matrix2d;
class quaternion;
class rotation;
/// \brief Atom position (different type name from rvector, to make
/// possible future PBC-transparent implementations)
typedef rvector atom_pos;
/// \brief 3x3 matrix of real numbers
class rmatrix;
// allow these classes to access protected data
class atom;
class atom_group;
friend class atom;
friend class atom_group;
typedef std::vector<atom>::iterator atom_iter;
typedef std::vector<atom>::const_iterator atom_const_iter;
/// Module-wide error state
/// see constants at the top of this file
static int errorCode;
static inline void set_error_bits(int code)
{
errorCode |= code;
}
static inline int get_error()
{
return errorCode;
}
static inline void clear_error()
{
errorCode = 0;
}
/// Current step number
static long it;
/// Starting step number for this run
static long it_restart;
/// Return the current step number from the beginning of this run
static inline long step_relative()
{
return it - it_restart;
}
/// Return the current step number from the beginning of the whole
/// calculation
static inline long step_absolute()
{
return it;
}
/// If true, get it_restart from the state file; if set to false,
/// the MD program is providing it
bool it_restart_from_state_file;
/// \brief Finite difference step size (if there is no dynamics, or
/// if gradients need to be tested independently from the size of
/// dt)
static real debug_gradients_step_size;
/// Prefix for all output files for this run
static std::string output_prefix;
/// input restart file name (determined from input_prefix)
static std::string restart_in_name;
/// Array of collective variables
static std::vector<colvar *> colvars;
/* TODO: implement named CVCs
/// Array of named (reusable) collective variable components
static std::vector<cvc *> cvcs;
/// Named cvcs register themselves at initialization time
inline void register_cvc(cvc *p) {
cvcs.push_back(p);
}
*/
/// Array of collective variable biases
static std::vector<colvarbias *> biases;
/// \brief Number of ABF biases initialized (in normal conditions
/// should be 1)
static size_t n_abf_biases;
/// \brief Number of metadynamics biases initialized (in normal
/// conditions should be 1)
static size_t n_meta_biases;
/// \brief Number of restraint biases initialized (no limit on the
/// number)
static size_t n_rest_biases;
/// \brief Number of histograms initialized (no limit on the
/// number)
static size_t n_histo_biases;
/// \brief Whether debug output should be enabled (compile-time option)
static inline bool debug()
{
return COLVARS_DEBUG;
}
/// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name
colvarmodule(colvarproxy *proxy);
/// Destructor
~colvarmodule();
/// Actual function called by the destructor
int reset();
/// Open a config file, load its contents, and pass it to config_string()
int read_config_file(char const *config_file_name);
/// \brief Parse a config string assuming it is a complete configuration
/// (i.e. calling all parse functions)
int read_config_string(std::string const &conf);
/// \brief Parse a "clean" config string (no comments)
int parse_config(std::string &conf);
// Parse functions (setup internal data based on a string)
/// Parse the few module's global parameters
int parse_global_params(std::string const &conf);
/// Parse and initialize collective variables
int parse_colvars(std::string const &conf);
/// Parse and initialize collective variable biases
int parse_biases(std::string const &conf);
/// Parse and initialize collective variable biases of a specific type
template <class bias_type>
int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count);
/// Test error condition and keyword parsing
/// on error, delete new bias
bool check_new_bias(std::string &conf, char const *key);
// "Setup" functions (change internal data based on related data
// from the proxy that may change during program execution)
// No additional parsing is done within these functions
/// (Re)initialize internal data (currently used by LAMMPS)
/// Also calls setup() member functions of colvars and biases
int setup();
/// (Re)initialize and (re)read the input state file calling read_restart()
int setup_input();
/// (Re)initialize the output trajectory and state file (does not write it yet)
int setup_output();
#ifdef NAMD_VERSION
typedef ofstream_namd ofstream;
#else
typedef std::ofstream ofstream;
#endif
/// Read the input restart file
std::istream & read_restart(std::istream &is);
/// Write the output restart file
std::ostream & write_restart(std::ostream &os);
/// Open a trajectory file if requested (and leave it open)
int open_traj_file(std::string const &file_name);
/// Close it
int close_traj_file();
/// Write in the trajectory file
std::ostream & write_traj(std::ostream &os);
/// Write explanatory labels in the trajectory file
std::ostream & write_traj_label(std::ostream &os);
/// Write all FINAL output files
int write_output_files();
/// Backup a file before writing it
static int backup_file(char const *filename);
/// Look up a bias by name; returns NULL if not found
static colvarbias * bias_by_name(std::string const &name);
/// Look up a colvar by name; returns NULL if not found
static colvar * colvar_by_name(std::string const &name);
/// Load new configuration for the given bias -
/// currently works for harmonic (force constant and/or centers)
int change_configuration(std::string const &bias_name, std::string const &conf);
/// Read a colvar value
std::string read_colvar(std::string const &name);
/// Calculate change in energy from using alt. config. for the given bias -
/// currently works for harmonic (force constant and/or centers)
real energy_difference(std::string const &bias_name, std::string const &conf);
/// Give the total number of bins for a given bias.
int bias_bin_num(std::string const &bias_name);
/// Calculate the bin index for a given bias.
int bias_current_bin(std::string const &bias_name);
//// Give the count at a given bin index.
int bias_bin_count(std::string const &bias_name, size_t bin_index);
//// Share among replicas.
int bias_share(std::string const &bias_name);
/// Calculate collective variables and biases
int calc();
/// Perform analysis
int analyze();
/// \brief Read a collective variable trajectory (post-processing
/// only, not called at runtime)
int read_traj(char const *traj_filename,
long traj_read_begin,
long traj_read_end);
/// Quick conversion of an object to a string
template<typename T> static std::string to_str(T const &x,
size_t const &width = 0,
size_t const &prec = 0);
/// Quick conversion of a vector of objects to a string
template<typename T> static std::string to_str(std::vector<T> const &x,
size_t const &width = 0,
size_t const &prec = 0);
/// Reduce the number of characters in a string
static inline std::string wrap_string(std::string const &s,
size_t const &nchars)
{
if (!s.size())
return std::string(nchars, ' ');
else
return ( (s.size() <= size_t(nchars)) ?
(s+std::string(nchars-s.size(), ' ')) :
(std::string(s, 0, nchars)) );
}
/// Number of characters to represent a time step
static size_t const it_width;
/// Number of digits to represent a collective variables value(s)
static size_t const cv_prec;
/// Number of characters to represent a collective variables value(s)
static size_t const cv_width;
/// Number of digits to represent the collective variables energy
static size_t const en_prec;
/// Number of characters to represent the collective variables energy
static size_t const en_width;
/// Line separator in the log output
static std::string const line_marker;
// proxy functions
/// \brief Value of the unit for atomic coordinates with respect to
/// angstroms (used by some variables for hard-coded default values)
static real unit_angstrom();
/// \brief Boltmann constant
static real boltzmann();
/// \brief Temperature of the simulation (K)
static real temperature();
/// \brief Time step of MD integrator (fs)
static real dt();
/// Request calculation of system force from MD engine
static void request_system_force();
/// Print a message to the main log
static void log(std::string const &message);
/// Print a message to the main log and exit with error code
static void fatal_error(std::string const &message);
/// Print a message to the main log and set global error code
static void error(std::string const &message, int code = GENERAL_ERROR);
/// Print a message to the main log and exit normally
static void exit(std::string const &message);
// Replica exchange commands.
static bool replica_enabled();
static int replica_index();
static int replica_num();
static void replica_comm_barrier();
static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
/// \brief Get the distance between two atomic positions with pbcs handled
/// correctly
static rvector position_distance(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the square distance between two positions (with
/// periodic boundary conditions handled transparently)
///
/// Note: in the case of periodic boundary conditions, this provides
/// an analytical square distance (while taking the square of
/// position_distance() would produce leads to a cusp)
static real position_dist2(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the closest periodic image to a reference position
/// \param pos The position to look for the closest periodic image
/// \param ref_pos (optional) The reference position
static void select_closest_image(atom_pos &pos,
atom_pos const &ref_pos);
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
static void select_closest_images(std::vector<atom_pos> &pos,
atom_pos const &ref_pos);
/// \brief Names of groups from a Gromacs .ndx file to be read at startup
static std::list<std::string> index_group_names;
/// \brief Groups from a Gromacs .ndx file read at startup
static std::list<std::vector<int> > index_groups;
/// \brief Read a Gromacs .ndx file
static int read_index_file(char const *filename);
/// \brief Create atoms from a file \param filename name of the file
/// (usually a PDB) \param atoms array of the atoms to be allocated
/// \param pdb_field (optiona) if "filename" is a PDB file, use this
/// field to determine which are the atoms to be set
static int load_atoms(char const *filename,
std::vector<atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (PDB or XYZ)
static int load_coords(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from an
/// XYZ file
static int load_coords_xyz(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices);
/// Frequency for collective variables trajectory output
static size_t cv_traj_freq;
/// \brief True if only analysis is performed and not a run
static bool b_analysis;
/// Frequency for saving output restarts
static size_t restart_out_freq;
/// Output restart file name
std::string restart_out_name;
/// Pseudo-random number with Gaussian distribution
static real rand_gaussian(void);
protected:
/// Configuration file
std::ifstream config_s;
/// Configuration file parser object
colvarparse *parse;
/// Name of the trajectory file
std::string cv_traj_name;
/// Collective variables output trajectory file
colvarmodule::ofstream cv_traj_os;
/// Appending to the existing trajectory file?
bool cv_traj_append;
/// Output restart file
colvarmodule::ofstream restart_out_os;
/// \brief Counter for the current depth in the object hierarchy (useg e.g. in output
static size_t depth;
/// Use scripted colvars forces?
static bool use_scripted_forces;
public:
/// \brief Pointer to the proxy object, used to retrieve atomic data
/// from the hosting program; it is static in order to be accessible
/// from static functions in the colvarmodule class
static colvarproxy *proxy;
/// Increase the depth (number of indentations in the output)
static void increase_depth();
/// Decrease the depth (number of indentations in the output)
static void decrease_depth();
static inline bool scripted_forces() { return use_scripted_forces; }
};
/// Shorthand for the frequently used type prefix
typedef colvarmodule cvm;
#include "colvartypes.h"
std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);
template<typename T> std::string cvm::to_str(T const &x,
size_t const &width,
size_t const &prec) {
std::ostringstream os;
if (width) os.width(width);
if (prec) {
os.setf(std::ios::scientific, std::ios::floatfield);
os.precision(prec);
}
os << x;
return os.str();
}
template<typename T> std::string cvm::to_str(std::vector<T> const &x,
size_t const &width,
size_t const &prec) {
if (!x.size()) return std::string("");
std::ostringstream os;
if (prec) {
os.setf(std::ios::scientific, std::ios::floatfield);
}
os << "{ ";
if (width) os.width(width);
if (prec) os.precision(prec);
os << x[0];
for (size_t i = 1; i < x.size(); i++) {
os << ", ";
if (width) os.width(width);
if (prec) os.precision(prec);
os << x[i];
}
os << " }";
return os.str();
}
#include "colvarproxy.h"
inline cvm::real cvm::unit_angstrom()
{
return proxy->unit_angstrom();
}
inline cvm::real cvm::boltzmann()
{
return proxy->boltzmann();
}
inline cvm::real cvm::temperature()
{
return proxy->temperature();
}
inline cvm::real cvm::dt()
{
return proxy->dt();
}
// Replica exchange commands
inline bool cvm::replica_enabled() {
return proxy->replica_enabled();
}
inline int cvm::replica_index() {
return proxy->replica_index();
}
inline int cvm::replica_num() {
return proxy->replica_num();
}
inline void cvm::replica_comm_barrier() {
return proxy->replica_comm_barrier();
}
inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
return proxy->replica_comm_recv(msg_data,buf_len,src_rep);
}
inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
return proxy->replica_comm_send(msg_data,msg_len,dest_rep);
}
inline void cvm::request_system_force()
{
proxy->request_system_force(true);
}
inline void cvm::select_closest_image(atom_pos &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_image(pos, ref_pos);
}
inline void cvm::select_closest_images(std::vector<atom_pos> &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_images(pos, ref_pos);
}
inline cvm::rvector cvm::position_distance(atom_pos const &pos1,
atom_pos const &pos2)
{
return proxy->position_distance(pos1, pos2);
}
inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return proxy->position_dist2(pos1, pos2);
}
inline cvm::real cvm::rand_gaussian(void)
{
return proxy->rand_gaussian();
}
#endif