ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

import colvars code with multi-thread support

Browse Source
This commit is contained in:
Axel Kohlmeyer
2016-04-08 09:30:28 -04:00
parent 247e55937a
commit 5a26b927bc
34 changed files with 2534 additions and 1453 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/colvars/Makefile.g++
View File

@ -18,7 +18,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
LIB = libcolvars.a

2
lib/colvars/Makefile.mingw32-cross
View File

@ -21,7 +21,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
DIR = Obj_mingw32/

2
lib/colvars/Makefile.mingw64-cross
View File

@ -21,7 +21,7 @@ SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \
colvarbias_histogram.cpp colvarbias_meta.cpp colvarbias_restraint.cpp \
colvarcomp_angles.cpp colvarcomp_coordnums.cpp colvarcomp.cpp \
colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \
colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvardeps.cpp colvar.cpp colvargrid.cpp colvarmodule.cpp colvarparse.cpp \
colvarscript.cpp colvartypes.cpp colvarvalue.cpp
DIR = Obj_mingw64/

911
lib/colvars/colvar.cpp
View File

File diff suppressed because it is too large Load Diff

162
lib/colvars/colvar.h
View File

@ -10,6 +10,7 @@
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvardeps.h"
/// \brief A collective variable (main class); to be defined, it needs
@ -37,7 +38,7 @@
/// \link colvarvalue \endlink type, you should also add its
/// initialization line in the \link colvar \endlink constructor.
class colvar : public colvarparse {
class colvar : public colvarparse, public cvm::deps {
public:
@ -63,8 +64,6 @@ public:
/// subtracted.
colvarvalue const & system_force() const;
/// \brief
/// \brief Typical fluctuation amplitude for this collective
/// variable (e.g. local width of a free energy basin)
///
@ -75,87 +74,15 @@ public:
/// variable.
cvm::real width;
/// \brief True if this \link colvar \endlink is a linear
/// combination of \link cvc \endlink elements
bool b_linear;
/// \brief Implementation of the feature list for colvar
static std::vector<feature *> cv_features;
/// \brief True if this \link colvar \endlink is a linear
/// combination of cvcs with coefficients 1 or -1
bool b_homogeneous;
/// \brief Implementation of the feature list accessor for colvar
std::vector<feature *> &features() {
return cv_features;
}
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating inverse gradients
bool b_inverse_gradients;
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating Jacobian forces
bool b_Jacobian_force;
/// \brief Options controlling the behaviour of the colvar during
/// the simulation, which are set from outside the cvcs
enum task {
/// \brief Gradients are calculated and temporarily stored, so
/// that external forces can be applied
task_gradients,
/// \brief Collect atomic gradient data from all cvcs into vector
/// atomic_gradients
task_collect_gradients,
/// \brief Calculate the velocity with finite differences
task_fdiff_velocity,
/// \brief The system force is calculated, projecting the atomic
/// forces on the inverse gradients
task_system_force,
/// \brief The variable has a harmonic restraint around a moving
/// center with fictitious mass; bias forces will be applied to
/// the center
task_extended_lagrangian,
/// \brief The extended system coordinate undergoes Langevin
/// dynamics
task_langevin,
/// \brief Output the potential and kinetic energies
/// (for extended Lagrangian colvars only)
task_output_energy,
/// \brief Compute analytically the "force" arising from the
/// geometric entropy component (for example, that from the angular
/// states orthogonal to a distance vector)
task_Jacobian_force,
/// \brief Report the Jacobian force as part of the system force
/// if disabled, apply a correction internally to cancel it
task_report_Jacobian_force,
/// \brief Output the value to the trajectory file (on by default)
task_output_value,
/// \brief Output the velocity to the trajectory file
task_output_velocity,
/// \brief Output the applied force to the trajectory file
task_output_applied_force,
/// \brief Output the system force to the trajectory file
task_output_system_force,
/// \brief A lower boundary is defined
task_lower_boundary,
/// \brief An upper boundary is defined
task_upper_boundary,
/// \brief Provide a discretization of the values of the colvar to
/// be used by the biases or in analysis (needs lower and upper
/// boundary)
task_grid,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes below the lower wall
task_lower_wall,
/// \brief Apply a restraining potential (|x-xb|^2) to the colvar
/// when it goes above the upper wall
task_upper_wall,
/// \brief Compute running average
task_runave,
/// \brief Compute time correlation function
task_corrfunc,
/// \brief Value and gradients computed by user script
task_scripted,
/// \brief Number of possible tasks
task_ntot
};
/// Tasks performed by this colvar
bool tasks[task_ntot];
int refresh_deps();
/// List of biases that depend on this colvar
std::vector<colvarbias *> biases;
@ -213,7 +140,7 @@ protected:
/// Cached reported velocity
colvarvalue v_reported;
// Options for task_extended_lagrangian
// Options for extended_lagrangian
/// Restraint center
colvarvalue xr;
/// Velocity of the restraint center
@ -230,7 +157,7 @@ protected:
/// \brief Harmonic restraint force
colvarvalue fr;
/// \brief Jacobian force, when task_Jacobian_force is enabled
/// \brief Jacobian force, when Jacobian_force is enabled
colvarvalue fj;
/// Cached reported system force
@ -243,7 +170,7 @@ public:
/// the biases are updated
colvarvalue fb;
/// \brief Total \em applied force; fr (if task_extended_lagrangian
/// \brief Total \em applied force; fr (if extended_lagrangian
/// is defined), fb (if biases are applied) and the walls' forces
/// (if defined) contribute to it
colvarvalue f;
@ -293,12 +220,6 @@ public:
/// Constructor
colvar(std::string const &conf);
/// Enable the specified task
int enable(colvar::task const &t);
/// Disable the specified task
void disable(colvar::task const &t);
/// Get ready for a run and possibly re-initialize internal data
void setup();
@ -306,16 +227,46 @@ public:
~colvar();
/// \brief Calculate the colvar value and all the other requested
/// quantities
void calc();
/// \brief Calculate the colvar's value and related quantities
int calc();
/// \brief Calculate a given subset of colvar components (CVCs) (default: all CVCs)
/// \brief Calculate a subset of the colvar components (CVCs) currently active
/// (default: all active CVCs)
/// Note: both arguments refer to the sect of *active* CVCs, not all CVCs
int calc_cvcs(int first = 0, size_t num_cvcs = 0);
/// Ensure that the selected range of CVCs is consistent
int check_cvc_range(int first_cvc, size_t num_cvcs);
/// \brief Calculate the values of the given subset of CVCs
int calc_cvc_values(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients
int calc_cvc_gradients(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for system forces
int calc_cvc_sys_forces(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces
int calc_cvc_Jacobians(int first, size_t num_cvcs);
/// \brief Collect quantities from CVCs and update aggregated data for the colvar
int collect_cvc_data();
/// \brief Collect the values of the CVCs
int collect_cvc_values();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients
int collect_cvc_gradients();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for system forces
int collect_cvc_sys_forces();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces
int collect_cvc_Jacobians();
/// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
int calc_colvar_properties();
/// Get the current biasing force
inline colvarvalue bias_force() const
{
return fb;
}
/// Set the total biasing force to zero
void reset_bias_force();
@ -326,7 +277,7 @@ public:
/// equations of motion of internal degrees of freedom; see also
/// colvar::communicate_forces()
/// return colvar energy if extended Lagrandian active
cvm::real update();
cvm::real update_forces_energy();
/// \brief Communicate forces (previously calculated in
/// colvar::update()) to the external degrees of freedom
@ -335,11 +286,19 @@ public:
/// \brief Enables and disables individual CVCs based on flags
int set_cvc_flags(std::vector<bool> const &flags);
/// \brief Updates the flags in the CVC objects
/// \brief Updates the flags in the CVC objects, and their number
int update_cvc_flags();
/// \brief Return the number of CVC objects with an active flag
int num_active_cvcs() const;
protected:
/// \brief Number of CVC objects with an active flag
size_t n_active_cvcs;
/// Sum of square coefficients for active cvcs
cvm::real active_cvc_square_norm;
public:
/// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
inline size_t num_active_cvcs() const { return n_active_cvcs; }
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
@ -528,7 +487,7 @@ protected:
std::vector<bool> cvc_flags;
/// \brief Initialize the sorted list of atom IDs for atoms involved
/// in all cvcs (called when enabling task_collect_gradients)
/// in all cvcs (called when enabling f_cv_collect_gradients)
void build_atom_list(void);
private:
@ -580,6 +539,9 @@ inline colvarvalue const & colvar::system_force() const
inline void colvar::add_bias_force(colvarvalue const &force)
{
if (cvm::debug()) {
cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
}
fb += force;
}

69
lib/colvars/colvaratoms.cpp
View File

@ -4,7 +4,6 @@
#include "colvarparse.h"
#include "colvaratoms.h"
cvm::atom::atom()
{
index = -1;
@ -68,8 +67,7 @@ cvm::atom_group::atom_group(std::string const &conf,
char const *key_in)
{
key = key_in;
cvm::log("Defining atom group \""+
std::string(key)+"\".\n");
cvm::log("Defining atom group \"" + key + "\".\n");
init();
// real work is done by parse
parse(conf);
@ -93,8 +91,9 @@ cvm::atom_group::atom_group()
cvm::atom_group::~atom_group()
{
if (index >= 0) {
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->clear_atom_group(index);
index = -1;
}
if (ref_pos_group) {
@ -151,7 +150,7 @@ int cvm::atom_group::add_atom_id(int aid)
int cvm::atom_group::remove_atom(cvm::atom_iter ai)
{
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
@ -173,10 +172,11 @@ int cvm::atom_group::remove_atom(cvm::atom_iter ai)
int cvm::atom_group::init()
{
if (!key.size()) key = "atoms";
atoms.clear();
b_scalable = false;
// TODO: check with proxy whether atom forces etc are available
init_ag_requires();
index = -1;
b_center = false;
@ -216,7 +216,7 @@ void cvm::atom_group::update_total_mass()
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
total_mass = (cvm::proxy)->get_atom_group_mass(index);
} else {
total_mass = 0.0;
@ -243,7 +243,7 @@ void cvm::atom_group::update_total_charge()
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
total_charge = (cvm::proxy)->get_atom_group_charge(index);
} else {
total_charge = 0.0;
@ -256,6 +256,7 @@ void cvm::atom_group::update_total_charge()
int cvm::atom_group::parse(std::string const &conf)
{
colvarproxy *proxy = cvm::proxy;
std::string group_conf;
// TODO move this to the cvc class constructor/init
@ -280,6 +281,8 @@ int cvm::atom_group::parse(std::string const &conf)
cvm::log("Initializing atom group \""+key+"\".\n");
description = "atom group " + key;
// whether or not to include messages in the log
// colvarparse::Parse_Mode mode = parse_silent;
// {
@ -291,9 +294,6 @@ int cvm::atom_group::parse(std::string const &conf)
int parse_error = COLVARS_OK;
// if the cvc allows it, the flag has been set to true by default
get_keyval(group_conf, "scalable", b_scalable, b_scalable);
{
std::string numbers_conf = "";
size_t pos = 0;
@ -373,9 +373,9 @@ int cvm::atom_group::parse(std::string const &conf)
}
// Catch any errors from all the initialization steps above
if (parse_error || cvm::get_error()) return COLVARS_ERROR;
if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
@ -449,7 +449,7 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
if (atom_indexes.size()) {
atoms_ids.reserve(atoms_ids.size()+atom_indexes.size());
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
for (size_t i = 0; i < atom_indexes.size(); i++) {
add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i]));
}
@ -490,7 +490,7 @@ int cvm::atom_group::add_index_group(std::string const &index_group_name)
atoms_ids.reserve(atoms_ids.size()+index_groups_i->size());
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
for (size_t i = 0; i < index_groups_i->size(); i++) {
add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i]));
}
@ -520,7 +520,7 @@ int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf)
atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
for (int anum = initial; anum <= final; anum++) {
add_atom_id((cvm::proxy)->check_atom_id(anum));
}
@ -558,7 +558,7 @@ int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
for (int resid = initial; resid <= final; resid++) {
add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid));
}
@ -754,7 +754,7 @@ int cvm::atom_group::calc_required_properties()
// TODO check if the com is needed?
calc_center_of_mass();
if (!b_scalable) {
if (!is_enabled(f_ag_scalable)) {
// TODO check if calc_center_of_geometry() is needed without a fit?
calc_center_of_geometry();
if (b_center || b_rotate) {
@ -805,7 +805,7 @@ void cvm::atom_group::apply_translation(cvm::rvector const &t)
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR);
return;
}
@ -822,7 +822,7 @@ void cvm::atom_group::apply_rotation(cvm::rotation const &rot)
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: cannot rotate the coordinates of a scalable atom group.\n", INPUT_ERROR);
return;
}
@ -852,6 +852,8 @@ void cvm::atom_group::read_velocities()
}
}
// TODO make this a calc function
void cvm::atom_group::read_system_forces()
{
if (b_dummy) return;
@ -891,7 +893,7 @@ int cvm::atom_group::calc_center_of_mass()
{
if (b_dummy) {
com = dummy_atom_pos;
} else if (b_scalable) {
} else if (is_enabled(f_ag_scalable)) {
com = (cvm::proxy)->get_atom_group_com(index);
} else {
com.reset();
@ -922,7 +924,7 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
{
if (b_dummy) return;
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
scalar_com_gradient = grad;
return;
}
@ -993,7 +995,7 @@ std::vector<cvm::atom_pos> cvm::atom_group::positions() const
"from a dummy atom group.\n", INPUT_ERROR);
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic positions are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
@ -1014,7 +1016,7 @@ std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted(cvm::rvector const
"from a dummy atom group.\n", INPUT_ERROR);
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic positions are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
@ -1035,7 +1037,7 @@ std::vector<cvm::rvector> cvm::atom_group::velocities() const
"from a dummy atom group.\n", INPUT_ERROR);
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic velocities are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
@ -1056,7 +1058,7 @@ std::vector<cvm::rvector> cvm::atom_group::system_forces() const
"from a dummy atom group.\n", INPUT_ERROR);
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
cvm::error("Error: atomic system forces are not available "
"from a scalable atom group.\n", INPUT_ERROR);
}
@ -1070,6 +1072,8 @@ std::vector<cvm::rvector> cvm::atom_group::system_forces() const
return f;
}
// TODO make this an accessor
cvm::rvector cvm::atom_group::system_force() const
{
if (b_dummy) {
@ -1077,7 +1081,7 @@ cvm::rvector cvm::atom_group::system_force() const
"from a dummy atom group.\n", INPUT_ERROR);
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
return (cvm::proxy)->get_atom_group_system_force(index);
}
@ -1104,7 +1108,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
return;
}
@ -1156,7 +1160,7 @@ void cvm::atom_group::apply_force(cvm::rvector const &force)
return;
}
if (b_scalable) {
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force);
}
@ -1209,3 +1213,8 @@ void cvm::atom_group::apply_forces(std::vector<cvm::rvector> const &forces)
}
}
// Static members
std::vector<cvm::deps::feature *> cvm::atom_group::ag_features;

