git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14830 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-04-15 16:07:28 +00:00
parent 212a955285
commit f50b03fcab
2 changed files with 764 additions and 0 deletions

499
lib/colvars/colvardeps.cpp Normal file
View File

@ -0,0 +1,499 @@
#include "colvarmodule.h"
#include "colvardeps.h"
cvm::deps::~deps() {
size_t i;
for (i=0; i<feature_states.size(); i++) {
if (feature_states[i] != NULL) delete feature_states[i];
}
// Do not delete features if it's static
// for (i=0; i<features.size(); i++) {
// if (features[i] != NULL) delete features[i];
// }
remove_all_children();
// Protest if we are deleting an object while a parent object may still depend on it
// Another possible strategy is to have the child unlist itself from the parent's children
if (parents.size()) {
cvm::log("Warning: destroying " + description + " before its parents objects:");
for (i=0; i<parents.size(); i++) {
cvm::log(parents[i]->description);
}
}
}
void cvm::deps::provide(int feature_id) {
feature_states[feature_id]->available = true;
}
int cvm::deps::enable(int feature_id,
bool dry_run /* default: false */,
// dry_run: fail silently, do not enable if available
// flag is passed recursively to deps of this feature
bool toplevel /* default: true */)
// toplevel: false if this is called as part of a chain of dependency resolution
// this is used to diagnose failed dependencies by displaying the full stack
// only the toplevel dependency will throw a fatal error
{
int res;
size_t i, j;
bool ok;
feature *f = features()[feature_id];
feature_state *fs = feature_states[feature_id];
if (cvm::debug()) {
cvm::log("DEPS: " + description +
(dry_run ? " testing " : " requiring ") +
"\"" + f->description);
}
if (fs->enabled) {
// Do not try to solve deps if already enabled
return COLVARS_OK;
}
if (!fs->available) {
if (!dry_run) {
if (toplevel) {
cvm::error("Error: Feature unavailable: \"" + f->description + "\" in " + description + ".");
} else {
cvm::log("Feature unavailable: \"" + f->description + "\" in " + description);
}
}
return COLVARS_ERROR;
}
// 1) enforce exclusions
for (i=0; i<f->requires_exclude.size(); i++) {
feature *g = features()[f->requires_exclude[i]];
if (cvm::debug())
cvm::log(f->description + " requires exclude " + g->description);
if (is_enabled(f->requires_exclude[i])) {
if (!dry_run) {
cvm::log("Features \"" + f->description + "\" is incompatible with \""
+ g->description + "\" in " + description);
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
}
return COLVARS_ERROR;
}
}
// 2) solve internal deps (self)
for (i=0; i<f->requires_self.size(); i++) {
if (cvm::debug())
cvm::log(f->description + " requires self " + features()[f->requires_self[i]]->description);
res = enable(f->requires_self[i], dry_run, false);
if (res != COLVARS_OK) {
if (!dry_run) {
cvm::log("...required by \"" + f->description + "\" in " + description);
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
}
return res;
}
}
// 3) solve internal alternate deps
for (i=0; i<f->requires_alt.size(); i++) {
// test if one is available; if yes, enable and exit w/ success
ok = false;
for (j=0; j<f->requires_alt[i].size(); j++) {
int g = f->requires_alt[i][j];
if (cvm::debug())
cvm::log(f->description + " requires alt " + features()[g]->description);
res = enable(g, true, false); // see if available
if (res == COLVARS_OK) {
ok = true;
if (!dry_run) enable(g, false, false); // Require again, for real
break;
}
}
if (!ok) {
if (!dry_run) {
cvm::log("No dependency satisfied among alternates:");
cvm::log("-----------------------------------------");
for (j=0; j<f->requires_alt[i].size(); j++) {
int g = f->requires_alt[i][j];
cvm::log(cvm::to_str(j+1) + ". " + features()[g]->description);
cvm::increase_depth();
enable(g, false, false); // Just for printing error output
cvm::decrease_depth();
}
cvm::log("-----------------------------------------");
cvm::log("for \"" + f->description + "\" in " + description);
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
}
return COLVARS_ERROR;
}
}
// 4) solve deps in children
for (i=0; i<f->requires_children.size(); i++) {
int g = f->requires_children[i];
if (cvm::debug())
cvm::log("requires children " + features()[g]->description);
// cvm::log("children " + cvm::to_str(g));
for (j=0; j<children.size(); j++) {
// cvm::log("child " + children[j]->description);
cvm::increase_depth();
res = children[j]->enable(g, dry_run, false);
cvm::decrease_depth();
if (res != COLVARS_OK) {
if (!dry_run) {
cvm::log("...required by \"" + f->description + "\" in " + description);
}
if (toplevel) {
cvm::error("Error: Failed dependency in " + description + ".");
}
return res;
}
}
// If we've just touched the features of child objects, refresh them
if (!dry_run && f->requires_children.size() != 0) {
for (j=0; j<children.size(); j++) {
children[j]->refresh_deps();
}
}
}
// Actually enable feature only once everything checks out
if (!dry_run) fs->enabled = true;
return COLVARS_OK;
}
// disable() {
//
// // we need refs to parents to walk up the deps tree!
// // or refresh
// }
// Shorthand macros for describing dependencies
#define f_description(f, d) features()[f]->description = d
#define f_req_self(f, g) features()[f]->requires_self.push_back(g)
// This macro ensures that exclusions are symmetric
#define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \
features()[g]->requires_exclude.push_back(f)
#define f_req_children(f, g) features()[f]->requires_children.push_back(g)
#define f_req_alt2(f, g, h) features()[f]->requires_alt.push_back(std::vector<int>(2));\
features()[f]->requires_alt.back()[0] = g; \
features()[f]->requires_alt.back()[1] = h
#define f_req_alt3(f, g, h, i) features()[f]->requires_alt.push_back(std::vector<int>(3));\
features()[f]->requires_alt.back()[0] = g; \
features()[f]->requires_alt.back()[1] = h; \
features()[f]->requires_alt.back()[2] = i
void cvm::deps::init_cvb_requires() {
int i;
if (features().size() == 0) {
for (i = 0; i < f_cv_ntot; i++) {
features().push_back(new feature);
}
}
f_description(f_cvb_active, "active");
f_req_children(f_cvb_active, f_cv_active);
f_description(f_cvb_apply_force, "apply force");
f_req_children(f_cvb_apply_force, f_cv_gradient);
f_description(f_cvb_get_system_force, "obtain system force");
f_req_children(f_cvb_get_system_force, f_cv_system_force);
// Initialize feature_states for each instance
feature_states.reserve(f_cvb_ntot);
for (i = 0; i < f_cvb_ntot; i++) {
feature_states.push_back(new feature_state(true, false));
// Most features are available, so we set them so
// and list exceptions below
}
}
void cvm::deps::init_cv_requires() {
size_t i;
if (features().size() == 0) {
for (i = 0; i < f_cv_ntot; i++) {
features().push_back(new feature);
}
f_description(f_cv_active, "active");
f_req_children(f_cv_active, f_cvc_active);
// Colvars must be either a linear combination, or scalar (and polynomial) or scripted
f_req_alt3(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted);
f_description(f_cv_gradient, "gradient");
f_req_children(f_cv_gradient, f_cvc_gradient);
f_description(f_cv_collect_gradient, "collect gradient");
f_req_self(f_cv_collect_gradient, f_cv_gradient);
f_req_self(f_cv_collect_gradient, f_cv_scalar);
f_description(f_cv_fdiff_velocity, "fdiff_velocity");
// System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
f_description(f_cv_system_force, "system force");
f_req_alt2(f_cv_system_force, f_cv_extended_Lagrangian, f_cv_system_force_calc);
// Deps for explicit system force calculation
f_description(f_cv_system_force_calc, "system force calculation");
f_req_self(f_cv_system_force_calc, f_cv_scalar);
f_req_self(f_cv_system_force_calc, f_cv_linear);
f_req_children(f_cv_system_force_calc, f_cvc_inv_gradient);
f_req_self(f_cv_system_force_calc, f_cv_Jacobian);
f_description(f_cv_Jacobian, "Jacobian derivative");
f_req_self(f_cv_Jacobian, f_cv_scalar);
f_req_self(f_cv_Jacobian, f_cv_linear);
f_req_children(f_cv_Jacobian, f_cvc_Jacobian);
f_description(f_cv_hide_Jacobian, "hide Jacobian force");
f_req_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
f_description(f_cv_extended_Lagrangian, "extended Lagrangian");
f_description(f_cv_Langevin, "Langevin dynamics");
f_description(f_cv_linear, "linear");
f_description(f_cv_scalar, "scalar");
f_description(f_cv_output_energy, "output energy");
f_description(f_cv_output_value, "output value");
f_description(f_cv_output_velocity, "output velocity");
f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity);
f_description(f_cv_output_applied_force, "output applied force");
f_description(f_cv_output_system_force, "output system force");
f_req_self(f_cv_output_system_force, f_cv_system_force);
f_description(f_cv_lower_boundary, "lower boundary");
f_req_self(f_cv_lower_boundary, f_cv_scalar);
f_description(f_cv_upper_boundary, "upper boundary");
f_req_self(f_cv_upper_boundary, f_cv_scalar);
f_description(f_cv_grid, "grid");
f_req_self(f_cv_grid, f_cv_lower_boundary);
f_req_self(f_cv_grid, f_cv_upper_boundary);
f_description(f_cv_lower_wall, "lower wall");
f_req_self(f_cv_lower_wall, f_cv_lower_boundary);
f_req_self(f_cv_lower_wall, f_cv_gradient);
f_description(f_cv_upper_wall, "upper wall");
f_req_self(f_cv_upper_wall, f_cv_upper_boundary);
f_req_self(f_cv_upper_wall, f_cv_gradient);
f_description(f_cv_runave, "running average");
f_description(f_cv_corrfunc, "correlation function");
// The features below are set programmatically
f_description(f_cv_scripted, "scripted");
f_description(f_cv_periodic, "periodic");
f_description(f_cv_scalar, "scalar");
f_description(f_cv_linear, "linear");
f_description(f_cv_homogeneous, "homogeneous");
}
// Initialize feature_states for each instance
feature_states.reserve(f_cv_ntot);
for (i = 0; i < f_cv_ntot; i++) {
feature_states.push_back(new feature_state(true, false));
// Most features are available, so we set them so
// and list exceptions below
}
// properties that may NOT be enabled as a dependency
int unavailable_deps[] = {
f_cv_lower_boundary,
f_cv_upper_boundary,
f_cv_extended_Lagrangian,
f_cv_Langevin,
f_cv_scripted,
f_cv_periodic,
f_cv_scalar,
f_cv_linear,
f_cv_homogeneous
};
for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) {
feature_states[unavailable_deps[i]]->available = false;
}
}
void cvm::deps::init_cvc_requires() {
size_t i;
// Initialize static array once and for all
if (features().size() == 0) {
for (i = 0; i < cvm::deps::f_cvc_ntot; i++) {
features().push_back(new feature);
}
f_description(f_cvc_active, "active");
// The dependency below may become useful if we use dynamic atom groups
// f_req_children(f_cvc_active, f_ag_active);
f_description(f_cvc_scalar, "scalar");
f_description(f_cvc_gradient, "gradient");
f_description(f_cvc_inv_gradient, "inverse gradient");
f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
f_description(f_cvc_Jacobian, "Jacobian");
f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
f_description(f_cvc_scalable, "scalable calculation");
f_description(f_cvc_scalable_com, "scalable calculation of centers of mass");
f_req_self(f_cvc_scalable, f_cvc_scalable_com);
// TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group
// f_req_children(f_cvc_scalable, f_ag_scalable);
// f_req_children(f_cvc_scalable_com, f_ag_scalable_com);
}
// Initialize feature_states for each instance
// default as unavailable, not enabled
feature_states.reserve(f_cvc_ntot);
for (i = 0; i < cvm::deps::f_cvc_ntot; i++) {
feature_states.push_back(new feature_state(false, false));
}
// Features that are implemented by all cvcs by default
feature_states[f_cvc_active]->available = true;
feature_states[f_cvc_gradient]->available = true;
// Each cvc specifies what other features are available
feature_states[f_cvc_scalable_com]->available = (proxy->scalable_group_coms() == COLVARS_OK);
feature_states[f_cvc_scalable]->available = feature_states[f_cvc_scalable_com]->available;
}
void cvm::deps::init_ag_requires() {
size_t i;
// Initialize static array once and for all
if (features().size() == 0) {
for (i = 0; i < f_ag_ntot; i++) {
features().push_back(new feature);
}
f_description(f_ag_active, "active");
f_description(f_ag_center, "translational fit");
f_description(f_ag_rotate, "rotational fit");
f_description(f_ag_ref_pos_group, "reference positions group");
f_description(f_ag_fit_gradient_group, "fit gradient for main group");
f_description(f_ag_fit_gradient_ref, "fit gradient for reference group");
f_description(f_ag_atom_forces, "atomic forces");
// parallel calculation implies that we have at least a scalable center of mass,
// but f_ag_scalable is kept as a separate feature to deal with future dependencies
f_description(f_ag_scalable, "scalable group calculation");
f_description(f_ag_scalable_com, "scalable group center of mass calculation");
f_req_self(f_ag_scalable, f_ag_scalable_com);
// f_description(f_ag_min_msd_fit, "minimum MSD fit")
// f_req_self(f_ag_min_msd_fit, f_ag_center)
// f_req_self(f_ag_min_msd_fit, f_ag_rotate)
// f_req_exclude(f_ag_min_msd_fit, f_ag_ref_pos_group)
}
// Initialize feature_states for each instance
// default as unavailable, not enabled
feature_states.reserve(f_ag_ntot);
for (i = 0; i < cvm::deps::f_ag_ntot; i++) {
feature_states.push_back(new feature_state(false, false));
}
// Features that are implemented (or not) by all atom groups
feature_states[f_ag_active]->available = true;
feature_states[f_ag_scalable_com]->available = (proxy->scalable_group_coms() == COLVARS_OK);
feature_states[f_ag_scalable]->available = feature_states[f_ag_scalable_com]->available;
}
void cvm::deps::print_state() {
size_t i;
cvm::log("Enabled features of " + description);
for (i = 0; i < feature_states.size(); i++) {
if (feature_states[i]->enabled)
cvm::log("- " + features()[i]->description);
}
for (i=0; i<children.size(); i++) {
cvm::log("* child " + cvm::to_str(i+1));
cvm::increase_depth();
children[i]->print_state();
cvm::decrease_depth();
}
}
void cvm::deps::add_child(deps *child) {
children.push_back(child);
child->parents.push_back((deps *)this);
}
void cvm::deps::remove_child(deps *child) {
int i;
bool found = false;
for (i = children.size()-1; i>=0; --i) {
if (children[i] == child) {
children.erase(children.begin() + i);
found = true;
break;
}
}
if (!found) {
cvm::error("Trying to remove missing child reference from " + description + "\n");
}
found = false;
for (i = child->parents.size()-1; i>=0; --i) {
if (child->parents[i] == this) {
child->parents.erase(child->parents.begin() + i);
found = true;
break;
}
}
if (!found) {
cvm::error("Trying to remove missing parent reference from " + child->description + "\n");
}
}
void cvm::deps::remove_all_children() {
size_t i;
int j;
bool found;
for (i = 0; i < children.size(); ++i) {
found = false;
for (j = children[i]->parents.size()-1; j>=0; --j) {
if (children[i]->parents[j] == this) {
children[i]->parents.erase(children[i]->parents.begin() + j);
found = true;
break;
}
}
if (!found) {
cvm::error("Trying to remove missing parent reference from " + children[i]->description + "\n");
}
}
children.clear();
}

