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

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

Browse Source
This commit is contained in:
sjplimp
2016-04-15 16:07:01 +00:00
parent 32509da721
commit 212a955285
31 changed files with 1713 additions and 1493 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/colvars/Makefile.fermi
View File

@ -19,7 +19,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.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/

845
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;
}

139
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;
@ -280,6 +280,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 +293,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,11 +372,7 @@ 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 (b_scalable) {
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
// checks of doubly-counted atoms have been handled by add_atom() already
@ -408,6 +403,10 @@ int cvm::atom_group::parse(std::string const &conf)
}
}
if (is_enabled(f_ag_scalable) && !b_dummy) {
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
parse_error |= parse_fitting_options(group_conf);
// TODO move this to colvarparse object
@ -449,7 +448,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 +489,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 +519,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 +557,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));
}
@ -750,15 +749,24 @@ void cvm::atom_group::read_positions()
int cvm::atom_group::calc_required_properties()
{
if (b_dummy) return COLVARS_OK;
// TODO check if the com is needed?
calc_center_of_mass();
if (!b_scalable) {
// TODO check if calc_center_of_geometry() is needed without a fit?
calc_center_of_geometry();
if (!is_enabled(f_ag_scalable)) {
if (b_center || b_rotate) {
if (ref_pos_group) {
ref_pos_group->calc_center_of_geometry();
}
calc_apply_roto_translation();
// update COM and COG after fitting
calc_center_of_geometry();
calc_center_of_mass();
if (ref_pos_group) {
ref_pos_group->calc_center_of_geometry();
}
}
}
@ -767,37 +775,49 @@ int cvm::atom_group::calc_required_properties()
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void cvm::atom_group::calc_apply_roto_translation()
{
atom_group *fit_group = ref_pos_group ? ref_pos_group : this;
if (b_center) {
// center on the origin first
cvm::atom_pos const cog = fit_group->center_of_geometry();
cvm::atom_pos const cog = ref_pos_group ?
ref_pos_group->center_of_geometry() : this->center_of_geometry();
apply_translation(-1.0 * cog);
if (ref_pos_group) {
ref_pos_group->apply_translation(-1.0 * cog);
}
}
if (b_rotate) {
// rotate the group (around the center of geometry if b_center is
// true, around the origin otherwise)
rot.calc_optimal_rotation(fit_group->positions(), ref_pos);
rot.calc_optimal_rotation(ref_pos_group ?
ref_pos_group->positions() :
this->positions(),
ref_pos);
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
cvm::atom_iter ai;
for (ai = this->begin(); ai != this->end(); ai++) {
ai->pos = rot.rotate(ai->pos);
}
if (ref_pos_group) {
for (ai = ref_pos_group->begin(); ai != ref_pos_group->end(); ai++) {
ai->pos = rot.rotate(ai->pos);
}
}
}
if (b_center) {
// align with the center of geometry of ref_pos
apply_translation(ref_pos_cog);
// update the center of geometry for external use
cog = ref_pos_cog;
if (ref_pos_group) {
ref_pos_group->apply_translation(ref_pos_cog);
}
// recalculate the center of mass
calc_center_of_mass();
}
// update of COM and COG is done from the calling routine
}
void cvm::atom_group::apply_translation(cvm::rvector const &t)
{
if (b_dummy) {
@ -805,7 +825,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;
}
@ -815,23 +835,6 @@ void cvm::atom_group::apply_translation(cvm::rvector const &t)
}
}
void cvm::atom_group::apply_rotation(cvm::rotation const &rot)
{
if (b_dummy) {
cvm::error("Error: cannot rotate the coordinates of a dummy atom group.\n", INPUT_ERROR);
return;
}
if (b_scalable) {
cvm::error("Error: cannot rotate the coordinates of a scalable atom group.\n", INPUT_ERROR);
return;
}
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->pos = rot.rotate(ai->pos - center_of_geometry()) + center_of_geometry();
}
}
void cvm::atom_group::read_velocities()
{
@ -852,6 +855,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 +896,10 @@ int cvm::atom_group::calc_center_of_mass()
{
if (b_dummy) {
com = dummy_atom_pos;
} else if (b_scalable) {
if (cvm::debug()) {
cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n");
}
} else if (is_enabled(f_ag_scalable)) {
com = (cvm::proxy)->get_atom_group_com(index);
} else {
com.reset();
@ -922,7 +930,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 +1001,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 +1022,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 +1043,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 +1064,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 +1078,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 +1087,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 +1114,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 +1166,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 +1219,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;

16
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();
@ -325,10 +331,6 @@ public:
/// \brief Move all positions
void apply_translation(cvm::rvector const &t);
/// \brief Rotate all positions around the center of geometry
void apply_rotation(cvm::rotation const &q);
/// \brief Get the current velocities; this must be called always
/// *after* read_positions(); if b_rotate is defined, the same
/// rotation applied to the coordinates will be used

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

39
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
colvars[i]->enable(f_cv_grid);
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_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++) {

11
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());
@ -53,6 +56,10 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf,
size_t i;
if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
if (colvar_centers.size() != colvars.size()) {
cvm::error("Error: number of target centers does not match "
"that of collective variables.\n");
}
b_chg_centers = true;
for (i = 0; i < target_centers.size(); i++) {
target_centers[i].apply_constraints();
@ -167,7 +174,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 +340,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);

170
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,
cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
char const *group_key,
cvm::atom_group &group,
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:

239
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");
}
@ -29,17 +28,19 @@ 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))
{
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,19 +252,17 @@ 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);
}
@ -287,10 +272,6 @@ 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))
{
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);
}
}

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");