12
lib/colvars/colvaratoms.h
View File

@ -5,6 +5,7 @@
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvardeps.h"
/// \brief Stores numeric id, mass and all mutable data for an atom,
@ -138,7 +139,7 @@ public:
/// \brief Group of \link atom \endlink objects, mostly used by a
/// \link cvc \endlink object to gather all atomic data
class colvarmodule::atom_group
: public colvarparse
: public colvarparse, public cvm::deps
{
public:
@ -187,8 +188,13 @@ public:
/// change atom masses after their initialization.
void reset_mass(std::string &name, int i, int j);
/// \brief Whether or not the properties of this group will be computed in parallel
bool b_scalable;
/// \brief Implementation of the feature list for atom group
static std::vector<feature *> ag_features;
/// \brief Implementation of the feature list accessor for atom group
virtual std::vector<feature *> &features() {
return ag_features;
}
/// \brief Default constructor
atom_group();

23
lib/colvars/colvarbias.cpp
View File

@ -10,6 +10,8 @@ colvarbias::colvarbias(std::string const &conf, char const *key)
{
cvm::log("Initializing a new \""+std::string(key)+"\" instance.\n");
init_cvb_requires();
size_t rank = 1;
std::string const key_str(key);
@ -35,6 +37,8 @@ colvarbias::colvarbias(std::string const &conf, char const *key)
return;
}
description = "bias " + name;
// lookup the associated colvars
std::vector<std::string> colvars_str;
if (get_keyval(conf, "colvars", colvars_str)) {
@ -46,6 +50,13 @@ colvarbias::colvarbias(std::string const &conf, char const *key)
cvm::error("Error: no collective variables specified.\n");
return;
}
for (size_t i=0; i<colvars.size(); i++) {
// All biases need at least the value of colvars
// although possibly not at all timesteps
add_child(colvars[i]);
}
// Start in active state by default
enable(f_cvb_active);
get_keyval(conf, "outputEnergy", b_output_energy, false);
}
@ -84,7 +95,8 @@ colvarbias::~colvarbias()
void colvarbias::add_colvar(std::string const &cv_name)
{
if (colvar *cv = cvm::colvar_by_name(cv_name)) {
cv->enable(colvar::task_gradients);
// Removed this as nor all biases apply forces eg histogram
// cv->enable(colvar::task_gradients);
if (cvm::debug())
cvm::log("Applying this bias to collective variable \""+
cv->name+"\".\n");
@ -100,10 +112,11 @@ void colvarbias::add_colvar(std::string const &cv_name)
}
cvm::real colvarbias::update()
int colvarbias::update()
{
// Note: if anything is added here, it should be added also in the SMP block of calc_biases()
has_data = true;
return 0.0;
return COLVARS_OK;
}
@ -172,3 +185,7 @@ std::ostream & colvarbias::write_traj(std::ostream &os)
<< bias_energy;
return os;
}
// Static members
std::vector<cvm::deps::feature *> colvarbias::cvb_features;

13
lib/colvars/colvarbias.h
View File