265
lib/colvars/colvardeps.h Normal file
View File

@ -0,0 +1,265 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#ifndef COLVARDEPS_H
#define COLVARDEPS_H
/// Parent class for a member object of a bias, cv or cvc etc. containing dependencies
/// (features) and handling dependency resolution
// Some features like colvar::f_linear have no dependencies, require() doesn't enable anything but fails if unavailable
// Policy: those features are unavailable at all times
// Other features are under user control
// They are unavailable unless requested by the user, then they may be enabled subject to
// satisfied dependencies
// It seems important to have available default to false (for safety) and enabled to false (for efficiency)
class cvm::deps {
public:
deps() {}
virtual ~deps();
// Subclasses should initialize the following members:
std::string description; // reference to object name (cv, cvc etc.)
/// This contains the current state of each feature for each object
struct feature_state {
feature_state(bool a, bool e)
: available(a), enabled(e) {}
/// Available means: supported, subject to dependencies as listed,
/// MAY BE ENABLED AS A RESULT OF DEPENDENCY SOLVING
/// Remains false for passive flags that are set based on other object properties,
/// eg. f_cv_linear
/// Is set to true upon user request for features that are implemented by the user
/// or under his/her direct control, e.g. f_cv_scripted or f_cv_extended_Lagrangian
bool available;
/// Currently enabled - this flag is subject to change dynamically
/// TODO consider implications for dependency solving: anyone who disables
/// it should trigger a refresh of parent objects
bool enabled; // see if this should be private depending on implementation
// bool enabledOnce; // this should trigger an update when object is evaluated
};
/// List of the state of all features
std::vector<feature_state *> feature_states;
/// Describes a feature and its dependecies
/// used in a static array within each subclass
class feature {
public:
feature() {}
~feature() {}
std::string description; // Set by derived object initializer
// features that this feature requires in the same object
// NOTE: we have no safety mechanism against circular dependencies, however, they would have to be internal to an object (ie. requires_self or requires_alt)
std::vector<int> requires_self;
// Features that are incompatible, ie. required disabled
// if enabled, they will cause a dependency failure (they will not be disabled)
// To enforce these dependencies regardless of the order in which they
// are enabled, they are always set in a symmetric way, so whichever is enabled
// second will cause the dependency to fail
std::vector<int> requires_exclude;
// sets of features that are required in an alternate way
// when parent feature is enabled, if none are enabled, the first one listed that is available will be enabled
std::vector<std::vector<int> > requires_alt;
// features that this feature requires in children
std::vector<int> requires_children;
};
// Accessor to array of all features with deps, static in most derived classes
// Subclasses with dynamic dependency trees may override this
// with a non-static array
// Intermediate classes (colvarbias and colvarcomp, which are also base classes)
// implement this as virtual to allow overriding
virtual std::vector<feature *>&features() = 0;
void add_child(deps *child);
void remove_child(deps *child);
/// Used before deleting an object, if not handled by that object's destructor
/// (useful for cvcs because their children are member objects)
void remove_all_children();
private:
// pointers to objects this object depends on
// list should be maintained by any code that modifies the object
// this could be secured by making lists of colvars / cvcs / atom groups private and modified through accessor functions
std::vector<deps *> children;
// pointers to objects that depend on this object
// the size of this array is in effect a reference counter
std::vector<deps *> parents;
public:
// disabling a feature f:
// if parents depend on f, tell them to refresh / check that they are ok?
// if children provide features to satisfy f ONLY, disable that
// When the state of this object has changed, recursively tell parents
// to enforce their dependencies
// void refresh_parents() {
//
// }
// std::vector<deps *> parents; // Needed to trigger a refresh if capabilities of this object change
// End of members to be initialized by subclasses
// Checks whether given feature is enabled
// Defaults to querying f_*_active
inline bool is_enabled(int f = f_cv_active) const {
return feature_states[f]->enabled;
}
// Checks whether given feature is available
// Defaults to querying f_*_active
inline bool is_available(int f = f_cv_active) const {
return feature_states[f]->available;
}
void provide(int feature_id); // set the feature's flag to available in local object
int enable(int f, bool dry_run = false, bool toplevel = true); // enable a feature and recursively solve its dependencies
// dry_run is set to true to recursively test if a feature is available, without enabling it
// int disable(int f);
/// This function is called whenever feature states are changed outside
/// of the object's control, that is, by parents
/// Eventually it may also be used when properties of children change
virtual int refresh_deps() { return COLVARS_OK; }
// NOTE that all feature enums should start with f_*_active
enum features_biases {
/// \brief Bias is active
f_cvb_active,
f_cvb_apply_force,
f_cvb_get_system_force,
f_cvb_ntot
};
enum features_colvar {
/// \brief Calculate colvar
f_cv_active,
/// \brief Gradients are calculated and temporarily stored, so
/// that external forces can be applied
f_cv_gradient,
/// \brief Collect atomic gradient data from all cvcs into vector
/// atomic_gradient
f_cv_collect_gradient,
/// \brief Calculate the velocity with finite differences
f_cv_fdiff_velocity,
/// \brief The system force is calculated, projecting the atomic
/// forces on the inverse gradient
f_cv_system_force,
/// \brief Calculate system force from atomic forces
f_cv_system_force_calc,
/// \brief Estimate Jacobian derivative
f_cv_Jacobian,
/// \brief Do not report the Jacobian force as part of the system force
/// instead, apply a correction internally to cancel it
f_cv_hide_Jacobian,
/// \brief The variable has a harmonic restraint around a moving
/// center with fictitious mass; bias forces will be applied to
/// the center
f_cv_extended_Lagrangian,
/// \brief The extended system coordinate undergoes Langevin dynamics
f_cv_Langevin,
/// \brief Output the potential and kinetic energies
/// (for extended Lagrangian colvars only)
f_cv_output_energy,
/// \brief Output the value to the trajectory file (on by default)
f_cv_output_value,
/// \brief Output the velocity to the trajectory file
f_cv_output_velocity,
/// \brief Output the applied force to the trajectory file
f_cv_output_applied_force,
/// \brief Output the system force to the trajectory file
f_cv_output_system_force,
/// \brief A lower boundary is defined
f_cv_lower_boundary,
/// \brief An upper boundary is defined
f_cv_upper_boundary,
/// \brief Provide a discretization of the values of the colvar to
/// be used by the biases or in analysis (needs lower and upper
/// boundary)
f_cv_grid,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes below the lower wall
f_cv_lower_wall,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes above the upper wall
f_cv_upper_wall,
/// \brief Compute running average
f_cv_runave,
/// \brief Compute time correlation function
f_cv_corrfunc,
/// \brief Value and gradient computed by user script
f_cv_scripted,
/// \brief Colvar is periodic
f_cv_periodic,
/// \brief is scalar
f_cv_scalar,
f_cv_linear,
f_cv_homogeneous,
/// \brief Number of colvar features
f_cv_ntot
};
enum features_cvc {
f_cvc_active,
f_cvc_scalar,
f_cvc_gradient,
f_cvc_inv_gradient,
/// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
f_cvc_debug_gradient,
f_cvc_Jacobian,
f_cvc_com_based,
f_cvc_scalable,
f_cvc_scalable_com,
f_cvc_ntot
};
enum features_atomgroup {
f_ag_active,
f_ag_center,
f_ag_rotate,
f_ag_ref_pos_group,
/// Perform a standard minimum msd fit for given atoms
/// ie. not using refpositionsgroup
// f_ag_min_msd_fit,
f_ag_fit_gradient_group,// TODO check that these are sometimes needed separately
// maybe for minimum RMSD?
f_ag_fit_gradient_ref,
f_ag_atom_forces,
f_ag_scalable,
f_ag_scalable_com,
f_ag_ntot
};
void init_cvb_requires();
void init_cv_requires();
void init_cvc_requires();
void init_ag_requires();
/// \brief print all enabled features and those of children, for debugging
void print_state();
};
#endif