diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 9788f3df9..fa38c4d33 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1,2032 +1,2035 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvarscript.h"
#include <algorithm>
/// Compare two cvcs using their names
/// Used to sort CVC array in scripted coordinates
bool compare(colvar::cvc *i, colvar::cvc *j) {
return i->name < j->name;
}
colvar::colvar(std::string const &conf)
: colvarparse(conf)
{
size_t i;
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()){ \
std::ostringstream s; \
s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;\
cvcs.back()->name = s.str(); \
/* pad cvc number for correct ordering when sorting by name */\
} \
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("dipole angle", "dipoleAngle", dipole_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 their names
// Note: default CVC names are in input order for same type of CVC
std::sort(cvcs.begin(), cvcs.end(), compare);
if(cvcs.size() > 1) {
cvm::log("Sorted list of components for this scripted colvar:");
for (i = 0; i < cvcs.size(); i++) {
cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name);
}
}
// Build ordered list of component values that will be
// passed to the script
for (i = 0; i < cvcs.size(); i++) {
sorted_cvc_values.push_back(&(cvcs[i]->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++) {
cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
temp_id_list.push_back(ag[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);
tasks[t] = true;
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");
tasks[t] = true;
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());
tasks[t] = true;
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());
tasks[t] = true;
break;
case task_output_applied_force:
case task_lower_wall:
case task_upper_wall:
// all of the above require gradients
enable(task_gradients);
tasks[t] = true;
break;
case task_fdiff_velocity:
x_old.type(value());
v_fdiff.type(value());
v_reported.type(value());
tasks[t] = true;
break;
case task_output_velocity:
enable(task_fdiff_velocity);
tasks[t] = true;
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);
}
tasks[t] = true;
break;
case task_extended_lagrangian:
enable(task_gradients);
v_reported.type(value());
tasks[t] = true;
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);
}
tasks[t] = true;
break;
case task_output_value:
case task_runave:
case task_corrfunc:
case task_langevin:
case task_output_energy:
case task_scripted:
case task_gradients:
tasks[t] = true;
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();
}
tasks[t] = true;
break;
default:
break;
}
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);
tasks[t] = false;
break;
case task_system_force:
disable(task_output_system_force);
tasks[t] = false;
break;
case task_Jacobian_force:
disable(task_report_Jacobian_force);
tasks[t] = false;
break;
case task_fdiff_velocity:
disable(task_output_velocity);
tasks[t] = false;
break;
case task_lower_boundary:
case task_upper_boundary:
disable(task_grid);
tasks[t] = false;
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_langevin:
case task_output_energy:
case task_collect_gradients:
case task_scripted:
tasks[t] = false;
break;
default:
break;
}
}
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.setup();
atoms.reset_mass(name,i,ig);
atoms.read_positions();
}
}
}
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()
{
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\".\n");
update_cvc_flags();
calc_cvcs();
calc_colvar_properties();
}
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
size_t i, ig;
size_t cvc_count;
size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
// prepare atom groups for calculation
if (cvm::debug())
cvm::log("Collecting data from atom groups.\n");
if (first_cvc >= cvcs.size()) {
cvm::error("Error: trying to address a component outside the "
"range defined for colvar \""+name+"\".\n", BUG_ERROR);
return BUG_ERROR;
}
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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();
atoms.calc_required_properties();
// 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 = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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, calculate component values
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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, using either a scripted function or a polynomial
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 COLVARS_NOT_IMPLEMENTED;
}
if (res != COLVARS_OK) {
cvm::error("Error running scripted colvar");
return COLVARS_OK;
}
} else if (x.type() == colvarvalue::type_scalar) {
// polynomial combination allowed
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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 = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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");
// calculate the gradients of each component
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
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.", COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
// Collect the atomic gradients inside colvar object
for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
atomic_gradients[a].reset();
}
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
// 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++) {
cvm::atom_group &ag = *(cvcs[i]->atom_groups[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++) {
size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
ag[k].id) - atom_ids.begin();
atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
}
} else {
for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
ag[k].id) - atom_ids.begin();
atomic_gradients[a] += coeff * ag[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.", COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
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 = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->b_enabled) continue;
cvc_count++;
(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");
}
+ return COLVARS_OK;
}
int colvar::calc_colvar_properties()
{
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");
+
+ return COLVARS_OK;
}
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_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: bias force on extended coordinate (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
// (note: wall potential is added to f after this block)
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);
}
// Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use
if (tasks[task_lower_wall] || tasks[task_upper_wall]) {
// Wall force
colvarvalue fw(x);
fw.reset();
if (cvm::debug())
cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
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;
f += fw;
if (cvm::debug())
cvm::log("Applying a lower wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
}
} else {
cvm::real const grad = this->dist2_lgrad(x, upper_wall);
if (grad > 0.0) {
fw = -0.5 * upper_wall_k * grad;
f += fw;
if (cvm::debug())
cvm::log("Applying an upper wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
}
}
}
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;
func_grads.reserve(cvcs.size());
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.", COLVARS_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[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)
{
if (flags.size() != cvcs.size()) {
cvm::error("ERROR: Wrong number of CVC flags provided.");
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;
}
int colvar::update_cvc_flags()
{
size_t i;
// 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 COLVARS_ERROR;
}
cvc_flags.resize(0);
}
return COLVARS_OK;
}
int colvar::num_active_cvcs() const
{
int result = 0;
for (size_t i = 0; i < cvcs.size(); i++) {
if (cvcs[i]->b_enabled) result++;
}
return result;
}
// ******************** 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/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 59ac775cb..f7c33c5e9 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -1,1253 +1,1251 @@
/// -*- c++ -*-
#include <sstream>
#include <string.h>
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias.h"
#include "colvarbias_abf.h"
#include "colvarbias_alb.h"
#include "colvarbias_histogram.h"
#include "colvarbias_meta.h"
#include "colvarbias_restraint.h"
#include "colvarscript.h"
colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
// pointer to the proxy object
if (proxy == NULL) {
proxy = proxy_in;
parse = new colvarparse();
} else {
// TODO relax this error to handle multiple molecules in VMD
// once the module is not static anymore
cvm::error("Error: trying to allocate the collective "
"variable module twice.\n");
return;
}
cvm::log(cvm::line_marker);
cvm::log("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n");
// set initial default values
// "it_restart" will be set by the input state file, if any;
// "it" should be updated by the proxy
colvarmodule::it = colvarmodule::it_restart = 0;
colvarmodule::it_restart_from_state_file = true;
colvarmodule::use_scripted_forces = false;
colvarmodule::b_analysis = false;
colvarmodule::debug_gradients_step_size = 1.0e-07;
colvarmodule::rotation::monitor_crossings = false;
colvarmodule::rotation::crossing_threshold = 1.0e-02;
colvarmodule::cv_traj_freq = 100;
colvarmodule::restart_out_freq = proxy->restart_frequency();
// by default overwrite the existing trajectory file
colvarmodule::cv_traj_append = false;
}
int colvarmodule::read_config_file(char const *config_filename)
{
cvm::log(cvm::line_marker);
cvm::log("Reading new configuration from file \""+
std::string(config_filename)+"\":\n");
// open the configfile
config_s.open(config_filename);
if (!config_s.is_open()) {
cvm::error("Error: in opening configuration file \""+
std::string(config_filename)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
}
// read the config file into a string
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
conf.append(line+"\n");
}
config_s.close();
return parse_config(conf);
}
int colvarmodule::read_config_string(std::string const &config_str)
{
cvm::log(cvm::line_marker);
cvm::log("Reading new configuration:\n");
std::istringstream config_s(config_str);
// strip the comments away
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
conf.append(line+"\n");
}
return parse_config(conf);
}
int colvarmodule::parse_config(std::string &conf)
{
- int error_code = COLVARS_OK;
-
// parse global options
if (catch_input_errors(parse_global_params(conf))) {
return get_error();
}
// parse the options for collective variables
if (catch_input_errors(parse_colvars(conf))) {
return get_error();
}
// parse the options for biases
if (catch_input_errors(parse_biases(conf))) {
return get_error();
}
// done parsing known keywords, check that all keywords found were valid ones
if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
return get_error();
}
cvm::log(cvm::line_marker);
cvm::log("Collective variables module (re)initialized.\n");
cvm::log(cvm::line_marker);
if (cv_traj_os.is_open()) {
// configuration might have changed, better redo the labels
write_traj_label(cv_traj_os);
}
return get_error();
}
int colvarmodule::parse_global_params(std::string const &conf)
{
std::string index_file_name;
if (parse->get_keyval(conf, "indexFile", index_file_name)) {
read_index_file(index_file_name.c_str());
}
parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
debug_gradients_step_size,
colvarparse::parse_silent);
parse->get_keyval(conf, "monitorEigenvalueCrossing",
colvarmodule::rotation::monitor_crossings,
colvarmodule::rotation::monitor_crossings,
colvarparse::parse_silent);
parse->get_keyval(conf, "eigenvalueCrossingThreshold",
colvarmodule::rotation::crossing_threshold,
colvarmodule::rotation::crossing_threshold,
colvarparse::parse_silent);
parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq);
parse->get_keyval(conf, "colvarsRestartFrequency",
restart_out_freq, restart_out_freq);
// if this is true when initializing, it means
// we are continuing after a reset(): default to true
parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append);
parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false,
colvarparse::parse_silent);
if (use_scripted_forces && !proxy->force_script_defined) {
cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::parse_colvars(std::string const &conf)
{
if (cvm::debug())
cvm::log("Initializing the collective variables.\n");
std::string colvar_conf = "";
size_t pos = 0;
while (parse->key_lookup(conf, "colvar", colvar_conf, pos)) {
if (colvar_conf.size()) {
cvm::log(cvm::line_marker);
cvm::increase_depth();
colvars.push_back(new colvar(colvar_conf));
if (cvm::get_error() ||
((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
cvm::log("Error while constructing colvar number " +
cvm::to_str(colvars.size()) + " : deleting.");
delete colvars.back(); // the colvar destructor updates the colvars array
return COLVARS_ERROR;
}
cvm::decrease_depth();
} else {
cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
cvm::decrease_depth();
colvar_conf = "";
}
if (!colvars.size()) {
cvm::log("Warning: no collective variables defined.\n");
}
if (colvars.size())
cvm::log(cvm::line_marker);
cvm::log("Collective variables initialized, "+
cvm::to_str(colvars.size())+
" in total.\n");
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
bool colvarmodule::check_new_bias(std::string &conf, char const *key)
{
if (cvm::get_error() ||
(biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
cvm::log("Error while constructing bias number " +
cvm::to_str(biases.size()) + " : deleting.\n");
delete biases.back(); // the bias destructor updates the biases array
return true;
}
return false;
}
template <class bias_type>
int colvarmodule::parse_biases_type(std::string const &conf,
char const *keyword,
size_t &bias_count)
{
std::string bias_conf = "";
size_t conf_saved_pos = 0;
while (parse->key_lookup(conf, keyword, bias_conf, conf_saved_pos)) {
if (bias_conf.size()) {
cvm::log(cvm::line_marker);
cvm::increase_depth();
biases.push_back(new bias_type(bias_conf, keyword));
if (cvm::check_new_bias(bias_conf, keyword)) {
return COLVARS_ERROR;
}
cvm::decrease_depth();
bias_count++;
} else {
cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n",
INPUT_ERROR);
return COLVARS_ERROR;
}
bias_conf = "";
}
return COLVARS_OK;
}
int colvarmodule::parse_biases(std::string const &conf)
{
if (cvm::debug())
cvm::log("Initializing the collective variables biases.\n");
/// initialize ABF instances
parse_biases_type<colvarbias_abf>(conf, "abf", n_abf_biases);
/// initialize adaptive linear biases
parse_biases_type<colvarbias_alb>(conf, "ALB", n_rest_biases);
/// initialize harmonic restraints
parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);
/// initialize histograms
parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
/// initialize linear restraints
parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);
/// initialize metadynamics instances
parse_biases_type<colvarbias_meta>(conf, "metadynamics", n_meta_biases);
if (use_scripted_forces) {
cvm::log(cvm::line_marker);
cvm::increase_depth();
cvm::log("User forces script will be run at each bias update.");
cvm::decrease_depth();
}
if (biases.size() || use_scripted_forces) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables biases initialized, "+
cvm::to_str(biases.size())+" in total.\n");
} else {
if (!use_scripted_forces) {
cvm::log("No collective variables biases were defined.\n");
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::catch_input_errors(int result)
{
if (result != COLVARS_OK || get_error()) {
set_error_bits(result | INPUT_ERROR);
parse->reset();
return get_error();
}
return COLVARS_OK;
}
colvarbias * colvarmodule::bias_by_name(std::string const &name) {
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ((*bi)->name == name) {
return (*bi);
}
}
return NULL;
}
colvar *colvarmodule::colvar_by_name(std::string const &name) {
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if ((*cvi)->name == name) {
return (*cvi);
}
}
return NULL;
}
int colvarmodule::change_configuration(std::string const &bias_name,
std::string const &conf)
{
// This is deprecated; supported strategy is to delete the bias
// and parse the new config
cvm::increase_depth();
colvarbias *b;
b = bias_by_name(bias_name);
if (b == NULL) { cvm::error("Error: bias not found: " + bias_name); }
b->change_configuration(conf);
cvm::decrease_depth();
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
std::string colvarmodule::read_colvar(std::string const &name)
{
cvm::increase_depth();
colvar *c;
std::stringstream ss;
c = colvar_by_name(name);
if (c == NULL) { cvm::fatal_error("Error: colvar not found: " + name); }
ss << c->value();
cvm::decrease_depth();
return ss.str();
}
cvm::real colvarmodule::energy_difference(std::string const &bias_name,
std::string const &conf)
{
cvm::increase_depth();
colvarbias *b;
cvm::real energy_diff = 0.;
b = bias_by_name(bias_name);
if (b == NULL) { cvm::fatal_error("Error: bias not found: " + bias_name); }
energy_diff = b->energy_difference(conf);
cvm::decrease_depth();
return energy_diff;
}
int colvarmodule::bias_current_bin(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->current_bin();
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_bin_num(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->bin_num();
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->bin_count(bin_index);
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_share(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
b->replica_share();
ret = COLVARS_OK;
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::calc() {
cvm::real total_bias_energy = 0.0;
cvm::real total_colvar_energy = 0.0;
std::vector<colvar *>::iterator cvi;
std::vector<colvarbias *>::iterator bi;
if (cvm::debug()) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables module, step no. "+
cvm::to_str(cvm::step_absolute())+"\n");
}
// calculate collective variables and their gradients
if (cvm::debug())
cvm::log("Calculating collective variables.\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->calc();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
// update the biases and communicate their forces to the collective
// variables
if (cvm::debug() && biases.size())
cvm::log("Updating collective variable biases.\n");
cvm::increase_depth();
for (bi = biases.begin(); bi != biases.end(); bi++) {
total_bias_energy += (*bi)->update();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
// sum the forces from all biases for each collective variable
if (cvm::debug() && biases.size())
cvm::log("Collecting forces from all biases.\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->reset_bias_force();
}
for (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
// Run user force script, if provided,
// potentially adding scripted forces to the colvars
if (use_scripted_forces) {
int res;
res = proxy->run_force_callback();
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Colvar forces scripts are not implemented.");
return COLVARS_ERROR;
}
if (res != COLVARS_OK) {
cvm::error("Error running user colvar forces script");
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
if (cvm::b_analysis) {
// perform runtime analysis of colvars and biases
if (cvm::debug() && biases.size())
cvm::log("Perform runtime analyses.\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->analyze();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
for (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->analyze();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
}
// sum up the forces for each colvar, including wall forces
// and integrate any internal
// equation of motion (extended system)
if (cvm::debug())
cvm::log("Updating the internal degrees of freedom "
"of colvars (if they have any).\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
total_colvar_energy += (*cvi)->update();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
proxy->add_energy(total_bias_energy + total_colvar_energy);
// make collective variables communicate their forces to their
// coupled degrees of freedom (i.e. atoms)
if (cvm::debug())
cvm::log("Communicating forces from the colvars to the atoms.\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
if ((*cvi)->tasks[colvar::task_gradients]) {
(*cvi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
}
cvm::decrease_depth();
// write restart file, if needed
if (restart_out_freq && restart_out_name.size()) {
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Writing the state file \""+
restart_out_name+"\".\n");
proxy->backup_file(restart_out_name.c_str());
restart_out_os.open(restart_out_name.c_str());
if (!restart_out_os.is_open() || !write_restart(restart_out_os))
cvm::error("Error: in writing restart file.\n");
restart_out_os.close();
}
}
// write trajectory file, if needed
if (cv_traj_freq && cv_traj_name.size()) {
if (!cv_traj_os.is_open()) {
open_traj_file(cv_traj_name);
}
// write labels in the traj file every 1000 lines and at first timestep
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
write_traj_label(cv_traj_os);
}
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
write_traj(cv_traj_os);
}
if (restart_out_freq && cv_traj_os.is_open()) {
// flush the trajectory file if we are at the restart frequency
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+
cv_traj_name+"\".\n");
cv_traj_os.flush();
}
}
} // end if (cv_traj_freq)
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::analyze()
{
if (cvm::debug()) {
cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n");
}
if (cvm::step_relative() == 0)
cvm::log("Performing analysis.\n");
// perform colvar-specific analysis
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
cvm::increase_depth();
(*cvi)->analyze();
cvm::decrease_depth();
}
// perform bias-specific analysis
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
cvm::increase_depth();
(*bi)->analyze();
cvm::decrease_depth();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::setup()
{
// loop over all components of all colvars to reset masses of all groups
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end(); cvi++) {
(*cvi)->setup();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
colvarmodule::~colvarmodule()
{
reset();
delete parse;
proxy = NULL;
}
int colvarmodule::reset()
{
parse->reset();
cvm::log("Resetting the Collective Variables Module.\n");
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
cvi != colvars.rend();
cvi++) {
delete *cvi; // the colvar destructor updates the colvars array
}
colvars.clear();
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
bi != biases.rend();
bi++) {
delete *bi; // the bias destructor updates the biases array
}
biases.clear();
index_groups.clear();
index_group_names.clear();
if (cv_traj_os.is_open()) {
// Do not close file here, as we might not be done with it yet.
cv_traj_os.flush();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::setup_input()
{
// name of input state file
restart_in_name = proxy->input_prefix().size() ?
std::string(proxy->input_prefix()+".colvars.state") :
std::string("") ;
// read the restart configuration, if available
if (restart_in_name.size()) {
// read the restart file
std::ifstream input_is(restart_in_name.c_str());
if (!input_is.good()) {
cvm::error("Error: in opening restart file \""+
std::string(restart_in_name)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
} else {
cvm::log("Restarting from file \""+restart_in_name+"\".\n");
read_restart(input_is);
if (cvm::get_error() != COLVARS_OK) {
return COLVARS_ERROR;
}
cvm::log(cvm::line_marker);
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::setup_output()
{
int error_code = 0;
// output state file (restart)
restart_out_name = proxy->restart_output_prefix().size() ?
std::string(proxy->restart_output_prefix()+".colvars.state") :
std::string("");
if (restart_out_name.size()) {
cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
}
output_prefix = proxy->output_prefix();
if (output_prefix.size()) {
cvm::log("The final output state file will be \""+
(output_prefix.size() ?
std::string(output_prefix+".colvars.state") :
std::string("colvars.state"))+"\".\n");
// cvm::log (cvm::line_marker);
}
cv_traj_name =
(output_prefix.size() ?
std::string(output_prefix+".colvars.traj") :
std::string(""));
if (cv_traj_freq && cv_traj_name.size()) {
error_code |= open_traj_file(cv_traj_name);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
error_code |= (*bi)->setup_output();
}
if (error_code != COLVARS_OK || cvm::get_error()) {
set_error_bits(FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
std::istream & colvarmodule::read_restart(std::istream &is)
{
{
// read global restart information
std::string restart_conf;
if (is >> colvarparse::read_block("configuration", restart_conf)) {
if (it_restart_from_state_file) {
parse->get_keyval(restart_conf, "step",
it_restart, (size_t) 0,
colvarparse::parse_silent);
it = it_restart;
}
}
is.clear();
parse->clear_keyword_registry();
}
// colvars restart
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if ( !((*cvi)->read_restart(is)) ) {
cvm::error("Error: in reading restart configuration for collective variable \""+
(*cvi)->name+"\".\n",
INPUT_ERROR);
}
}
// biases restart
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if (!((*bi)->read_restart(is))) {
cvm::error("Error: in reading restart configuration for bias \""+
(*bi)->name+"\".\n",
INPUT_ERROR);
}
}
cvm::decrease_depth();
return is;
}
int colvarmodule::backup_file(char const *filename)
{
return proxy->backup_file(filename);
}
int colvarmodule::write_output_files()
{
// if this is a simulation run (i.e. not a postprocessing), output data
// must be written to be able to restart the simulation
std::string const out_name =
(output_prefix.size() ?
std::string(output_prefix+".colvars.state") :
std::string("colvars.state"));
cvm::log("Saving collective variables state to \""+out_name+"\".\n");
std::ostream * os = proxy->output_stream(out_name);
os->setf(std::ios::scientific, std::ios::floatfield);
this->write_restart(*os);
proxy->close_output_stream(out_name);
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_output_files();
}
cvm::decrease_depth();
cvm::increase_depth();
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_output_files();
}
cvm::decrease_depth();
if (cv_traj_os.is_open()) {
// do not close to avoid problems with multiple NAMD runs
cv_traj_os.flush();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::read_traj(char const *traj_filename,
long traj_read_begin,
long traj_read_end)
{
cvm::log("Opening trajectory file \""+
std::string(traj_filename)+"\".\n");
std::ifstream traj_is(traj_filename);
while (true) {
while (true) {
std::string line("");
do {
if (!colvarparse::getline_nocomments(traj_is, line)) {
cvm::log("End of file \""+std::string(traj_filename)+
"\" reached, or corrupted file.\n");
traj_is.close();
return false;
}
} while (line.find_first_not_of(colvarparse::white_space) == std::string::npos);
std::istringstream is(line);
if (!(is >> it)) return false;
if ( (it < traj_read_begin) ) {
if ((it % 1000) == 0)
std::cerr << "Skipping trajectory step " << it
<< " \r";
continue;
} else {
if ((it % 1000) == 0)
std::cerr << "Reading from trajectory, step = " << it
<< " \r";
if ( (traj_read_end > traj_read_begin) &&
(it > traj_read_end) ) {
std::cerr << "\n";
cvm::error("Reached the end of the trajectory, "
"read_end = "+cvm::to_str(traj_read_end)+"\n",
FILE_ERROR);
return COLVARS_ERROR;
}
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if (!(*cvi)->read_traj(is)) {
cvm::error("Error: in reading colvar \""+(*cvi)->name+
"\" from trajectory file \""+
std::string(traj_filename)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
}
}
break;
}
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
std::ostream & colvarmodule::write_restart(std::ostream &os)
{
os.setf(std::ios::scientific, std::ios::floatfield);
os << "configuration {\n"
<< " step " << std::setw(it_width)
<< it << "\n"
<< " dt " << dt() << "\n"
<< "}\n\n";
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_restart(os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_restart(os);
}
cvm::decrease_depth();
return os;
}
int colvarmodule::open_traj_file(std::string const &file_name)
{
if (cv_traj_os.is_open()) {
return COLVARS_OK;
}
// (re)open trajectory file
if (cv_traj_append) {
cvm::log("Appending to colvar trajectory file \""+file_name+
"\".\n");
cv_traj_os.open(file_name.c_str(), std::ios::app);
} else {
cvm::log("Writing to colvar trajectory file \""+file_name+
"\".\n");
proxy->backup_file(file_name.c_str());
cv_traj_os.open(file_name.c_str());
}
if (!cv_traj_os.is_open()) {
cvm::error("Error: cannot write to file \""+file_name+"\".\n",
FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::close_traj_file()
{
if (cv_traj_os.good()) {
cv_traj_os.close();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
std::ostream & colvarmodule::write_traj_label(std::ostream &os)
{
if (!os.good()) {
cvm::error("Cannot write to trajectory file.");
return os;
}
os.setf(std::ios::scientific, std::ios::floatfield);
os << "# " << cvm::wrap_string("step", cvm::it_width-2)
<< " ";
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_traj_label(os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_traj_label(os);
}
os << "\n";
if (cvm::debug()) {
os.flush();
}
cvm::decrease_depth();
return os;
}
std::ostream & colvarmodule::write_traj(std::ostream &os)
{
os.setf(std::ios::scientific, std::ios::floatfield);
os << std::setw(cvm::it_width) << it
<< " ";
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_traj(os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_traj(os);
}
os << "\n";
if (cvm::debug()) {
os.flush();
}
cvm::decrease_depth();
return os;
}
void cvm::log(std::string const &message)
{
if (depth > 0)
proxy->log((std::string(2*depth, ' '))+message);
else
proxy->log(message);
}
void cvm::increase_depth()
{
depth++;
}
void cvm::decrease_depth()
{
if (depth) depth--;
}
void cvm::error(std::string const &message, int code)
{
set_error_bits(code);
proxy->error(message);
}
void cvm::fatal_error(std::string const &message)
{
// TODO once all non-fatal errors have been set to be handled by error(),
// set DELETE_COLVARS here for VMD to handle it
set_error_bits(FATAL_ERROR);
proxy->fatal_error(message);
}
void cvm::exit(std::string const &message)
{
proxy->exit(message);
}
int cvm::read_index_file(char const *filename)
{
std::ifstream is(filename, std::ios::binary);
if (!is.good()) {
cvm::error("Error: in opening index file \""+
std::string(filename)+"\".\n",
FILE_ERROR);
}
while (is.good()) {
char open, close;
std::string group_name;
if ( (is >> open) && (open == '[') &&
(is >> group_name) &&
(is >> close) && (close == ']') ) {
for (std::list<std::string>::iterator names_i = index_group_names.begin();
names_i != index_group_names.end();
names_i++) {
if (*names_i == group_name) {
cvm::error("Error: the group name \""+group_name+
"\" appears in multiple index files.\n",
FILE_ERROR);
}
}
cvm::index_group_names.push_back(group_name);
cvm::index_groups.push_back(std::vector<int> ());
} else {
cvm::error("Error: in parsing index file \""+
std::string(filename)+"\".\n",
INPUT_ERROR);
}
int atom_number = 1;
size_t pos = is.tellg();
while ( (is >> atom_number) && (atom_number > 0) ) {
(cvm::index_groups.back()).push_back(atom_number);
pos = is.tellg();
}
is.clear();
is.seekg(pos, std::ios::beg);
std::string delim;
if ( (is >> delim) && (delim == "[") ) {
// new group
is.clear();
is.seekg(pos, std::ios::beg);
} else {
break;
}
}
cvm::log("The following index groups were read from the index file \""+
std::string(filename)+"\":\n");
std::list<std::string>::iterator names_i = index_group_names.begin();
std::list<std::vector<int> >::iterator lists_i = index_groups.begin();
for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) {
cvm::log(" "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n");
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int cvm::load_atoms(char const *file_name,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value)
{
return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
}
int cvm::load_coords(char const *file_name,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value)
{
// Differentiate between PDB and XYZ files
// for XYZ files, use CVM internal parser
// otherwise call proxy function for PDB
std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) {
if ( pdb_field.size() > 0 ) {
cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
return cvm::load_coords_xyz(file_name, pos, indices);
} else {
return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value);
}
}
int cvm::load_coords_xyz(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices)
{
std::ifstream xyz_is(filename);
unsigned int natoms;
char symbol[256];
std::string line;
if ( ! (xyz_is >> natoms) ) {
cvm::error("Error: cannot parse XYZ file "
+ std::string(filename) + ".\n", INPUT_ERROR);
}
// skip comment line
std::getline(xyz_is, line);
std::getline(xyz_is, line);
xyz_is.width(255);
std::vector<atom_pos>::iterator pos_i = pos.begin();
if (pos.size() != natoms) { // Use specified indices
int next = 0; // indices are zero-based
std::vector<int>::const_iterator index = indices.begin();
for ( ; pos_i != pos.end() ; pos_i++, index++) {
while ( next < *index ) {
std::getline(xyz_is, line);
next++;
}
xyz_is >> symbol;
xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
}
} else { // Use all positions
for ( ; pos_i != pos.end() ; pos_i++) {
xyz_is >> symbol;
xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
// static pointers
std::vector<colvar *> colvarmodule::colvars;
std::vector<colvarbias *> colvarmodule::biases;
size_t colvarmodule::n_abf_biases = 0;
size_t colvarmodule::n_rest_biases = 0;
size_t colvarmodule::n_histo_biases = 0;
size_t colvarmodule::n_meta_biases = 0;
colvarproxy *colvarmodule::proxy = NULL;
// static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03;
int colvarmodule::errorCode = 0;
long colvarmodule::it = 0;
long colvarmodule::it_restart = 0;
size_t colvarmodule::restart_out_freq = 0;
size_t colvarmodule::cv_traj_freq = 0;
size_t colvarmodule::depth = 0;
bool colvarmodule::b_analysis = false;
std::list<std::string> colvarmodule::index_group_names;
std::list<std::vector<int> > colvarmodule::index_groups;
bool colvarmodule::use_scripted_forces = false;
// file name prefixes
std::string colvarmodule::output_prefix = "";
std::string colvarmodule::restart_in_name = "";
// i/o constants
size_t const colvarmodule::it_width = 12;
size_t const colvarmodule::cv_prec = 14;
size_t const colvarmodule::cv_width = 21;
size_t const colvarmodule::en_prec = 14;
size_t const colvarmodule::en_width = 21;
std::string const colvarmodule::line_marker =
"----------------------------------------------------------------------\n";
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 337180d2a..ba3be8453 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -1,631 +1,631 @@
/// -*- c++ -*-
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-02-24"
+#define COLVARS_VERSION "2016-02-28"
#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.
// Error codes
#define COLVARS_OK 0
#define COLVARS_ERROR -1
#define GENERAL_ERROR -1 // TODO this can be simply merged with COLVARS_ERROR
#define COLVARS_NOT_IMPLEMENTED (-1<<1)
#define INPUT_ERROR (-1<<2) // out of bounds or inconsistent input
#define BUG_ERROR (-1<<3) // Inconsistent state indicating bug
#define FILE_ERROR (-1<<4)
#define MEMORY_ERROR (-1<<5)
#define FATAL_ERROR (-1<<6) // Should be set, or not, together with other bits
#define DELETE_COLVARS (-1<<7) // Instruct the caller to delete cvm
#define COLVARS_NO_SUCH_FRAME (-1<<8) // Cannot load the requested frame
#include <iostream>
#include <iomanip>
#include <string>
#include <cstring>
#include <sstream>
#include <fstream>
#include <cmath>
#include <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;
errorCode |= COLVARS_ERROR;
}
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);
private:
/// Useful wrapper to interrupt parsing if any error occurs
int catch_input_errors(int result);
public:
// "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,
atom_group &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