@ -5,10 +5,11 @@
#include "colvar.h"
#include "colvarparse.h"
#include "colvardeps.h"
/// \brief Collective variable bias, base class
class colvarbias : public colvarparse {
class colvarbias : public colvarparse, public cvm::deps {
public:
/// Name of this bias
@ -19,7 +20,7 @@ public:
/// Retrieve colvar values and calculate their biasing forces
/// Return bias energy
virtual cvm::real update();
virtual int update();
// TODO: move update_bias here (share with metadynamics)
@ -77,6 +78,14 @@ public:
inline cvm::real get_energy() {
return bias_energy;
}
/// \brief Implementation of the feature list for colvarbias
static std::vector<feature *> cvb_features;
/// \brief Implementation of the feature list accessor for colvarbias
virtual std::vector<feature *> &features() {
return cvb_features;
}
protected:
/// \brief Pointers to collective variables to which the bias is

41
lib/colvars/colvarbias_abf.cpp
View File

@ -58,37 +58,22 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
cvm::error("Error: no collective variables specified for the ABF bias.\n");
}
if (update_bias) {
// Request calculation of system force (which also checks for availability)
enable(f_cvb_get_system_force);
}
if (apply_bias) {
enable(f_cvb_apply_force);
}
for (size_t i = 0; i < colvars.size(); i++) {
if (colvars[i]->value().type() != colvarvalue::type_scalar) {
cvm::error("Error: ABF bias can only use scalar-type variables.\n");
}
colvars[i]->enable(colvar::task_gradients);
if (update_bias) {
// Request calculation of system force (which also checks for availability)
colvars[i]->enable(colvar::task_system_force);
if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) {
// request computation of Jacobian force
// ultimately, will be done regardless of extended Lagrangian
// and colvar should then just report zero Jacobian force
colvars[i]->enable(colvar::task_Jacobian_force);
// request Jacobian force as part as system force
// except if the user explicitly requires the "silent" Jacobian
// correction AND the colvar has a single component
if (hide_Jacobian) {
if (colvars[i]->n_components() > 1) {
cvm::log("WARNING: colvar \"" + colvars[i]->name
+ "\" has multiple components; reporting its Jacobian forces\n");
colvars[i]->enable(colvar::task_report_Jacobian_force);
}
} else {
colvars[i]->enable(colvar::task_report_Jacobian_force);
}
}
colvars[i]->enable(f_cv_grid);
if (hide_Jacobian) {
colvars[i]->enable(f_cv_hide_Jacobian);
}
// Here we could check for orthogonality of the Cartesian coordinates
@ -176,7 +161,7 @@ colvarbias_abf::~colvarbias_abf()
/// Update the FE gradient, compute and apply biasing force
/// also output data to disk if needed
cvm::real colvarbias_abf::update()
int colvarbias_abf::update()
{
if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
@ -283,7 +268,7 @@ cvm::real colvarbias_abf::update()
cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
}
return 0.0;
return COLVARS_OK;
}
int colvarbias_abf::replica_share() {

2
lib/colvars/colvarbias_abf.h
View File

@ -22,7 +22,7 @@ public:
colvarbias_abf(std::string const &conf, char const *key);
~colvarbias_abf();
cvm::real update();
int update();
private:

6
lib/colvars/colvarbias_alb.cpp
View File

@ -41,6 +41,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
current_coupling.resize(colvars.size());
coupling_rate.resize(colvars.size());
enable(f_cvb_apply_force);
for (i = 0; i < colvars.size(); i++) {
colvar_centers[i].type(colvars[i]->value());
@ -120,7 +121,7 @@ colvarbias_alb::~colvarbias_alb() {
}
cvm::real colvarbias_alb::update() {
int colvarbias_alb::update() {
bias_energy = 0.0;
update_calls++;
@ -224,8 +225,7 @@ cvm::real colvarbias_alb::update() {
}
return bias_energy;
return COLVARS_OK;
}

2
lib/colvars/colvarbias_alb.h
View File

@ -13,7 +13,7 @@ public:
virtual ~colvarbias_alb();
virtual cvm::real update();
virtual int update();
/// Read the bias configuration from a restart file
virtual std::istream & read_restart(std::istream &is);

14
lib/colvars/colvarbias_histogram.cpp
View File

@ -10,6 +10,8 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
: colvarbias(conf, key),
grid(NULL), out_name("")
{
size_t i;
get_keyval(conf, "outputFile", out_name, std::string(""));
get_keyval(conf, "outputFileDX", out_name_dx, std::string(""));
get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
@ -21,7 +23,6 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
colvar_array_size = 0;
{
size_t i;
bool colvar_array = false;
get_keyval(conf, "gatherVectorColvars", colvar_array, colvar_array);
@ -59,6 +60,10 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent);
}
for (i = 0; i < colvars.size(); i++) {
colvars[i]->enable(f_cv_grid);
}
grid = new colvar_grid_scalar();
{
@ -86,10 +91,11 @@ colvarbias_histogram::~colvarbias_histogram()
}
/// Update the grid
cvm::real colvarbias_histogram::update()
int colvarbias_histogram::update()
{
int error_code = COLVARS_OK;
// update base class
colvarbias::update();
error_code |= colvarbias::update();
if (cvm::debug()) {
cvm::log("Updating histogram bias " + this->name);
@ -142,7 +148,7 @@ cvm::real colvarbias_histogram::update()
write_output_files();
}
return 0.0; // no bias energy for histogram
return error_code | cvm::get_error();
}

2
lib/colvars/colvarbias_histogram.h
View File

@ -19,7 +19,7 @@ public:
colvarbias_histogram(std::string const &conf, char const *key);
~colvarbias_histogram();
cvm::real update();
int update();
int write_output_files();

12
lib/colvars/colvarbias_meta.cpp
View File

@ -59,6 +59,9 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
comm = single_replica;
}
// This implies gradients for all colvars
enable(f_cvb_apply_force);
get_keyval(conf, "useGrids", use_grids, true);
if (use_grids) {
@ -68,6 +71,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
expand_grids = false;
size_t i;
for (i = 0; i < colvars.size(); i++) {
colvars[i]->enable(f_cv_grid);
if (colvars[i]->expand_boundaries) {
expand_grids = true;
cvm::log("Metadynamics bias \""+this->name+"\""+
@ -82,10 +86,6 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false);
for (i = 0; i < colvars.size(); i++) {
colvars[i]->enable(colvar::task_grid);
}
hills_energy = new colvar_grid_scalar(colvars);
hills_energy_gradients = new colvar_grid_gradient(colvars);
} else {
@ -253,7 +253,7 @@ colvarbias_meta::delete_hill(hill_iter &h)
}
cvm::real colvarbias_meta::update()
int colvarbias_meta::update()
{
if (cvm::debug())
cvm::log("Updating the metadynamics bias \""+this->name+"\""+
@ -544,7 +544,7 @@ cvm::real colvarbias_meta::update()
", hills forces = "+cvm::to_str(colvar_forces)+".\n");
}
return bias_energy;
return COLVARS_OK;
}

4
lib/colvars/colvarbias_meta.h
View File

@ -36,7 +36,7 @@ public:
/// Destructor
virtual ~colvarbias_meta();
virtual cvm::real update();
virtual int update();
virtual std::istream & read_restart(std::istream &is);
@ -278,7 +278,7 @@ public:
W(W_in),
centers(cv.size()),
widths(cv.size()),
it(cvm::it),
it(cvm::step_absolute()),
replica(replica_in)
{
for (size_t i = 0; i < cv.size(); i++) {

7
lib/colvars/colvarbias_restraint.cpp
View File

@ -18,6 +18,9 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf,
colvar_centers.resize(colvars.size());
colvar_centers_raw.resize(colvars.size());
size_t i;
enable(f_cvb_apply_force);
for (i = 0; i < colvars.size(); i++) {
colvar_centers[i].type(colvars[i]->value());
colvar_centers_raw[i].type(colvars[i]->value());
@ -167,7 +170,7 @@ cvm::real colvarbias_restraint::energy_difference(std::string const &conf)
}
cvm::real colvarbias_restraint::update()
int colvarbias_restraint::update()
{
bias_energy = 0.0;
@ -333,7 +336,7 @@ cvm::real colvarbias_restraint::update()
cvm::log("Current forces for the restraint bias \""+
this->name+"\": "+cvm::to_str(colvar_forces)+".\n");
return bias_energy;
return COLVARS_OK;
}

2
lib/colvars/colvarbias_restraint.h
View File

@ -12,7 +12,7 @@ class colvarbias_restraint : public colvarbias {
public:
/// Retrieve colvar values and calculate their biasing forces
virtual cvm::real update();
virtual int update();
/// Load new configuration - force constant and/or centers only
virtual void change_configuration(std::string const &conf);

174
lib/colvars/colvarcomp.cpp
View File

@ -10,26 +10,22 @@
colvar::cvc::cvc()
: sup_coeff(1.0),
sup_np(1),
b_enabled(true),
b_periodic(false),
b_inverse_gradients(false),
b_Jacobian_derivative(false),
b_debug_gradients(false)
b_try_scalable(true)
{}
colvar::cvc::cvc(std::string const &conf)
: sup_coeff(1.0),
sup_np(1),
b_enabled(true),
b_periodic(false),
b_inverse_gradients(false),
b_Jacobian_derivative(false),
b_debug_gradients(false)
b_try_scalable(true)
{
if (cvm::debug())
cvm::log("Initializing cvc base object.\n");
init_cvc_requires();
get_keyval(conf, "name", this->name, std::string(""), parse_silent);
get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
@ -38,38 +34,109 @@ colvar::cvc::cvc(std::string const &conf)
get_keyval(conf, "period", period, 0.0);
get_keyval(conf, "wrapAround", wrap_center, 0.0);
get_keyval(conf, "debugGradients", b_debug_gradients, false, parse_silent);
{
bool b_debug_gradient;
get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent);
if (b_debug_gradient) enable(f_cvc_debug_gradient);
}
// Attempt scalable calculations when in parallel? (By default yes, if available)
get_keyval(conf, "scalable", b_try_scalable, true);
if (cvm::debug())
cvm::log("Done initializing cvc base object.\n");
}
void colvar::cvc::parse_group(std::string const &conf,
char const *group_key,
cvm::atom_group &group,
bool optional)
cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
char const *group_key,
bool optional)
{
cvm::atom_group *group = NULL;
if (key_lookup(conf, group_key)) {
// TODO turn on scalable flag for group objects in cvc init function
group.key = group_key;
if (group.parse(conf) != COLVARS_OK) {
group = new cvm::atom_group;
group->key = group_key;
if (b_try_scalable) {
if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
enable(f_cvc_scalable_com);
enable(f_cvc_scalable);
group->enable(f_ag_scalable_com);
group->enable(f_ag_scalable);
}
// TODO check for other types of parallelism here
if (is_enabled(f_cvc_scalable)) {
cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
} else {
cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n");
}
}
if (group->parse(conf) == COLVARS_OK) {
atom_groups.push_back(group);
} else {
cvm::error("Error parsing definition for atom group \""+
std::string(group_key)+"\".\n");
return;
}
} else {
if (! optional) {
cvm::error("Error: definition for atom group \""+
std::string(group_key)+"\" not found.\n");
return;
}
}
return group;
}
int colvar::cvc::setup()
{
size_t i;
description = "cvc " + name;
for (i = 0; i < atom_groups.size(); i++) {
add_child((cvm::deps *) atom_groups[i]);
}
if (b_try_scalable && is_available(f_cvc_scalable)) {
enable(f_cvc_scalable);
}
return COLVARS_OK;
}
colvar::cvc::~cvc()
{}
{
remove_all_children();
for (size_t i = 0; i < atom_groups.size(); i++) {
if (atom_groups[i] != NULL) delete atom_groups[i];
}
}
void colvar::cvc::read_data()
{
size_t ig;
for (ig = 0; ig < atom_groups.size(); ig++) {
cvm::atom_group &atoms = *(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();
// }
// }
// }
}
void colvar::cvc::calc_force_invgrads()
@ -86,15 +153,15 @@ void colvar::cvc::calc_Jacobian_derivative()
}
void colvar::cvc::debug_gradients(cvm::atom_group &group)
void colvar::cvc::debug_gradients(cvm::atom_group *group)
{
// this function should work for any scalar variable:
// the only difference will be the name of the atom group (here, "group")
if (group.b_dummy) return;
if (group->b_dummy) return;
cvm::rotation const rot_0 = group.rot;
cvm::rotation const rot_inv = group.rot.inverse();
cvm::rotation const rot_0 = group->rot;
cvm::rotation const rot_inv = group->rot.inverse();
cvm::real x_0 = x.real_value;
if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0];
@ -103,45 +170,45 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group)
// it only makes sense to debug the fit gradients
// when the fitting group is the same as this group
if (group.b_rotate || group.b_center)
if (group.b_fit_gradients && (group.ref_pos_group == NULL)) {
group.calc_fit_gradients();
if (group.b_rotate) {
if (group->b_rotate || group->b_center)
if (group->b_fit_gradients && (group->ref_pos_group == NULL)) {
group->calc_fit_gradients();
if (group->b_rotate) {
// fit_gradients are in the original frame, we should print them in the rotated frame
for (size_t j = 0; j < group.fit_gradients.size(); j++) {
group.fit_gradients[j] = rot_0.rotate(group.fit_gradients[j]);
for (size_t j = 0; j < group->fit_gradients.size(); j++) {
group->fit_gradients[j] = rot_0.rotate(group->fit_gradients[j]);
}
}
cvm::log("fit_gradients = "+cvm::to_str(group.fit_gradients)+"\n");
if (group.b_rotate) {
for (size_t j = 0; j < group.fit_gradients.size(); j++) {
group.fit_gradients[j] = rot_inv.rotate(group.fit_gradients[j]);
cvm::log("fit_gradients = "+cvm::to_str(group->fit_gradients)+"\n");
if (group->b_rotate) {
for (size_t j = 0; j < group->fit_gradients.size(); j++) {
group->fit_gradients[j] = rot_inv.rotate(group->fit_gradients[j]);
}
}
}
for (size_t ia = 0; ia < group.size(); ia++) {
for (size_t ia = 0; ia < group->size(); ia++) {
// tests are best conducted in the unrotated (simulation) frame
cvm::rvector const atom_grad = group.b_rotate ?
rot_inv.rotate(group[ia].grad) :
group[ia].grad;
cvm::rvector const atom_grad = group->b_rotate ?
rot_inv.rotate((*group)[ia].grad) :
(*group)[ia].grad;
for (size_t id = 0; id < 3; id++) {
// (re)read original positions
group.read_positions();
group->read_positions();
// change one coordinate
group[ia].pos[id] += cvm::debug_gradients_step_size;
group.calc_required_properties();
(*group)[ia].pos[id] += cvm::debug_gradients_step_size;
group->calc_required_properties();
calc_value();
cvm::real x_1 = x.real_value;
if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n");
cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0,
21, 14)+"\n");
//cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ?
cvm::real const dx_pred = (group.fit_gradients.size()) ?
(cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) :
//cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
cvm::real const dx_pred = (group->fit_gradients.size()) ?
(cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
(cvm::debug_gradients_step_size * atom_grad[id]);
cvm::log("dx(interp) = "+cvm::to_str(dx_pred,
21, 14)+"\n");
@ -154,27 +221,27 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group)
/*
* The code below is WIP
*/
// if (group.ref_pos_group != NULL) {
// cvm::atom_group &ref = *group.ref_pos_group;
// group.calc_fit_gradients();
// if (group->ref_pos_group != NULL) {
// cvm::atom_group &ref = *group->ref_pos_group;
// group->calc_fit_gradients();
//
// for (size_t ia = 0; ia < ref.size(); ia++) {
//
// for (size_t id = 0; id < 3; id++) {
// // (re)read original positions
// group.read_positions();
// group->read_positions();
// ref.read_positions();
// // change one coordinate
// ref[ia].pos[id] += cvm::debug_gradients_step_size;
// group.update();
// group->update();
// calc_value();
// cvm::real const x_1 = x.real_value;
// cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");
// cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0,
// 21, 14)+"\n");
// //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ?
// // cvm::real const dx_pred = (group.fit_gradients.size()) ?
// // (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) :
// //cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
// // cvm::real const dx_pred = (group->fit_gradients.size()) ?
// // (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
// // (cvm::debug_gradients_step_size * atom_grad[id]);
// cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][id];
// cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
@ -190,3 +257,8 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group)
return;
}
// Static members
std::vector<cvm::deps::feature *> colvar::cvc::cvc_features;

83
lib/colvars/colvarcomp.h
View File

@ -61,7 +61,7 @@
/// call to e.g. apply_force().
class colvar::cvc
: public colvarparse
: public colvarparse, public cvm::deps
{
public:
@ -84,10 +84,6 @@ public:
/// \brief Exponent in the polynomial combination (default: 1)
int sup_np;
/// \brief This defaults to true; setting it to false disables
/// update of this cvc to save compute time (useful with scriptedFunction)
bool b_enabled;
/// \brief Is this a periodic component?
bool b_periodic;
@ -107,11 +103,14 @@ public:
/// \brief Within the constructor, make a group parse its own
/// options from the provided configuration string
void parse_group(std::string const &conf,
/// Returns reference to new group
cvm::atom_group *parse_group(std::string const &conf,
char const *group_key,
cvm::atom_group &group,
bool optional = false);
/// \brief After construction, set data related to dependency handling
int setup();
/// \brief Default constructor (used when \link cvc \endlink
/// objects are declared within other ones)
cvc();
@ -119,18 +118,16 @@ public:
/// Destructor
virtual ~cvc();
/// \brief Implementation of the feature list for colvar
static std::vector<feature *> cvc_features;
/// \brief If this flag is false (default), inverse gradients
/// (derivatives of atom coordinates with respect to x) are
/// unavailable; it should be set to true by the constructor of each
/// derived object capable of calculating them
bool b_inverse_gradients;
/// \brief Implementation of the feature list accessor for colvar
virtual std::vector<feature *> &features() {
return cvc_features;
}
/// \brief If this flag is false (default), the Jacobian derivative
/// (divergence of the inverse gradients) is unavailable; it should
/// be set to true by the constructor of each derived object capable
/// of calculating it
bool b_Jacobian_derivative;
/// \brief Obtain data needed for the calculation for the backend
void read_data();
/// \brief Calculate the variable
virtual void calc_value() = 0;
@ -139,11 +136,8 @@ public:
/// order to apply forces
virtual void calc_gradients() = 0;
/// \brief If true, calc_gradients() will call debug_gradients() for every group needed
bool b_debug_gradients;
/// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
virtual void debug_gradients(cvm::atom_group &group);
virtual void debug_gradients(cvm::atom_group *group);
/// \brief Calculate the total force from the system using the
/// inverse atomic gradients
@ -227,6 +221,9 @@ public:
/// e.g. atomic gradients
std::vector<cvm::atom_group *> atom_groups;
/// \brief Whether or not this CVC will be computed in parallel whenever possible
bool b_try_scalable;
protected:
/// \brief Cached value
@ -307,9 +304,9 @@ class colvar::distance
{
protected:
/// First atom group
cvm::atom_group group1;
cvm::atom_group *group1;
/// Second atom group
cvm::atom_group group2;
cvm::atom_group *group2;
/// Vector distance, cached to be recycled
cvm::rvector dist_v;
/// Use absolute positions, ignoring PBCs when present
@ -389,11 +386,11 @@ class colvar::distance_z
{
protected:
/// Main atom group
cvm::atom_group main;
cvm::atom_group *main;
/// Reference atom group
cvm::atom_group ref1;
cvm::atom_group *ref1;
/// Optional, second ref atom group
cvm::atom_group ref2;
cvm::atom_group *ref2;
/// Use absolute positions, ignoring PBCs when present
bool b_no_PBC;
/// Compute system force on one site only to avoid unwanted
@ -486,9 +483,9 @@ class colvar::distance_pairs
{
protected:
/// First atom group
cvm::atom_group group1;
cvm::atom_group *group1;
/// Second atom group
cvm::atom_group group2;
cvm::atom_group *group2;
/// Use absolute positions, ignoring PBCs when present
bool b_no_PBC;
public:
@ -515,7 +512,7 @@ class colvar::gyration
{
protected:
/// Atoms involved
cvm::atom_group atoms;
cvm::atom_group *atoms;
public:
/// Constructor
gyration(std::string const &conf);
@ -590,7 +587,7 @@ class colvar::eigenvector
protected:
/// Atom group
cvm::atom_group atoms;
cvm::atom_group * atoms;
/// Reference coordinates
std::vector<cvm::atom_pos> ref_pos;
@ -632,11 +629,11 @@ class colvar::angle
protected:
/// Atom group
cvm::atom_group group1;
cvm::atom_group *group1;
/// Atom group
cvm::atom_group group2;
cvm::atom_group *group2;
/// Atom group
cvm::atom_group group3;
cvm::atom_group *group3;
/// Inter site vectors
cvm::rvector r21, r23;
@ -678,11 +675,11 @@ class colvar::dipole_angle
protected:
/// Dipole atom group
cvm::atom_group group1;
cvm::atom_group *group1;
/// Atom group
cvm::atom_group group2;
cvm::atom_group *group2;
/// Atom group
cvm::atom_group group3;
cvm::atom_group *group3;
/// Inter site vectors
cvm::rvector r21, r23;
@ -722,13 +719,13 @@ class colvar::dihedral
protected:
/// Atom group
cvm::atom_group group1;
cvm::atom_group *group1;
/// Atom group
cvm::atom_group group2;
cvm::atom_group *group2;
/// Atom group
cvm::atom_group group3;
cvm::atom_group *group3;
/// Atom group
cvm::atom_group group4;
cvm::atom_group *group4;
/// Inter site vectors
cvm::rvector r12, r23, r34;
@ -1009,7 +1006,7 @@ class colvar::orientation
protected:
/// Atom group
cvm::atom_group atoms;
cvm::atom_group * atoms;
/// Center of geometry of the group
cvm::atom_pos atoms_cog;
@ -1155,7 +1152,7 @@ class colvar::rmsd
protected:
/// Atom group
cvm::atom_group atoms;
cvm::atom_group *atoms;
/// Reference coordinates (for RMSD calculation only)
std::vector<cvm::atom_pos> ref_pos;
@ -1186,7 +1183,7 @@ class colvar::cartesian
{
protected:
/// Atom group
cvm::atom_group atoms;
cvm::atom_group *atoms;
/// Which Cartesian coordinates to include
std::vector<size_t> axes;
public:

249
lib/colvars/colvarcomp_angles.cpp
View File

@ -11,14 +11,13 @@ colvar::angle::angle(std::string const &conf)
: cvc(conf)
{
function_type = "angle";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
parse_group(conf, "group1", group1);
parse_group(conf, "group2", group2);
parse_group(conf, "group3", group3);
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
provide(f_cvc_com_based);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
group3 = parse_group(conf, "group3");
if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
cvm::log("Computing system force on group 1 only");
}
@ -27,19 +26,21 @@ colvar::angle::angle(std::string const &conf)
colvar::angle::angle(cvm::atom const &a1,
cvm::atom const &a2,
cvm::atom const &a3)
: group1(std::vector<cvm::atom>(1, a1)),
group2(std::vector<cvm::atom>(1, a2)),
group3(std::vector<cvm::atom>(1, a3))
cvm::atom const &a2,
cvm::atom const &a3)
{
function_type = "angle";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
provide(f_cvc_com_based);
b_1site_force = false;
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
x.type(colvarvalue::type_scalar);
}
@ -54,9 +55,9 @@ colvar::angle::angle()
void colvar::angle::calc_value()
{
cvm::atom_pos const g1_pos = group1.center_of_mass();
cvm::atom_pos const g2_pos = group2.center_of_mass();
cvm::atom_pos const g3_pos = group3.center_of_mass();
cvm::atom_pos const g1_pos = group1->center_of_mass();
cvm::atom_pos const g2_pos = group2->center_of_mass();
cvm::atom_pos const g3_pos = group3->center_of_mass();
r21 = cvm::position_distance(g2_pos, g1_pos);
r21l = r21.norm();
@ -71,7 +72,6 @@ void colvar::angle::calc_value()
void colvar::angle::calc_gradients()
{
size_t i;
cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
@ -81,22 +81,11 @@ void colvar::angle::calc_gradients()
dxdr3 = (180.0/PI) * dxdcos *
(1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
for (i = 0; i < group1.size(); i++) {
group1[i].grad = (group1[i].mass/group1.total_mass) *
(dxdr1);
}
group1->set_weighted_gradient(dxdr1);
group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
group3->set_weighted_gradient(dxdr3);
for (i = 0; i < group2.size(); i++) {
group2[i].grad = (group2[i].mass/group2.total_mass) *
(dxdr1 + dxdr3) * (-1.0);
}
for (i = 0; i < group3.size(); i++) {
group3[i].grad = (group3[i].mass/group3.total_mass) *
(dxdr3);
}
if (b_debug_gradients) {
if (is_enabled(f_cvc_debug_gradient)) {
debug_gradients(group1);
debug_gradients(group2);
debug_gradients(group3);
@ -112,15 +101,15 @@ void colvar::angle::calc_force_invgrads()
// when propagating changes in the angle)
if (b_1site_force) {
group1.read_system_forces();
group1->read_system_forces();
cvm::real norm_fact = 1.0 / dxdr1.norm2();
ft.real_value = norm_fact * dxdr1 * group1.system_force();
ft.real_value = norm_fact * dxdr1 * group1->system_force();
} else {
group1.read_system_forces();
group3.read_system_forces();
group1->read_system_forces();
group3->read_system_forces();
cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());
ft.real_value = norm_fact * ( dxdr1 * group1.system_force()
+ dxdr3 * group3.system_force());
ft.real_value = norm_fact * ( dxdr1 * group1->system_force()
+ dxdr3 * group3->system_force());
}
return;
}
@ -136,14 +125,14 @@ void colvar::angle::calc_Jacobian_derivative()
void colvar::angle::apply_force(colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force(force.real_value);
if (!group1->noforce)
group1->apply_colvar_force(force.real_value);
if (!group2.noforce)
group2.apply_colvar_force(force.real_value);
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
if (!group3.noforce)
group3.apply_colvar_force(force.real_value);
if (!group3->noforce)
group3->apply_colvar_force(force.real_value);
}
@ -153,13 +142,10 @@ colvar::dipole_angle::dipole_angle(std::string const &conf)
: cvc(conf)
{
function_type = "dipole_angle";
parse_group(conf, "group1", group1);
parse_group(conf, "group2", group2);
parse_group(conf, "group3", group3);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
group3 = parse_group(conf, "group3");
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
cvm::log("Computing system force on group 1 only");
}
@ -170,15 +156,16 @@ colvar::dipole_angle::dipole_angle(std::string const &conf)
colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
cvm::atom const &a2,
cvm::atom const &a3)
: group1(std::vector<cvm::atom>(1, a1)),
group2(std::vector<cvm::atom>(1, a2)),
group3(std::vector<cvm::atom>(1, a3))
{
function_type = "dipole_angle";
b_1site_force = false;
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
x.type(colvarvalue::type_scalar);
}
@ -193,13 +180,13 @@ colvar::dipole_angle::dipole_angle()
void colvar::dipole_angle::calc_value()
{
cvm::atom_pos const g1_pos = group1.center_of_mass();
cvm::atom_pos const g2_pos = group2.center_of_mass();
cvm::atom_pos const g3_pos = group3.center_of_mass();
cvm::atom_pos const g1_pos = group1->center_of_mass();
cvm::atom_pos const g2_pos = group2->center_of_mass();
cvm::atom_pos const g3_pos = group3->center_of_mass();
group1.calc_dipole(g1_pos);
group1->calc_dipole(g1_pos);
r21 = group1.dipole();
r21 = group1->dipole();
r21l = r21.norm();
r23 = cvm::position_distance(g2_pos, g3_pos);
r23l = r23.norm();
@ -225,35 +212,35 @@ void colvar::dipole_angle::calc_gradients()
(1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
//this auxiliar variables are to avoid numerical errors inside "for"
double aux1 = group1.total_charge/group1.total_mass;
// double aux2 = group2.total_charge/group2.total_mass;
// double aux3 = group3.total_charge/group3.total_mass;
double aux1 = group1->total_charge/group1->total_mass;
// double aux2 = group2->total_charge/group2->total_mass;
// double aux3 = group3->total_charge/group3->total_mass;
size_t i;
for (i = 0; i < group1.size(); i++) {
group1[i].grad =(group1[i].charge + (-1)* group1[i].mass * aux1) * (dxdr1);
for (i = 0; i < group1->size(); i++) {
(*group1)[i].grad =((*group1)[i].charge + (-1)* (*group1)[i].mass * aux1) * (dxdr1);
}
for (i = 0; i < group2.size(); i++) {
group2[i].grad = (group2[i].mass/group2.total_mass)* dxdr3 * (-1.0);
for (i = 0; i < group2->size(); i++) {
(*group2)[i].grad = ((*group2)[i].mass/group2->total_mass)* dxdr3 * (-1.0);
}
for (i = 0; i < group3.size(); i++) {
group3[i].grad =(group3[i].mass/group3.total_mass) * (dxdr3);
for (i = 0; i < group3->size(); i++) {
(*group3)[i].grad =((*group3)[i].mass/group3->total_mass) * (dxdr3);
}
}
void colvar::dipole_angle::apply_force(colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force(force.real_value);
if (!group1->noforce)
group1->apply_colvar_force(force.real_value);
if (!group2.noforce)
group2.apply_colvar_force(force.real_value);
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
if (!group3.noforce)
group3.apply_colvar_force(force.real_value);
if (!group3->noforce)
group3->apply_colvar_force(force.real_value);
}
@ -265,32 +252,26 @@ colvar::dihedral::dihedral(std::string const &conf)
function_type = "dihedral";
period = 360.0;
b_periodic = true;
b_inverse_gradients = true;
b_Jacobian_derivative = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
provide(f_cvc_com_based);
if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
cvm::log("Computing system force on group 1 only");
}
parse_group(conf, "group1", group1);
parse_group(conf, "group2", group2);
parse_group(conf, "group3", group3);
parse_group(conf, "group4", group4);
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
atom_groups.push_back(&group4);
group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2");
group3 = parse_group(conf, "group3");
group4 = parse_group(conf, "group4");
x.type(colvarvalue::type_scalar);
}
colvar::dihedral::dihedral(cvm::atom const &a1,
cvm::atom const &a2,
cvm::atom const &a3,
cvm::atom const &a4)
: group1(std::vector<cvm::atom>(1, a1)),
group2(std::vector<cvm::atom>(1, a2)),
group3(std::vector<cvm::atom>(1, a3)),
group4(std::vector<cvm::atom>(1, a4))
cvm::atom const &a2,
cvm::atom const &a3,
cvm::atom const &a4)
{
if (cvm::debug())
cvm::log("Initializing dihedral object from atom groups.\n");
@ -298,14 +279,20 @@ colvar::dihedral::dihedral(cvm::atom const &a1,
function_type = "dihedral";
period = 360.0;
b_periodic = true;
b_inverse_gradients = true;
b_Jacobian_derivative = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
provide(f_cvc_com_based);
b_1site_force = false;
atom_groups.push_back(&group1);
atom_groups.push_back(&group2);
atom_groups.push_back(&group3);
atom_groups.push_back(&group4);
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
group3 = new cvm::atom_group(std::vector<cvm::atom>(1, a3));
group4 = new cvm::atom_group(std::vector<cvm::atom>(1, a4));
atom_groups.push_back(group1);
atom_groups.push_back(group2);
atom_groups.push_back(group3);
atom_groups.push_back(group4);
x.type(colvarvalue::type_scalar);
@ -319,18 +306,18 @@ colvar::dihedral::dihedral()
function_type = "dihedral";
period = 360.0;
b_periodic = true;
b_inverse_gradients = true;
b_Jacobian_derivative = true;
provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian);
x.type(colvarvalue::type_scalar);
}
void colvar::dihedral::calc_value()
{
cvm::atom_pos const g1_pos = group1.center_of_mass();
cvm::atom_pos const g2_pos = group2.center_of_mass();
cvm::atom_pos const g3_pos = group3.center_of_mass();
cvm::atom_pos const g4_pos = group4.center_of_mass();
cvm::atom_pos const g1_pos = group1->center_of_mass();
cvm::atom_pos const g2_pos = group2->center_of_mass();
cvm::atom_pos const g3_pos = group3->center_of_mass();
cvm::atom_pos const g4_pos = group4->center_of_mass();
// Usual sign convention: r12 = r2 - r1
r12 = cvm::position_distance(g1_pos, g2_pos);
@ -415,18 +402,10 @@ void colvar::dihedral::calc_gradients()
+dsindB.y*r34.x - dsindB.x*r34.y);
}
size_t i;
for (i = 0; i < group1.size(); i++)
group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1);
for (i = 0; i < group2.size(); i++)
group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1);
for (i = 0; i < group3.size(); i++)
group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2);
for (i = 0; i < group4.size(); i++)
group4[i].grad = (group4[i].mass/group4.total_mass) * (f3);
group1->set_weighted_gradient(-f1);
group2->set_weighted_gradient(-f2 + f1);
group3->set_weighted_gradient(-f3 + f2);
group4->set_weighted_gradient(f3);
}
@ -448,15 +427,15 @@ void colvar::dihedral::calc_force_invgrads()
cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1);
cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);
group1.read_system_forces();
group1->read_system_forces();
if ( b_1site_force ) {
// This is only measuring the force on group 1
ft.real_value = PI/180.0 * fact1 * (cross1 * group1.system_force());
ft.real_value = PI/180.0 * fact1 * (cross1 * group1->system_force());
} else {
// Default case: use groups 1 and 4
group4.read_system_forces();
ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1.system_force())
+ fact4 * (cross4 * group4.system_force()));
group4->read_system_forces();
ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1->system_force())
+ fact4 * (cross4 * group4->system_force()));
}
}
@ -470,17 +449,17 @@ void colvar::dihedral::calc_Jacobian_derivative()
void colvar::dihedral::apply_force(colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force(force.real_value);
if (!group1->noforce)
group1->apply_colvar_force(force.real_value);
if (!group2.noforce)
group2.apply_colvar_force(force.real_value);
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
if (!group3.noforce)
group3.apply_colvar_force(force.real_value);
if (!group3->noforce)
group3->apply_colvar_force(force.real_value);
if (!group4.noforce)
group4.apply_colvar_force(force.real_value);
if (!group4->noforce)
group4->apply_colvar_force(force.real_value);
}

63
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -79,13 +79,13 @@ colvar::coordnum::coordnum(std::string const &conf)
x.type(colvarvalue::type_scalar);
// group1 and group2 are already initialized by distance()
if (group1.b_dummy)
if (group1->b_dummy)
cvm::fatal_error("Error: only group2 is allowed to be a dummy atom\n");
// need to specify this explicitly because the distance() constructor
// has set it to true
b_inverse_gradients = false;
feature_states[f_cvc_inv_gradient]->available = false;
bool const b_scale = get_keyval(conf, "cutoff", r0,
cvm::real(4.0 * cvm::unit_angstrom()));
@ -110,7 +110,7 @@ colvar::coordnum::coordnum(std::string const &conf)
cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
}
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy);
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
}
@ -130,26 +130,26 @@ void colvar::coordnum::calc_value()
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2.center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
}
} else {
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
}
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
}
}
@ -163,29 +163,29 @@ void colvar::coordnum::calc_gradients()
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2.center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
}
group2.set_weighted_gradient(group2_com_atom.grad);
group2->set_weighted_gradient(group2_com_atom.grad);
} else {
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
}
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
switching_function<true>(r0, en, ed, *ai1, *ai2);
}
}
@ -194,11 +194,11 @@ void colvar::coordnum::calc_gradients()
void colvar::coordnum::apply_force(colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force(force.real_value);
if (!group1->noforce)
group1->apply_colvar_force(force.real_value);
if (!group2.noforce)
group2.apply_colvar_force(force.real_value);
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
}
@ -298,7 +298,8 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
// need to specify this explicitly because the distance() constructor
// has set it to true
b_inverse_gradients = false;
feature_states[f_cvc_inv_gradient]->available = false;
get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
get_keyval(conf, "expNumer", en, int(6) );
@ -320,9 +321,9 @@ colvar::selfcoordnum::selfcoordnum()
void colvar::selfcoordnum::calc_value()
{
x.real_value = 0.0;
for (size_t i = 0; i < group1.size() - 1; i++) {
for (size_t j = i + 1; j < group1.size(); j++) {
x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, group1[i], group1[j]);
for (size_t i = 0; i < group1->size() - 1; i++) {
for (size_t j = i + 1; j < group1->size(); j++) {
x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, (*group1)[i], (*group1)[j]);
}
}
}
@ -330,17 +331,17 @@ void colvar::selfcoordnum::calc_value()
void colvar::selfcoordnum::calc_gradients()
{
for (size_t i = 0; i < group1.size() - 1; i++) {
for (size_t j = i + 1; j < group1.size(); j++) {
colvar::coordnum::switching_function<true>(r0, en, ed, group1[i], group1[j]);
for (size_t i = 0; i < group1->size() - 1; i++) {
for (size_t j = i + 1; j < group1->size(); j++) {
colvar::coordnum::switching_function<true>(r0, en, ed, (*group1)[i], (*group1)[j]);
}
}
}
void colvar::selfcoordnum::apply_force(colvarvalue const &force)
{
if (!group1.noforce) {
group1.apply_colvar_force(force.real_value);
if (!group1->noforce) {
group1->apply_colvar_force(force.real_value);
}
}

587
lib/colvars/colvarcomp_distances.cpp
View File

File diff suppressed because it is too large Load Diff

85
lib/colvars/colvarcomp_rotations.cpp
View File

@ -14,15 +14,14 @@ colvar::orientation::orientation(std::string const &conf)
: cvc(conf)
{
function_type = "orientation";
parse_group(conf, "atoms", atoms);
atom_groups.push_back(&atoms);
atoms = parse_group(conf, "atoms");
x.type(colvarvalue::type_quaternion);
ref_pos.reserve(atoms.size());
ref_pos.reserve(atoms->size());
if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
cvm::log("Using reference positions from input file.\n");
if (ref_pos.size() != atoms.size()) {
if (ref_pos.size() != atoms->size()) {
cvm::fatal_error("Error: reference positions do not "
"match the number of requested atoms.\n");
}
@ -42,10 +41,10 @@ colvar::orientation::orientation(std::string const &conf)
"if provided, must be non-zero.\n");
} else {
// if not, use atom indices
atoms.create_sorted_ids();
atoms->create_sorted_ids();
}
ref_pos.resize(atoms.size());
cvm::load_coords(file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
ref_pos.resize(atoms->size());
cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
}
}
@ -71,8 +70,8 @@ colvar::orientation::orientation(std::string const &conf)
get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0));
// initialize rot member data
if (!atoms.noforce) {
rot.request_group2_gradients(atoms.size());
if (!atoms->noforce) {
rot.request_group2_gradients(atoms->size());
}
}
@ -88,9 +87,9 @@ colvar::orientation::orientation()
void colvar::orientation::calc_value()
{
atoms_cog = atoms.center_of_geometry();
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
if ((rot.q).inner(ref_quat) >= 0.0) {
x.quaternion_value = rot.q;
@ -113,10 +112,10 @@ void colvar::orientation::apply_force(colvarvalue const &force)
{
cvm::quaternion const &FQ = force.quaternion_value;
if (!atoms.noforce) {
for (size_t ia = 0; ia < atoms.size(); ia++) {
if (!atoms->noforce) {
for (size_t ia = 0; ia < atoms->size(); ia++) {
for (size_t i = 0; i < 4; i++) {
atoms[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]);
(*atoms)[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]);
}
}
}
@ -142,9 +141,9 @@ colvar::orientation_angle::orientation_angle()
void colvar::orientation_angle::calc_value()
{
atoms_cog = atoms.center_of_geometry();
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
if ((rot.q).q0 >= 0.0) {
x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0);
@ -161,10 +160,10 @@ void colvar::orientation_angle::calc_gradients()
((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) :
0.0 );
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
}
if (b_debug_gradients) {
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging orientationAngle component gradients:\n");
debug_gradients(atoms);
}
@ -174,8 +173,8 @@ void colvar::orientation_angle::calc_gradients()
void colvar::orientation_angle::apply_force(colvarvalue const &force)
{
cvm::real const &fw = force.real_value;
if (!atoms.noforce) {
atoms.apply_colvar_force(fw);
if (!atoms->noforce) {
atoms->apply_colvar_force(fw);
}
}
@ -199,8 +198,8 @@ colvar::orientation_proj::orientation_proj()
void colvar::orientation_proj::calc_value()
{
atoms_cog = atoms.center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0;
}
@ -208,10 +207,10 @@ void colvar::orientation_proj::calc_value()
void colvar::orientation_proj::calc_gradients()
{
cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
}
if (b_debug_gradients) {
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging orientationProj component gradients:\n");
debug_gradients(atoms);
}
@ -222,8 +221,8 @@ void colvar::orientation_proj::apply_force(colvarvalue const &force)
{
cvm::real const &fw = force.real_value;
if (!atoms.noforce) {
atoms.apply_colvar_force(fw);
if (!atoms->noforce) {
atoms->apply_colvar_force(fw);
}
}
@ -255,9 +254,9 @@ colvar::tilt::tilt()
void colvar::tilt::calc_value()
{
atoms_cog = atoms.center_of_geometry();
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
x.real_value = rot.cos_theta(axis);
}
@ -267,14 +266,14 @@ void colvar::tilt::calc_gradients()
{
cvm::quaternion const dxdq = rot.dcos_theta_dq(axis);
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
for (size_t iq = 0; iq < 4; iq++) {
atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
(*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
}
}
if (b_debug_gradients) {
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging tilt component gradients:\n");
debug_gradients(atoms);
}
@ -285,8 +284,8 @@ void colvar::tilt::apply_force(colvarvalue const &force)
{
cvm::real const &fw = force.real_value;
if (!atoms.noforce) {
atoms.apply_colvar_force(fw);
if (!atoms->noforce) {
atoms->apply_colvar_force(fw);
}
}
@ -322,9 +321,9 @@ colvar::spin_angle::spin_angle()
void colvar::spin_angle::calc_value()
{
atoms_cog = atoms.center_of_geometry();
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
x.real_value = rot.spin_angle(axis);
this->wrap(x);
@ -335,10 +334,10 @@ void colvar::spin_angle::calc_gradients()
{
cvm::quaternion const dxdq = rot.dspin_angle_dq(axis);
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
for (size_t iq = 0; iq < 4; iq++) {
atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
(*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
}
}
}
@ -348,7 +347,7 @@ void colvar::spin_angle::apply_force(colvarvalue const &force)
{
cvm::real const &fw = force.real_value;
if (!atoms.noforce) {
atoms.apply_colvar_force(fw);
if (!atoms->noforce) {
atoms->apply_colvar_force(fw);
}
}

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

9
lib/colvars/colvargrid.h
View File

@ -251,15 +251,6 @@ public:
size_t i;
for (i = 0; i < cv.size(); i++) {
if (!cv[i]->tasks[colvar::task_lower_boundary] ||
!cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::error("Tried to initialize a grid on a "
"variable with undefined boundaries.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
}
if (cvm::debug()) {
cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+
" collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");

409
lib/colvars/colvarmodule.cpp
View File

@ -15,6 +15,7 @@
#include "colvarbias_restraint.h"
#include "colvarscript.h"
colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
// pointer to the proxy object
@ -100,6 +101,7 @@ int colvarmodule::read_config_string(std::string const &config_str)
return parse_config(conf);
}
int colvarmodule::parse_config(std::string &conf)
{
// parse global options
@ -165,8 +167,9 @@ int colvarmodule::parse_global_params(std::string const &conf)
// 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);
parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);
parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true);
if (use_scripted_forces && !proxy->force_script_defined) {
cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
@ -218,6 +221,7 @@ int colvarmodule::parse_colvars(std::string const &conf)
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
bool colvarmodule::check_new_bias(std::string &conf, char const *key)
{
if (cvm::get_error() ||
@ -289,6 +293,12 @@ int colvarmodule::parse_biases(std::string const &conf)
cvm::decrease_depth();
}
for (size_t i = 0; i < biases.size(); i++) {
biases[i]->enable(cvm::deps::f_cvb_active);
if (cvm::debug())
biases[i]->print_state();
}
if (biases.size() || use_scripted_forces) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables biases initialized, "+
@ -448,12 +458,9 @@ int colvarmodule::bias_share(std::string const &bias_name)
}
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;
int colvarmodule::calc()
{
int error_code = COLVARS_OK;
if (cvm::debug()) {
cvm::log(cvm::line_marker);
@ -461,75 +468,87 @@ int colvarmodule::calc() {
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++) {
error_code |= calc_colvars();
// set biasing forces to zero before biases are calculated and summed over
for (std::vector<colvar *>::iterator 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();
error_code |= calc_biases();
error_code |= update_colvar_forces();
if (cvm::b_analysis) {
// perform runtime analysis of colvars and biases
if (cvm::debug() && biases.size())
cvm::log("Perform runtime analyses.\n");
error_code |= analyze();
}
// write trajectory files, if needed
if (cv_traj_freq && cv_traj_name.size()) {
error_code |= write_traj_files();
}
// write restart files, if needed
if (restart_out_freq && restart_out_name.size()) {
error_code |= write_restart_files();
}
return error_code;
}
int colvarmodule::calc_colvars()
{
if (cvm::debug())
cvm::log("Calculating collective variables.\n");
// calculate collective variables and their gradients
int error_code = COLVARS_OK;
std::vector<colvar *>::iterator cvi;
// if SMP support is available, split up the work
if (proxy->smp_enabled() == COLVARS_OK) {
// first, calculate how much work (currently, how many active CVCs) each colvar has
colvars_smp.resize(0);
colvars_smp_items.resize(0);
colvars_smp.reserve(colvars.size());
colvars_smp_items.reserve(colvars.size());
// set up a vector containing all components
size_t num_colvar_items = 0;
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->analyze();
if (cvm::get_error()) {
return COLVARS_ERROR;
error_code |= (*cvi)->update_cvc_flags();
size_t num_items = (*cvi)->num_active_cvcs();
colvars_smp.reserve(colvars_smp.size() + num_items);
colvars_smp_items.reserve(colvars_smp_items.size() + num_items);
for (size_t icvc = 0; icvc < num_items; icvc++) {
colvars_smp.push_back(*cvi);
colvars_smp_items.push_back(icvc);
}
num_colvar_items += num_items;
}
for (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->analyze();
cvm::decrease_depth();
// calculate colvar components in parallel
error_code |= proxy->smp_colvars_loop();
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
error_code |= (*cvi)->collect_cvc_data();
}
cvm::decrease_depth();
} else {
// calculate colvars one at a time
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
error_code |= (*cvi)->calc();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
@ -537,6 +556,81 @@ int colvarmodule::calc() {
cvm::decrease_depth();
}
return error_code | cvm::get_error();
}
int colvarmodule::calc_biases()
{
// update the biases and communicate their forces to the collective
// variables
if (cvm::debug() && biases.size())
cvm::log("Updating collective variable biases.\n");
std::vector<colvarbias *>::iterator bi;
int error_code = COLVARS_OK;
// if SMP support is available, split up the work
if (proxy->smp_enabled() == COLVARS_OK) {
if (use_scripted_forces && !scripting_after_biases) {
// calculate biases and scripted forces in parallel
error_code |= proxy->smp_biases_script_loop();
} else {
// calculate biases in parallel
error_code |= proxy->smp_biases_loop();
}
} else {
if (use_scripted_forces && !scripting_after_biases) {
error_code |= calc_scripted_forces();
}
cvm::increase_depth();
for (bi = biases.begin(); bi != biases.end(); bi++) {
error_code |= (*bi)->update();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
}
cvm::real total_bias_energy = 0.0;
for (bi = biases.begin(); bi != biases.end(); bi++) {
total_bias_energy += (*bi)->get_energy();
}
proxy->add_energy(total_bias_energy);
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::update_colvar_forces()
{
int error_code = COLVARS_OK;
std::vector<colvar *>::iterator cvi;
std::vector<colvarbias *>::iterator bi;
// 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 (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
if (use_scripted_forces && scripting_after_biases) {
error_code |= calc_scripted_forces();
}
cvm::real total_colvar_energy = 0.0;
// sum up the forces for each colvar, including wall forces
// and integrate any internal
// equation of motion (extended system)
@ -545,13 +639,13 @@ int colvarmodule::calc() {
"of colvars (if they have any).\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
total_colvar_energy += (*cvi)->update();
total_colvar_energy += (*cvi)->update_forces_energy();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
proxy->add_energy(total_bias_energy + total_colvar_energy);
proxy->add_energy(total_colvar_energy);
// make collective variables communicate their forces to their
// coupled degrees of freedom (i.e. atoms)
@ -559,7 +653,7 @@ int colvarmodule::calc() {
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]) {
if ((*cvi)->is_enabled(cvm::deps::f_cv_gradient)) {
(*cvi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
@ -568,46 +662,68 @@ int colvarmodule::calc() {
}
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();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::calc_scripted_forces()
{
// Run user force script, if provided,
// potentially adding scripted forces to the colvars
int res;
res = proxy->run_force_callback();
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Colvar forces scripts are not implemented.");
return COLVARS_NOT_IMPLEMENTED;
}
if (res != COLVARS_OK) {
cvm::error("Error running user colvar forces script");
return COLVARS_ERROR;
}
}
int colvarmodule::write_restart_files()
{
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()) {
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
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);
}
int colvarmodule::write_traj_files()
{
if (!cv_traj_os.is_open()) {
open_traj_file(cv_traj_name);
}
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
write_traj(cv_traj_os);
}
// 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 (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();
}
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);
}
@ -643,6 +759,7 @@ int colvarmodule::analyze()
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::setup()
{
// loop over all components of all colvars to reset masses of all groups
@ -653,25 +770,23 @@ int colvarmodule::setup()
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
colvarmodule::~colvarmodule()
{
reset();
delete parse;
proxy = NULL;
if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) ||
(proxy->smp_thread_id() == 0)) {
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();
@ -681,6 +796,14 @@ int colvarmodule::reset()
}
biases.clear();
// 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();
index_groups.clear();
index_group_names.clear();
@ -1048,28 +1171,71 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
void cvm::log(std::string const &message)
{
if (depth > 0)
proxy->log((std::string(2*depth, ' '))+message);
size_t const d = depth();
if (d > 0)
proxy->log((std::string(2*d, ' '))+message);
else
proxy->log(message);
}
void cvm::increase_depth()
{
depth++;
(depth())++;
}
void cvm::decrease_depth()
{
if (depth) depth--;
if (depth() > 0) {
(depth())--;
}
}
size_t & cvm::depth()
{
// NOTE: do not call log() or error() here, to avoid recursion
size_t const nt = proxy->smp_num_threads();
if (proxy->smp_enabled() == COLVARS_OK) {
if (depth_v.size() != nt) {
// update array of depths
proxy->smp_lock();
if (depth_v.size() > 0) { depth_s = depth_v[0]; }
depth_v.clear();
depth_v.assign(nt, depth_s);
proxy->smp_unlock();
}
return depth_v[proxy->smp_thread_id()];
}
return depth_s;
}
void colvarmodule::set_error_bits(int code)
{
proxy->smp_lock();
errorCode |= code;
errorCode |= COLVARS_ERROR;
proxy->smp_unlock();
}
void colvarmodule::clear_error()
{
proxy->smp_lock();
errorCode = COLVARS_OK;
proxy->smp_unlock();
}
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(),
@ -1078,6 +1244,7 @@ void cvm::fatal_error(std::string const &message)
proxy->fatal_error(message);
}
void cvm::exit(std::string const &message)
{
proxy->exit(message);
@ -1232,11 +1399,13 @@ 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;
size_t colvarmodule::depth_s = 0;
std::vector<size_t> colvarmodule::depth_v(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;
bool colvarmodule::scripting_after_biases = true;
// file name prefixes
std::string colvarmodule::output_prefix = "";

81
lib/colvars/colvarmodule.h
View File

@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2016-03-08"
#define COLVARS_VERSION "2016-03-30"
#endif
#ifndef COLVARS_DEBUG
@ -74,6 +74,9 @@ private:
public:
/// Base class to handle mutual dependencies of most objects
class deps;
friend class colvarproxy;
// TODO colvarscript should be unaware of colvarmodule's internals
friend class colvarscript;
@ -106,20 +109,21 @@ public:
/// Module-wide error state
/// see constants at the top of this file
protected:
static int errorCode;
static inline void set_error_bits(int code)
{
errorCode |= code;
errorCode |= COLVARS_ERROR;
}
public:
static void set_error_bits(int code);
static inline int get_error()
{
return errorCode;
}
static inline void clear_error()
{
errorCode = 0;
}
static void clear_error();
/// Current step number
static long it;
@ -167,6 +171,12 @@ public:
}
*/
/// Collective variables to be calculated on different threads;
/// colvars with multple items (e.g. multiple active CVCs) are duplicated
std::vector<colvar *> colvars_smp;
/// Indexes of the items to calculate for each colvar
std::vector<int> colvars_smp_items;
/// Array of collective variable biases
static std::vector<colvarbias *> biases;
/// \brief Number of ABF biases initialized (in normal conditions
@ -269,6 +279,10 @@ public:
/// Write explanatory labels in the trajectory file
std::ostream & write_traj_label(std::ostream &os);
/// Write all trajectory files
int write_traj_files();
/// Write all restart files
int write_restart_files();
/// Write all FINAL output files
int write_output_files();
/// Backup a file before writing it
@ -300,11 +314,21 @@ public:
//// Share among replicas.
int bias_share(std::string const &bias_name);
/// Calculate collective variables and biases
/// Main worker function
int calc();
/// Calculate collective variables
int calc_colvars();
/// Calculate biases
int calc_biases();
/// Integrate bias and restraint forces, send colvar forces to atoms
int update_colvar_forces();
/// Perform analysis
int analyze();
/// \brief Read a collective variable trajectory (post-processing
/// only, not called at runtime)
int read_traj(char const *traj_filename,
@ -480,18 +504,18 @@ protected:
/// 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;
protected:
/// Use scripted colvars forces?
static bool use_scripted_forces;
/// Counter for the current depth in the object hierarchy (useg e.g. in output)
static size_t depth_s;
/// Thread-specific depth
static std::vector<size_t> depth_v;
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;
/// Get the current object depth in the hierarchy
static size_t & depth();
/// Increase the depth (number of indentations in the output)
static void increase_depth();
@ -499,7 +523,24 @@ public:
/// Decrease the depth (number of indentations in the output)
static void decrease_depth();
static inline bool scripted_forces() { return use_scripted_forces; }
static inline bool scripted_forces()
{
return use_scripted_forces;
}
/// Use scripted colvars forces?
static bool use_scripted_forces;
/// Wait for all biases before calculating scripted forces?
static bool scripting_after_biases;
/// Calculate the energy and forces of scripted biases
int calc_scripted_forces();
/// \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;
};

75
lib/colvars/colvarproxy.h
View File

@ -113,6 +113,62 @@ protected:
public:
// ***************** SHARED-MEMORY PARALLELIZATION *****************
/// Whether or not threaded parallelization is available
virtual int smp_enabled()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Distribute calculation of colvars (and their components) across threads
virtual int smp_colvars_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Distribute calculation of biases across threads
virtual int smp_biases_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread
virtual int smp_biases_script_loop()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Index of this thread
virtual int smp_thread_id()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Number of threads sharing this address space
virtual int smp_num_threads()
{
return COLVARS_NOT_IMPLEMENTED;
}
/// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing
virtual int smp_lock()
{
return COLVARS_OK;
}
/// Attempt to lock the proxy's shared data
virtual int smp_trylock()
{
return COLVARS_OK;
}
/// Release the lock
virtual int smp_unlock()
{
return COLVARS_OK;
}
// **************** MULTIPLE REPLICAS COMMUNICATION ****************
// Replica exchange commands:
@ -466,9 +522,9 @@ protected:
public:
/// \brief Whether this proxy implementation has capability for scalable groups
virtual bool has_scalable_groups() const
virtual int scalable_group_coms()
{
return false;
return COLVARS_NOT_IMPLEMENTED;
}
/// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
@ -476,6 +532,7 @@ public:
inline int add_atom_group_slot(int atom_group_id)
{
atom_groups_ids.push_back(atom_group_id);
atom_groups_ncopies.push_back(1);
atom_groups_masses.push_back(1.0);
atom_groups_charges.push_back(0.0);
atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0));
@ -496,18 +553,18 @@ public:
/// \brief Used by the atom_group class destructor
virtual void clear_atom_group(int index)
{
if (cvm::debug()) {
log("Trying to remove/disable atom group number "+cvm::to_str(index)+"\n");
}
if (((size_t) index) >= atom_groups_ids.size()) {
cvm::error("Error: trying to disable an atom group that was not previously requested.\n",
INPUT_ERROR);
}
atom_groups_ids.erase(atom_groups_ids.begin()+index);
atom_groups_masses.erase(atom_groups_masses.begin()+index);
atom_groups_charges.erase(atom_groups_charges.begin()+index);
atom_groups_coms.erase(atom_groups_coms.begin()+index);
atom_groups_total_forces.erase(atom_groups_total_forces.begin()+index);
atom_groups_applied_forces.erase(atom_groups_applied_forces.begin()+index);
atom_groups_new_colvar_forces.erase(atom_groups_new_colvar_forces.begin()+index);
if (atom_groups_ncopies[index] > 0) {
atom_groups_ncopies[index] -= 1;
}
}
/// Get the numeric ID of the given atom group (for the MD program)

12
lib/colvars/colvarscript.cpp
View File

@ -221,7 +221,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
if (subcmd == "update") {
cv->calc();
cv->update();
cv->update_forces_energy();
result = (cv->value()).to_simple_string();
return COLVARS_OK;
}
@ -243,6 +243,16 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
return COLVARS_OK;
}
if (subcmd == "getappliedforce") {
result = (cv->bias_force()).to_simple_string();
return COLVARS_OK;
}
if (subcmd == "getsystemforce") {
result = (cv->system_force()).to_simple_string();
return COLVARS_OK;
}
if (subcmd == "addforce") {
if (argc < 4) {
result = "addforce: missing parameter: force value\n" + help_string();

85
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -9,6 +9,8 @@
#include "fix_colvars.h"
#include "colvarmodule.h"
#include "colvar.h"
#include "colvarbias.h"
#include "colvaratoms.h"
#include "colvarproxy.h"
#include "colvarproxy_lammps.h"
@ -109,11 +111,20 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
if (restart_output_prefix_str.rfind(".*") != std::string::npos)
restart_output_prefix_str.erase(restart_output_prefix_str.rfind(".*"),2);
#if defined(_OPENMP)
if (smp_thread_id() == 0) {
omp_init_lock(&smp_lock_state);
}
#endif
// initialize multi-replica support, if available
if (replica_enabled()) {
MPI_Comm_rank(inter_comm, &inter_me);
MPI_Comm_size(inter_comm, &inter_num);
}
if (cvm::debug())
log("Done initializing the colvars proxy object.\n");
}
@ -320,6 +331,80 @@ int colvarproxy_lammps::backup_file(char const *filename)
}
}
#if defined(_OPENMP)
// SMP support
int colvarproxy_lammps::smp_enabled()
{
return COLVARS_OK;
}
int colvarproxy_lammps::smp_colvars_loop()
{
colvarmodule *cv = this->colvars;
#pragma omp parallel for
for (size_t i = 0; i < cv->colvars_smp.size(); i++) {
if (cvm::debug()) {
cvm::log("Calculating colvar \""+cv->colvars_smp[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n");
}
cv->colvars_smp[i]->calc_cvcs(cv->colvars_smp_items[i], 1);
}
return cvm::get_error();
}
int colvarproxy_lammps::smp_biases_loop()
{
colvarmodule *cv = this->colvars;
#pragma omp parallel for
for (size_t i = 0; i < cv->biases.size(); i++) {
if (cvm::debug()) {
cvm::log("Calculating bias \""+cv->biases[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n");
}
cv->biases[i]->update();
}
return cvm::get_error();
}
int colvarproxy_lammps::smp_thread_id()
{
return omp_get_thread_num();
}
int colvarproxy_lammps::smp_num_threads()
{
return omp_get_num_threads();
}
int colvarproxy_lammps::smp_lock()
{
omp_set_lock(&smp_lock_state);
return COLVARS_OK;
}
int colvarproxy_lammps::smp_trylock()
{
return omp_test_lock(&smp_lock_state) ? COLVARS_OK : COLVARS_ERROR;
}
int colvarproxy_lammps::smp_unlock()
{
omp_unset_lock(&smp_lock_state);
return COLVARS_OK;
}
#endif
// multi-replica support
void colvarproxy_lammps::replica_comm_barrier() {

23
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -5,6 +5,7 @@
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarvalue.h"
#include "lammps.h"
#include "domain.h"
@ -16,8 +17,12 @@
#include <vector>
#include <iostream>
#if defined(_OPENMP)
#include <omp.h>
#endif
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2016-02-28"
#define COLVARPROXY_VERSION "2016-03-24"
#endif
/* struct for packed data communication of coordinates and forces. */
@ -130,6 +135,22 @@ class colvarproxy_lammps : public colvarproxy {
// implementation of optional methods from base class
public:
#if defined(_OPENMP)
// SMP support
int smp_enabled();
int smp_colvars_loop();
int smp_biases_loop();
int smp_thread_id();
int smp_num_threads();
protected:
omp_lock_t smp_lock_state;
public:
int smp_lock();
int smp_trylock();
int smp_unlock();
#endif
// Multi-replica support
// Indicate if multi-replica support is available and active
virtual bool replica_enabled() { return (inter_comm != MPI_COMM_NULL); }