310
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,82 +468,169 @@ int colvarmodule::calc() {
cvm::to_str(cvm::step_absolute())+"\n");
}
// calculate collective variables and their gradients
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();
}
error_code |= calc_biases();
error_code |= update_colvar_forces();
if (cvm::b_analysis) {
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)->calc();
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;
}
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;
}
}
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++) {
total_bias_energy += (*bi)->update();
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 (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->reset_bias_force();
}
for (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
// Run user force script, if provided,
// potentially adding scripted forces to the colvars
if (use_scripted_forces) {
int res;
res = proxy->run_force_callback();
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Colvar forces scripts are not implemented.");
return COLVARS_ERROR;
}
if (res != COLVARS_OK) {
cvm::error("Error running user colvar forces script");
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
if (cvm::b_analysis) {
// perform runtime analysis of colvars and biases
if (cvm::debug() && biases.size())
cvm::log("Perform runtime analyses.\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
(*cvi)->analyze();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
for (bi = biases.begin(); bi != biases.end(); bi++) {
(*bi)->analyze();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
cvm::decrease_depth();
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,8 +662,30 @@ int colvarmodule::calc() {
}
cvm::decrease_depth();
// write restart file, if needed
if (restart_out_freq && restart_out_name.size()) {
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;
}
return COLVARS_OK;
}
int colvarmodule::write_restart_files()
{
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Writing the state file \""+
@ -580,11 +696,13 @@ int colvarmodule::calc() {
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);
}
int colvarmodule::write_traj_files()
{
if (!cv_traj_os.is_open()) {
open_traj_file(cv_traj_name);
}
@ -607,7 +725,6 @@ int colvarmodule::calc() {
cv_traj_os.flush();
}
}
} // end if (cv_traj_freq)
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@ -643,6 +760,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 +771,23 @@ int colvarmodule::setup()
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
colvarmodule::~colvarmodule()
{
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 +797,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();
@ -983,7 +1107,7 @@ int colvarmodule::open_traj_file(std::string const &file_name)
int colvarmodule::close_traj_file()
{
if (cv_traj_os.good()) {
if (cv_traj_os.is_open()) {
cv_traj_os.close();
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -1048,28 +1172,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 +1245,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 +1400,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-04-14"
#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();