CVS diff colvar.C

Difference for src/colvar.C from version 1.25 to 1.26

version 1.25

version 1.26

Line 1

/// -*- c++ -*-

// -*- c++ -*-

#include "colvarmodule.h"

#include "colvarvalue.h"

Line 19

colvar::colvar(std::string const &conf)

: colvarparse(conf)

{

size_t i;

cvm::log("Initializing a new collective variable.\n");

int error_code = COLVARS_OK;

get_keyval(conf, "name", this->name,

(std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));

Line 32

return;

}

// all tasks disabled by default

// Initialize dependency members

for (i = 0; i < task_ntot; i++) {

// Could be a function defined in a different source file, for space?

tasks[i] = false;

}

this->description = "colvar " + this->name;

// Initialize static array once and for all

init_cv_requires();

kinetic_energy = 0.0;

potential_energy = 0.0;

// read the configuration and set up corresponding instances, for

error_code |= init_components(conf);

// each type of component implemented

if (error_code != COLVARS_OK) return;

#define initialize_components(def_desc,def_config_key,def_class_name) \

{ \

size_t def_count = 0; \

std::string def_conf = ""; \

size_t pos = 0; \

while ( this->key_lookup(conf, \

def_config_key, \

def_conf, \

pos) ) { \

if (!def_conf.size()) continue; \

cvm::log("Initializing " \

"a new \""+std::string(def_config_key)+"\" component"+ \

(cvm::debug() ? ", with configuration:\n"+def_conf \

: ".\n")); \

cvm::increase_depth(); \

cvc *cvcp = new colvar::def_class_name(def_conf); \

if (cvcp != NULL) { \

cvcs.push_back(cvcp); \

cvcp->check_keywords(def_conf, def_config_key); \

if (cvm::get_error()) { \

cvm::error("Error: in setting up component \"" \

def_config_key"\".\n"); \

return; \

} \

cvm::decrease_depth(); \

} else { \

cvm::error("Error: in allocating component \"" \

def_config_key"\".\n", \

MEMORY_ERROR); \

return; \

} \

if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) { \

if ( (cvcp->function_type != std::string("distance_z")) && \

(cvcp->function_type != std::string("dihedral")) && \

(cvcp->function_type != std::string("spin_angle")) ) { \

cvm::error("Error: invalid use of period and/or " \

"wrapAround in a \""+ \

std::string(def_config_key)+ \

"\" component.\n"+ \

"Period: "+cvm::to_str(cvcp->period) + \

" wrapAround: "+cvm::to_str(cvcp->wrap_center), \

INPUT_ERROR); \

return; \

} \

if ( ! cvcs.back()->name.size()){ \

std::ostringstream s; \

s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;\

cvcs.back()->name = s.str(); \

/* pad cvc number for correct ordering when sorting by name */\

} \

if (cvm::debug()) \

cvm::log("Done initializing a \""+ \

std::string(def_config_key)+ \

"\" component"+ \

(cvm::debug() ? \

", named \""+cvcs.back()->name+"\"" \

: "")+".\n"); \

def_conf = ""; \

if (cvm::debug()) \

cvm::log("Parsed "+cvm::to_str(cvcs.size())+ \

" components at this time.\n"); \

} \

}

initialize_components("distance", "distance", distance);

initialize_components("distance vector", "distanceVec", distance_vec);

initialize_components("Cartesian coordinates", "cartesian", cartesian);

initialize_components("distance vector "

"direction", "distanceDir", distance_dir);

initialize_components("distance projection "

"on an axis", "distanceZ", distance_z);

initialize_components("distance projection "

"on a plane", "distanceXY", distance_xy);

initialize_components("average distance weighted by inverse power",

"distanceInv", distance_inv);

initialize_components("N1xN2-long vector of pairwise distances",

"distancePairs", distance_pairs);

initialize_components("coordination "

"number", "coordNum", coordnum);

initialize_components("self-coordination "

"number", "selfCoordNum", selfcoordnum);

initialize_components("angle", "angle", angle);

initialize_components("dihedral", "dihedral", dihedral);

initialize_components("hydrogen bond", "hBond", h_bond);

// initialize_components ("alpha helix", "alphaDihedrals", alpha_dihedrals);

initialize_components("alpha helix", "alpha", alpha_angles);

initialize_components("dihedral principal "

"component", "dihedralPC", dihedPC);

initialize_components("orientation", "orientation", orientation);

initialize_components("orientation "

"angle", "orientationAngle",orientation_angle);

initialize_components("orientation "

"projection", "orientationProj",orientation_proj);

initialize_components("tilt", "tilt", tilt);

initialize_components("spin angle", "spinAngle", spin_angle);

initialize_components("RMSD", "rmsd", rmsd);

// initialize_components ("logarithm of MSD", "logmsd", logmsd);

initialize_components("radius of "

"gyration", "gyration", gyration);

initialize_components("moment of "

"inertia", "inertia", inertia);

initialize_components("moment of inertia around an axis",

"inertiaZ", inertia_z);

initialize_components("eigenvector", "eigenvector", eigenvector);

if (!cvcs.size()) {

size_t i;

cvm::error("Error: no valid components were provided "

"for this collective variable.\n",

INPUT_ERROR);

return;

}

cvm::log("All components initialized.\n");

// Setup colvar as scripted function of components

if (get_keyval(conf, "scriptedFunction", scripted_function,

"", colvarparse::parse_silent)) {

enable(task_scripted);

// Make feature available only on user request

provide(f_cv_scripted);

enable(f_cv_scripted);

cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");

std::string type_str;

Line 189

Line 72

cvm::error("Could not parse scripted colvar type.");

return;

}

x_reported.type(x.type());

cvm::log(std::string("Expecting colvar value of type ")

+ colvarvalue::type_desc(x.type()));

Line 202

Line 85

x.vector1d_value.resize(size);

}

x_reported.type(x);

// Sort array of cvcs based on their names

// Note: default CVC names are in input order for same type of CVC

std::sort(cvcs.begin(), cvcs.end(), compare);

Line 218

Line 103

for (i = 0; i < cvcs.size(); i++) {

sorted_cvc_values.push_back(&(cvcs[i]->value()));

}

b_homogeneous = false;

// Scripted functions are deemed non-periodic

b_periodic = false;

period = 0.0;

b_inverse_gradients = false;

b_Jacobian_force = false;

}

if (!tasks[task_scripted]) {

if (!is_enabled(f_cv_scripted)) {

colvarvalue const &cvc_value = (cvcs[0])->value();

if (cvm::debug())

cvm::log ("This collective variable is a "+

Line 241

Line 119

// If using scripted biases, any colvar may receive bias forces

// and will need its gradient

if (cvm::scripted_forces()) {

enable(task_gradients);

enable(f_cv_gradient);

}

// check for linear combinations

{

b_linear = !tasks[task_scripted];

bool lin = !is_enabled(f_cv_scripted);

for (i = 0; i < cvcs.size(); i++) {

if ((cvcs[i])->b_debug_gradients)

enable(task_gradients);

// FIXME this is a reverse dependency, ie. cv feature depends on cvc flag

// need to clarify this case

// if ((cvcs[i])->b_debug_gradients)

// enable(task_gradients);

if ((cvcs[i])->sup_np != 1) {

if (cvm::debug() && b_linear)

if (cvm::debug() && lin)

cvm::log("Warning: You are using a non-linear polynomial "

"combination to define this collective variable, "

"some biasing methods may be unavailable.\n");

b_linear = false;

lin = false;

if ((cvcs[i])->sup_np < 0) {

cvm::log("Warning: you chose a negative exponent in the combination; "

Line 266

Line 147

}

feature_states[f_cv_linear]->enabled = lin;

}

// Colvar is homogeneous iff:

// Colvar is homogeneous if:

// - it is not scripted

// - it is linear (hence not scripted)

// - it is linear

// - all cvcs have coefficient 1 or -1

// i.e. sum or difference of cvcs

{

b_homogeneous = !tasks[task_scripted] && b_linear;

bool homogeneous = is_enabled(f_cv_linear);

for (i = 0; i < cvcs.size(); i++) {

if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {

b_homogeneous = false;

homogeneous = false;

}

feature_states[f_cv_homogeneous]->enabled = homogeneous;

}

// Colvar is deemed periodic iff:

// - it is homogeneous

// - all cvcs are periodic

// - all cvcs have the same period

b_periodic = cvcs[0]->b_periodic && b_homogeneous;

b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous);

period = cvcs[0]->period;

for (i = 1; i < cvcs.size(); i++) {

if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {

Line 293

Line 176

period = 0.0;

}

feature_states[f_cv_periodic]->enabled = b_periodic;

// these will be set to false if any of the cvcs has them false

// check that cvcs are compatible

b_inverse_gradients = !tasks[task_scripted];

b_Jacobian_force = !tasks[task_scripted];

// check the available features of each cvc

for (i = 0; i < cvcs.size(); i++) {

if ((cvcs[i])->b_periodic && !b_periodic) {

cvm::log("Warning: although this component is periodic, the colvar will "

Line 307

Line 188

"you know what you are doing!");

}

if (! (cvcs[i])->b_inverse_gradients)

b_inverse_gradients = false;

if (! (cvcs[i])->b_Jacobian_derivative)

b_Jacobian_force = false;

// components may have different types only for scripted functions

if (!tasks[task_scripted] && (colvarvalue::check_types(cvcs[i]->value(),

if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),

cvcs[0]->value())) ) {

cvm::error("ERROR: you are definining this collective variable "

"by using components of different types. "

Line 324

Line 199

}

active_cvc_square_norm = 0.;

for (i = 0; i < cvcs.size(); i++) {

active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;

}

// at this point, the colvar's type is defined

f.type(value());

f_accumulated.type(value());

fb.type(value());

get_keyval(conf, "width", width, 1.0);

Line 333

Line 214

cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);

}

// NOTE: not porting wall stuff to new deps, as this will change to a separate bias

// the grid functions will wait a little as well

lower_boundary.type(value());

lower_wall.type(value());

upper_boundary.type(value());

upper_wall.type(value());

if (value().type() == colvarvalue::type_scalar) {

feature_states[f_cv_scalar]->enabled = (value().type() == colvarvalue::type_scalar);

if (is_enabled(f_cv_scalar)) {

if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {

enable(task_lower_boundary);

provide(f_cv_lower_boundary);

enable(f_cv_lower_boundary);

}

get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0);

if (lower_wall_k > 0.0) {

get_keyval(conf, "lowerWall", lower_wall, lower_boundary);

enable(task_lower_wall);

enable(f_cv_lower_wall);

}

if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {

enable(task_upper_boundary);

provide(f_cv_upper_boundary);

enable(f_cv_upper_boundary);

}

get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0);

if (upper_wall_k > 0.0) {

get_keyval(conf, "upperWall", upper_wall, upper_boundary);

enable(task_upper_wall);

enable(f_cv_upper_wall);

}

if (tasks[task_lower_boundary]) {

if (is_enabled(f_cv_lower_boundary)) {

get_keyval(conf, "hardLowerBoundary", hard_lower_boundary, false);

}

if (tasks[task_upper_boundary]) {

if (is_enabled(f_cv_upper_boundary)) {

get_keyval(conf, "hardUpperBoundary", hard_upper_boundary, false);

}

// consistency checks for boundaries and walls

if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) {

if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) {

if (lower_boundary >= upper_boundary) {

cvm::error("Error: the upper boundary, "+

cvm::to_str(upper_boundary)+

Line 380

Line 268

}

if (tasks[task_lower_wall] && tasks[task_upper_wall]) {

if (is_enabled(f_cv_lower_wall) && is_enabled(f_cv_upper_wall)) {

if (lower_wall >= upper_wall) {

cvm::error("Error: the upper wall, "+

cvm::to_str(upper_wall)+

Line 388

Line 276

cvm::to_str(lower_wall)+".\n",

INPUT_ERROR);

}

if (dist2(lower_wall, upper_wall) < 1.0E-12) {

cvm::log("Lower wall and upper wall are equal "

"in the periodic domain of the colvar: disabling walls.\n");

disable(task_lower_wall);

disable(task_upper_wall);

}

get_keyval(conf, "expandBoundaries", expand_boundaries, false);

Line 410

Line 291

INPUT_ERROR);

}

get_keyval(conf, "timeStepFactor", time_step_factor, 1);

{

bool b_extended_lagrangian;

bool b_extended_Lagrangian;

get_keyval(conf, "extendedLagrangian", b_extended_lagrangian, false);

get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);

if (b_extended_lagrangian) {

if (b_extended_Lagrangian) {

cvm::real temp, tolerance, period;

cvm::log("Enabling the extended Lagrangian term for colvar \""+

this->name+"\".\n");

enable(task_extended_lagrangian);

// Make feature available only on user request

provide(f_cv_extended_Lagrangian);

enable(f_cv_extended_Lagrangian);

provide(f_cv_Langevin);

xr.type(value());

vr.type(value());

Line 455

Line 341

bool b_output_energy;

get_keyval(conf, "outputEnergy", b_output_energy, false);

if (b_output_energy) {

enable(task_output_energy);

enable(f_cv_output_energy);

}

Line 464

Line 350

cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);

}

if (ext_gamma != 0.0) {

enable(task_langevin);

enable(f_cv_Langevin);

ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1

ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());

}

Line 475

Line 361

bool b_output_value;

get_keyval(conf, "outputValue", b_output_value, true);

if (b_output_value) {

enable(task_output_value);

enable(f_cv_output_value);

}

Line 483

Line 369

bool b_output_velocity;

get_keyval(conf, "outputVelocity", b_output_velocity, false);

if (b_output_velocity) {

enable(task_output_velocity);

enable(f_cv_output_velocity);

}

{

bool b_output_system_force;

bool temp;

get_keyval(conf, "outputSystemForce", b_output_system_force, false);

if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {

if (b_output_system_force) {

cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"

enable(task_output_system_force);

"The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);

return;

}

{

get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);

bool b_output_applied_force;

get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);

get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);

get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);

if (b_output_applied_force) {

enable(task_output_applied_force);

// Start in active state by default

}

enable(f_cv_active);

}

// Make sure dependency side-effects are correct

refresh_deps();

x_old.type(value());

v_fdiff.type(value());

v_reported.type(value());

fj.type(value());

ft.type(value());

ft_reported.type(value());

f_old.type(value());

f_old.reset();

if (cvm::b_analysis)

parse_analysis(conf);

Line 512

Line 409

// read the configuration and set up corresponding instances, for

// each type of component implemented

template<typename def_class_name> int colvar::init_components_type(std::string const &conf,

char const *def_desc,

char const *def_config_key)

{

size_t def_count = 0;

std::string def_conf = "";

size_t pos = 0;

while ( this->key_lookup(conf,

def_config_key,

def_conf,

pos) ) {

if (!def_conf.size()) continue;

cvm::log("Initializing "

"a new \""+std::string(def_config_key)+"\" component"+

(cvm::debug() ? ", with configuration:\n"+def_conf

: ".\n"));

cvm::increase_depth();

cvc *cvcp = new def_class_name(def_conf);

if (cvcp != NULL) {

cvcs.push_back(cvcp);

cvcp->check_keywords(def_conf, def_config_key);

if (cvm::get_error()) {

cvm::error("Error: in setting up component \""+

std::string(def_config_key)+"\".\n", INPUT_ERROR);

return INPUT_ERROR;

}

cvm::decrease_depth();

} else {

cvm::error("Error: in allocating component \""+

std::string(def_config_key)+"\".\n",

MEMORY_ERROR);

return MEMORY_ERROR;

}

if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) {

if ( (cvcp->function_type != std::string("distance_z")) &&

(cvcp->function_type != std::string("dihedral")) &&

(cvcp->function_type != std::string("spin_angle")) ) {

cvm::error("Error: invalid use of period and/or "

"wrapAround in a \""+

std::string(def_config_key)+

"\" component.\n"+

"Period: "+cvm::to_str(cvcp->period) +

" wrapAround: "+cvm::to_str(cvcp->wrap_center),

INPUT_ERROR);

return INPUT_ERROR;

}

if ( ! cvcs.back()->name.size()) {

std::ostringstream s;

s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;

cvcs.back()->name = s.str();

/* pad cvc number for correct ordering when sorting by name */

}

cvcs.back()->setup();

if (cvm::debug()) {

cvm::log("Done initializing a \""+

std::string(def_config_key)+

"\" component"+

(cvm::debug() ?

", named \""+cvcs.back()->name+"\""

: "")+".\n");

}

def_conf = "";

if (cvm::debug()) {

cvm::log("Parsed "+cvm::to_str(cvcs.size())+

" components at this time.\n");

}

return COLVARS_OK;

}

int colvar::init_components(std::string const &conf)

{

int error_code = COLVARS_OK;

error_code |= init_components_type<distance>(conf, "distance", "distance");

error_code |= init_components_type<distance_vec>(conf, "distance vector", "distanceVec");

error_code |= init_components_type<cartesian>(conf, "Cartesian coordinates", "cartesian");

error_code |= init_components_type<distance_dir>(conf, "distance vector "

"direction", "distanceDir");

error_code |= init_components_type<distance_z>(conf, "distance projection "

"on an axis", "distanceZ");

error_code |= init_components_type<distance_xy>(conf, "distance projection "

"on a plane", "distanceXY");

error_code |= init_components_type<distance_inv>(conf, "average distance "

"weighted by inverse power", "distanceInv");

error_code |= init_components_type<distance_pairs>(conf, "N1xN2-long vector "

"of pairwise distances", "distancePairs");

error_code |= init_components_type<coordnum>(conf, "coordination "

"number", "coordNum");

error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "

"number", "selfCoordNum");

error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "

"number", "groupCoord");

error_code |= init_components_type<angle>(conf, "angle", "angle");

error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");

error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");

error_code |= init_components_type<h_bond>(conf, "hydrogen bond", "hBond");

error_code |= init_components_type<alpha_angles>(conf, "alpha helix", "alpha");

error_code |= init_components_type<dihedPC>(conf, "dihedral "

"principal component", "dihedralPC");

error_code |= init_components_type<orientation>(conf, "orientation", "orientation");

error_code |= init_components_type<orientation_angle>(conf, "orientation "

"angle", "orientationAngle");

error_code |= init_components_type<orientation_proj>(conf, "orientation "

"projection", "orientationProj");

error_code |= init_components_type<tilt>(conf, "tilt", "tilt");

error_code |= init_components_type<spin_angle>(conf, "spin angle", "spinAngle");

error_code |= init_components_type<rmsd>(conf, "RMSD", "rmsd");

error_code |= init_components_type<gyration>(conf, "radius of "

"gyration", "gyration");

error_code |= init_components_type<inertia>(conf, "moment of "

"inertia", "inertia");

error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");

error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");

if (!cvcs.size() || (error_code != COLVARS_OK)) {

cvm::error("Error: no valid components were provided "

"for this collective variable.\n",

INPUT_ERROR);

return INPUT_ERROR;

}

n_active_cvcs = cvcs.size();

cvm::log("All components initialized.\n");

// Store list of children cvcs for dependency checking purposes

for (size_t i = 0; i < cvcs.size(); i++) {

add_child(cvcs[i]);

}

return COLVARS_OK;

}

int colvar::refresh_deps()

{

// If enabled features are changed upstream, the features below should be refreshed

if (is_enabled(f_cv_total_force_calc)) {

cvm::request_total_force();

}

if (is_enabled(f_cv_collect_gradient) && atom_ids.size() == 0) {

build_atom_list();

}

return COLVARS_OK;

}

void colvar::build_atom_list(void)

{

// If atomic gradients are requested, build full list of atom ids from all cvcs

Line 519

Line 573

for (size_t i = 0; i < cvcs.size(); i++) {

for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {

for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {

cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);

temp_id_list.push_back(cvcs[i]->atom_groups[j]->at(k).id);

for (size_t k = 0; k < ag.size(); k++) {

temp_id_list.push_back(ag[k].id);

}

Line 559

Line 614

bool b_runave = false;

if (get_keyval(conf, "runAve", b_runave) && b_runave) {

enable(task_runave);

enable(f_cv_runave);

get_keyval(conf, "runAveLength", runave_length, 1000);

get_keyval(conf, "runAveStride", runave_stride, 1);

Line 574

Line 629

this->name+".runave.traj"));

size_t const this_cv_width = x.output_width(cvm::cv_width);

cvm::backup_file(runave_outfile.c_str());

runave_os.open(runave_outfile.c_str());

runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)

<< " "

Line 587

Line 643

bool b_acf = false;

if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {

enable(task_corrfunc);

enable(f_cv_corrfunc);

std::string acf_colvar_name;

get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);

Line 604

Line 660

acf_type = acf_coor;

} else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {

acf_type = acf_vel;

enable(task_fdiff_velocity);

enable(f_cv_fdiff_velocity);

if (acf_colvar_name.size())

(cvm::colvar_by_name(acf_colvar_name))->enable(task_fdiff_velocity);

(cvm::colvar_by_name(acf_colvar_name))->enable(f_cv_fdiff_velocity);

} else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {

acf_type = acf_p2coor;

} else {

Line 632

Line 688

}

int colvar::enable(colvar::task const &t)

{

switch (t) {

case task_output_system_force:

enable(task_system_force);

break;

case task_report_Jacobian_force:

enable(task_Jacobian_force);

enable(task_system_force);

if (cvm::debug())

cvm::log("Adding the Jacobian force to the system force, "

"rather than applying its opposite silently.\n");

break;

case task_Jacobian_force:

// checks below do not apply to extended-system colvars

if ( !tasks[task_extended_lagrangian] ) {

enable(task_gradients);

if (!b_Jacobian_force) {

cvm::error("Error: colvar \""+this->name+

"\" does not have Jacobian forces implemented.\n",

INPUT_ERROR);

}

if (!b_linear) {

cvm::error("Error: colvar \""+this->name+

"\" must be defined as a linear combination "

"to calculate the Jacobian force.\n",

INPUT_ERROR);

}

if (cvm::debug())

cvm::log("Enabling calculation of the Jacobian force "

"on this colvar.\n");

}

fj.type(value());

break;

case task_system_force:

if (!tasks[task_extended_lagrangian]) {

if (!b_inverse_gradients) {

cvm::error("Error: one or more of the components of "

"colvar \""+this->name+

"\" does not support system force calculation.\n",

INPUT_ERROR);

}

cvm::request_system_force();

}

ft.type(value());

ft_reported.type(value());

break;

case task_output_applied_force:

case task_lower_wall:

case task_upper_wall:

// all of the above require gradients

enable(task_gradients);

break;

case task_fdiff_velocity:

x_old.type(value());

v_fdiff.type(value());

v_reported.type(value());

break;

case task_output_velocity:

enable(task_fdiff_velocity);

break;

case task_grid:

if (value().type() != colvarvalue::type_scalar) {

cvm::error("Cannot calculate a grid for collective variable, \""+

this->name+"\", because its value is not a scalar number.\n",

INPUT_ERROR);

}

break;

case task_extended_lagrangian:

enable(task_gradients);

v_reported.type(value());

break;

case task_lower_boundary:

case task_upper_boundary:

if (value().type() != colvarvalue::type_scalar) {

cvm::error("Error: this colvar is not a scalar value "

"and cannot produce a grid.\n",

INPUT_ERROR);

}

break;

case task_output_value:

case task_runave:

case task_corrfunc:

case task_ntot:

case task_langevin:

case task_output_energy:

case task_scripted:

case task_gradients:

break;

case task_collect_gradients:

if (value().type() != colvarvalue::type_scalar) {

cvm::error("Collecting atomic gradients for non-scalar collective variable \""+

this->name+"\" is not yet implemented.\n",

INPUT_ERROR);

}

enable(task_gradients);

if (atom_ids.size() == 0) {

build_atom_list();

}

break;

}

tasks[t] = true;

return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);

}

void colvar::disable(colvar::task const &t)

{

// check dependencies

switch (t) {

case task_gradients:

disable(task_upper_wall);

disable(task_lower_wall);

disable(task_output_applied_force);

disable(task_system_force);

disable(task_Jacobian_force);

break;

case task_system_force:

disable(task_output_system_force);

break;

case task_Jacobian_force:

disable(task_report_Jacobian_force);

break;

case task_fdiff_velocity:

disable(task_output_velocity);

break;

case task_lower_boundary:

case task_upper_boundary:

disable(task_grid);

break;

case task_extended_lagrangian:

case task_report_Jacobian_force:

case task_output_value:

case task_output_velocity:

case task_output_applied_force:

case task_output_system_force:

case task_runave:

case task_corrfunc:

case task_grid:

case task_lower_wall:

case task_upper_wall:

case task_ntot:

case task_langevin:

case task_output_energy:

case task_collect_gradients:

case task_scripted:

break;

}

tasks[t] = false;

}

void colvar::setup() {

// loop over all components to reset masses of all groups

for (size_t i = 0; i < cvcs.size(); i++) {

for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);

atoms.read_positions();

atoms.setup();

atoms.reset_mass(name,i,ig);

atoms.read_positions();

}

Line 819

Line 703

colvar::~colvar()

{

for (size_t i = 0; i < cvcs.size(); i++) {

// Clear references to this colvar's cvcs as children

delete cvcs[i];

// for dependency purposes

remove_all_children();

for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();

ci != cvcs.rend();

++ci) {

// clear all children of this cvc (i.e. its atom groups)

// because the cvc base class destructor can't do it early enough

// and we don't want to have each cvc derived class do it separately

(*ci)->remove_all_children();

delete *ci;

}

// remove reference to this colvar from the CVM

Line 839

Line 733

// ******************** CALC FUNCTIONS ********************

void colvar::calc()

// Default schedule (everything is serialized)

int colvar::calc()

{

size_t i, ig;

// Note: if anything is added here, it should be added also in the SMP block of calc_colvars()

if (cvm::debug())

int error_code = COLVARS_OK;

cvm::log("Calculating colvar \""+this->name+"\".\n");

if (is_enabled(f_cv_active)) {

error_code |= update_cvc_flags();

if (error_code != COLVARS_OK) return error_code;

error_code |= calc_cvcs();

if (error_code != COLVARS_OK) return error_code;

error_code |= collect_cvc_data();

}

return error_code;

}

// prepare atom groups for calculation

int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)

{

int error_code = COLVARS_OK;

if (cvm::debug())

cvm::log("Collecting data from atom groups.\n");

cvm::log("Calculating colvar \""+this->name+"\", components "+

cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");

// Update the enabled/disabled status of cvcs if necessary

error_code |= check_cvc_range(first_cvc, num_cvcs);

if (cvc_flags.size()) {

if (error_code != COLVARS_OK) {

bool any = false;

return error_code;

for (i = 0; i < cvcs.size(); i++) {

cvcs[i]->b_enabled = cvc_flags[i];

any = any || cvc_flags[i];

}

if (!any) {

cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");

return;

}

cvc_flags.resize(0);

}

for (i = 0; i < cvcs.size(); i++) {

error_code |= calc_cvc_values(first_cvc, num_cvcs);

if (!cvcs[i]->b_enabled) continue;

error_code |= calc_cvc_gradients(first_cvc, num_cvcs);

for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

error_code |= calc_cvc_total_force(first_cvc, num_cvcs);

cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);

error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);

atoms.reset_atoms_data();

atoms.read_positions();

if (cvm::debug())

if (atoms.b_center || atoms.b_rotate) {

cvm::log("Done calculating colvar \""+this->name+"\".\n");

atoms.calc_apply_roto_translation();

}

return error_code;

// each atom group will take care of its own ref_pos_group, if defined

}

//// Don't try to get atom velocities, as no back-end currently implements it

// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {

// for (i = 0; i < cvcs.size(); i++) {

// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

// cvcs[i]->atom_groups[ig]->read_velocities();

// }

if (tasks[task_system_force]) {

int colvar::collect_cvc_data()

for (i = 0; i < cvcs.size(); i++) {

{

for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

if (cvm::debug())

cvcs[i]->atom_groups[ig]->read_system_forces();

cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");

int error_code = COLVARS_OK;

error_code |= collect_cvc_values();

error_code |= collect_cvc_gradients();

error_code |= collect_cvc_total_forces();

error_code |= collect_cvc_Jacobians();

error_code |= calc_colvar_properties();

if (cvm::debug())

cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");

return error_code;

}

int colvar::check_cvc_range(int first_cvc, size_t num_cvcs)

{

if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) {

cvm::error("Error: trying to address a component outside the "

"range defined for colvar \""+name+"\".\n", BUG_ERROR);

return BUG_ERROR;

}

return COLVARS_OK;

}

int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)

{

size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();

size_t i, cvc_count;

// calculate the value of the colvar

if (cvm::debug())

cvm::log("Calculating colvar components.\n");

x.reset();

// First, update component values

// First, calculate component values

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

cvm::increase_depth();

for (i = first_cvc, cvc_count = 0;

(i < cvcs.size()) && (cvc_count < cvc_max_count);

i++) {

if (!cvcs[i]->is_enabled()) continue;

cvc_count++;

(cvcs[i])->read_data();

(cvcs[i])->calc_value();

cvm::decrease_depth();

if (cvm::debug())

cvm::log("Colvar component no. "+cvm::to_str(i+1)+

" within colvar \""+this->name+"\" has value "+

cvm::to_str((cvcs[i])->value(),

cvm::cv_width, cvm::cv_prec)+".\n");

}

cvm::decrease_depth();

return COLVARS_OK;

}

int colvar::collect_cvc_values()

{

x.reset();

size_t i;

// Then combine them appropriately

// combine them appropriately, using either a scripted function or a polynomial

if (tasks[task_scripted]) {

if (is_enabled(f_cv_scripted)) {

// cvcs combined by user script

int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);

if (res == COLVARS_NOT_IMPLEMENTED) {

cvm::error("Scripted colvars are not implemented.");

return;

return COLVARS_NOT_IMPLEMENTED;

}

if (res != COLVARS_OK) {

cvm::error("Error running scripted colvar");

return;

return COLVARS_OK;

}

} else if (x.type() == colvarvalue::type_scalar) {

// polynomial combination allowed

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

x += (cvcs[i])->sup_coeff *

( ((cvcs[i])->sup_np != 1) ?

std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :

(cvcs[i])->value().real_value );

}

} else {

// only linear combination allowed

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

x += (cvcs[i])->sup_coeff * (cvcs[i])->value();

}

Line 945

Line 872

cvm::log("Colvar \""+this->name+"\" has value "+

cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");

if (tasks[task_gradients]) {

return COLVARS_OK;

}

int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)

{

size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();

size_t i, cvc_count;

if (cvm::debug())

cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");

for (i = 0; i < cvcs.size(); i++) {

// calculate the gradients of each component

if (!cvcs[i]->b_enabled) continue;

cvm::increase_depth();

for (i = first_cvc, cvc_count = 0;

(i < cvcs.size()) && (cvc_count < cvc_max_count);

i++) {

if (!cvcs[i]->is_enabled()) continue;

cvc_count++;

if ((cvcs[i])->is_enabled(f_cvc_gradient)) {

(cvcs[i])->calc_gradients();

// if requested, propagate (via chain rule) the gradients above

// to the atoms used to define the roto-translation

for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

// This could be integrated in the CVC base class

for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {

if (cvcs[i]->atom_groups[ig]->b_fit_gradients)

cvcs[i]->atom_groups[ig]->calc_fit_gradients();

if (cvcs[i]->is_enabled(f_cvc_debug_gradient)) {

cvm::log("Debugging gradients for " + cvcs[i]->description);

cvcs[i]->debug_gradients(cvcs[i]->atom_groups[ig]);

}

cvm::decrease_depth();

}

cvm::decrease_depth();

if (cvm::debug())

cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");

}

return COLVARS_OK;

}

if (tasks[task_collect_gradients]) {

if (tasks[task_scripted]) {

int colvar::collect_cvc_gradients()

{

size_t i;

if (is_enabled(f_cv_collect_gradient)) {

if (is_enabled(f_cv_scripted)) {

cvm::error("Collecting atomic gradients is not implemented for "

"scripted colvars.");

"scripted colvars.", COLVARS_NOT_IMPLEMENTED);

return;

return COLVARS_NOT_IMPLEMENTED;

}

// Collect the atomic gradients inside colvar object

Line 980

Line 935

atomic_gradients[a].reset();

}

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

// Coefficient: d(a * x^n) = a * n * x^(n-1) * dx

cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *

std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);

for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {

cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);

// If necessary, apply inverse rotation to get atomic

// gradient in the laboratory frame

if (cvcs[i]->atom_groups[j]->b_rotate) {

if (ag.b_rotate) {

cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse();

cvm::rotation const rot_inv = ag.rot.inverse();

for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {

for (size_t k = 0; k < ag.size(); k++) {

int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),

size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),

cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();

ag[k].id) - atom_ids.begin();

atomic_gradients[a] += coeff *

atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);

rot_inv.rotate(cvcs[i]->atom_groups[j]->at(k).grad);

}

} else {

for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {

for (size_t k = 0; k < ag.size(); k++) {

int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),

size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),

cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();

ag[k].id) - atom_ids.begin();

atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad;

atomic_gradients[a] += coeff * ag[k].grad;

}

return COLVARS_OK;

}

if (tasks[task_system_force] && !tasks[task_extended_lagrangian]) {

// If extended Lagrangian is enabled, system force calculation is trivial

// and done together with integration of the extended coordinate.

if (tasks[task_scripted]) {

int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)

// TODO see if this could reasonably be done in a generic way

{

// from generic inverse gradients

size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();

cvm::error("System force is not implemented for "

size_t i, cvc_count;

"scripted colvars.");

return;

}

if (cvm::debug())

cvm::log("Calculating system force of colvar \""+this->name+"\".\n");

ft.reset();

if (is_enabled(f_cv_total_force_calc)) {

if (cvm::debug())

cvm::log("Calculating total force of colvar \""+this->name+"\".\n");

// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {

// Disabled check to allow for explicit system force calculation

// Disabled check to allow for explicit total force calculation

// even with extended Lagrangian

if (cvm::step_relative() > 0) {

// get from the cvcs the system forces from the PREVIOUS step

for (i = 0; i < cvcs.size(); i++) {

(cvcs[i])->calc_force_invgrads();

// linear combination is assumed

cvm::increase_depth();

ft += (cvcs[i])->system_force() / ((cvcs[i])->sup_coeff * cvm::real(cvcs.size()));

// get from the cvcs the total forces from the PREVIOUS step

for (i = first_cvc, cvc_count = 0;

(i < cvcs.size()) && (cvc_count < cvc_max_count);

i++) {

if (!cvcs[i]->is_enabled()) continue;

cvc_count++;

(cvcs[i])->calc_force_invgrads();

}

cvm::decrease_depth();

}

if (cvm::debug())

cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");

}

return COLVARS_OK;

}

int colvar::collect_cvc_total_forces()

{

if (is_enabled(f_cv_total_force_calc)) {

ft.reset();

if (cvm::step_relative() > 0) {

// get from the cvcs the total forces from the PREVIOUS step

for (size_t i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->is_enabled()) continue;

// linear combination is assumed

ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;

}

if (tasks[task_report_Jacobian_force]) {

if (!is_enabled(f_cv_hide_Jacobian)) {

// add the Jacobian force to the system force, and don't apply any silent

// add the Jacobian force to the total force, and don't apply any silent

// correction internally: biases such as colvarbias_abf will handle it

ft += fj;

}

if (cvm::debug())

return COLVARS_OK;

cvm::log("Done calculating system force of colvar \""+this->name+"\".\n");

}

if (tasks[task_fdiff_velocity]) {

int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs)

{

size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();

if (is_enabled(f_cv_Jacobian)) {

cvm::increase_depth();

size_t i, cvc_count;

for (i = first_cvc, cvc_count = 0;

(i < cvcs.size()) && (cvc_count < cvc_max_count);

i++) {

if (!cvcs[i]->is_enabled()) continue;

cvc_count++;

(cvcs[i])->calc_Jacobian_derivative();

}

cvm::decrease_depth();

}

return COLVARS_OK;

}

int colvar::collect_cvc_Jacobians()

{

if (is_enabled(f_cv_Jacobian)) {

fj.reset();

for (size_t i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->is_enabled()) continue;

// linear combination is assumed

fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;

}

fj *= cvm::boltzmann() * cvm::temperature();

}

return COLVARS_OK;

}

int colvar::calc_colvar_properties()

{

if (is_enabled(f_cv_fdiff_velocity)) {

// calculate the velocity by finite differences

if (cvm::step_relative() == 0)

x_old = x;

Line 1063

Line 1078

}

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

// initialize the restraint center in the first step to the value

// just calculated from the cvcs

Line 1076

Line 1091

// report the restraint center as "value"

x_reported = xr;

v_reported = vr;

// the "system force" with the extended Lagrangian is just the

// the "total force" with the extended Lagrangian is

// harmonic term acting on the extended coordinate

// calculated in update_forces_energy() below

// Note: this is the force for current timestep

ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);

} else {

if (is_enabled(f_cv_subtract_applied_force)) {

// correct the total force only if it has been measured

// TODO add a specific test instead of relying on sq norm

if (ft.norm2() > 0.0) {

ft -= f_old;

}

x_reported = x;

ft_reported = ft;

}

if (cvm::debug())

return COLVARS_OK;

cvm::log("Done calculating colvar \""+this->name+"\".\n");

}

cvm::real colvar::update()

cvm::real colvar::update_forces_energy()

{

if (cvm::debug())

cvm::log("Updating colvar \""+this->name+"\".\n");

Line 1105

Line 1125

// been summed over each bias using this colvar

f += fb;

if (is_enabled(f_cv_Jacobian)) {

if (tasks[task_Jacobian_force]) {

// the instantaneous Jacobian force was not included in the reported total force;

size_t i;

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

cvm::increase_depth();

(cvcs[i])->calc_Jacobian_derivative();

cvm::decrease_depth();

}

size_t ncvc = 0;

fj.reset();

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

// linear combination is assumed

fj += 1.0 / cvm::real((cvcs[i])->sup_coeff) *

(cvcs[i])->Jacobian_derivative();

ncvc++;

}

fj *= (1.0/cvm::real(ncvc)) * cvm::boltzmann() * cvm::temperature();

// the instantaneous Jacobian force was not included in the reported system force;

// instead, it is subtracted from the applied force (silent Jacobian correction)

if (! tasks[task_report_Jacobian_force])

if (is_enabled(f_cv_hide_Jacobian))

f -= fj;

}

if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {

if (tasks[task_extended_lagrangian]) {

cvm::real dt = cvm::dt();

cvm::real f_ext;

// the total force is applied to the fictitious mass, while the

// atoms only feel the harmonic force

// fr: bias force on extended coordinate (without harmonic spring), for output in trajectory

// f_ext: total force on extended coordinate (including harmonic spring)

// f: - initially, external biasing force

// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force

// (note: wall potential is added to f after this block)

fr = f;

f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);

f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);

// leapfrog: starting from x_i, f_i, v_(i-1/2)

vr += (0.5 * dt) * f_ext / ext_mass;

// Because of leapfrog, kinetic energy at time i is approximate

kinetic_energy = 0.5 * ext_mass * vr * vr;

potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);

// leap to v_(i+1/2)

if (tasks[task_langevin]) {

vr -= dt * ext_gamma * vr.real_value;

vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;

}

vr += (0.5 * dt) * f_ext / ext_mass;

xr += dt * vr;

xr.apply_constraints();

if (this->b_periodic) this->wrap(xr);

}

// Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use

if (tasks[task_lower_wall] || tasks[task_upper_wall]) {

// Wall force

colvarvalue fw(x);

Line 1178

Line 1143

// For a periodic colvar, both walls may be applicable at the same time

// in which case we pick the closer one

if ( (!tasks[task_upper_wall]) ||

if ( (!is_enabled(f_cv_upper_wall)) ||

(this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {

cvm::real const grad = this->dist2_lgrad(x, lower_wall);

Line 1203

Line 1168

}

if (is_enabled(f_cv_extended_Lagrangian)) {

cvm::real dt = cvm::dt();

cvm::real f_ext;

// the total force is applied to the fictitious mass, while the

// atoms only feel the harmonic force

// fr: bias force on extended variable (without harmonic spring), for output in trajectory

// f_ext: total force on extended variable (including harmonic spring)

// f: - initially, external biasing force

// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force

fr = f;

f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);

f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);

// The total force acting on the extended variable is f_ext

// This will be used in the next timestep

ft_reported = f_ext;

// leapfrog: starting from x_i, f_i, v_(i-1/2)

vr += (0.5 * dt) * f_ext / ext_mass;

// Because of leapfrog, kinetic energy at time i is approximate

kinetic_energy = 0.5 * ext_mass * vr * vr;

potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);

// leap to v_(i+1/2)

if (is_enabled(f_cv_Langevin)) {

vr -= dt * ext_gamma * vr.real_value;

vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;

}

vr += (0.5 * dt) * f_ext / ext_mass;

xr += dt * vr;

xr.apply_constraints();

if (this->b_periodic) this->wrap(xr);

}

f_accumulated += f;

if (tasks[task_fdiff_velocity]) {

if (is_enabled(f_cv_fdiff_velocity)) {

// set it for the next step

x_old = x;

}

if (is_enabled(f_cv_subtract_applied_force)) {

f_old = f;

}

if (cvm::debug())

cvm::log("Done updating colvar \""+this->name+"\".\n");

return (potential_energy + kinetic_energy);

Line 1221

Line 1226

if (cvm::debug())

cvm::log("Communicating forces from colvar \""+this->name+"\".\n");

if (tasks[task_scripted]) {

if (is_enabled(f_cv_scripted)) {

std::vector<cvm::matrix2d<cvm::real> > func_grads;

func_grads.reserve(cvcs.size());

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),

cvcs[i]->value().size()));

}

Line 1232

Line 1238

if (res != COLVARS_OK) {

if (res == COLVARS_NOT_IMPLEMENTED) {

cvm::error("Colvar gradient scripts are not implemented.");

cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);

} else {

cvm::error("Error running colvar gradient script");

}

Line 1241

Line 1247

int grad_index = 0; // index in the scripted gradients, to account for some components being disabled

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

cvm::increase_depth();

// cvc force is colvar force times colvar/cvc Jacobian

// (vector-matrix product)

(cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],

(cvcs[i])->apply_force(colvarvalue(f_accumulated.as_vector() * func_grads[grad_index++],

cvcs[i]->value().type()));

cvm::decrease_depth();

}

} else if (x.type() == colvarvalue::type_scalar) {

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

cvm::increase_depth();

(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff *

(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *

cvm::real((cvcs[i])->sup_np) *

(std::pow((cvcs[i])->value().real_value,

(cvcs[i])->sup_np-1)) );

cvm::decrease_depth();

}

} else {

for (i = 0; i < cvcs.size(); i++) {

if (!cvcs[i]->b_enabled) continue;

if (!cvcs[i]->is_enabled()) continue;

cvm::increase_depth();

(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff);

(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);

cvm::decrease_depth();

}

// Accumulated forces have been applied, impulse-style

// Reset to start accumulating again

f_accumulated.reset();

if (cvm::debug())

cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");

}

int colvar::set_cvc_flags(std::vector<bool> const &flags) {

int colvar::set_cvc_flags(std::vector<bool> const &flags)

{

if (flags.size() != cvcs.size()) {

cvm::error("ERROR: Wrong number of CVC flags provided.");

return COLVARS_ERROR;

Line 1289

Line 1293

}

int colvar::update_cvc_flags()

{

// Update the enabled/disabled status of cvcs if necessary

if (cvc_flags.size()) {

n_active_cvcs = 0;

active_cvc_square_norm = 0.;

for (size_t i = 0; i < cvcs.size(); i++) {

cvcs[i]->feature_states[f_cvc_active]->enabled = cvc_flags[i];

if (cvcs[i]->is_enabled()) {

n_active_cvcs++;

active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;

}

if (!n_active_cvcs) {

cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");

return COLVARS_ERROR;

}

cvc_flags.resize(0);

}

return COLVARS_OK;

}

// ******************** METRIC FUNCTIONS ********************

// Use the metrics defined by \link cvc \endlink objects

bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const

{

if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) {

if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {

cvm::log("Error: requesting to histogram the "

cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "

"collective variable \""+this->name+"\", but a "

"requires lower and upper boundaries to be defined.\n");

"pair of lower and upper boundaries must be "

"defined.\n");

cvm::set_error_bits(INPUT_ERROR);

}

Line 1315

Line 1342

bool colvar::periodic_boundaries() const

{

if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) {

if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {

cvm::error("Error: requesting to histogram the "

cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "

"collective variable \""+this->name+"\", but a "

"requires lower and upper boundaries to be defined.\n");

"pair of lower and upper boundaries must be "

"defined.\n");

}

return periodic_boundaries(lower_boundary, upper_boundary);

Line 1329

Line 1354

cvm::real colvar::dist2(colvarvalue const &x1,

colvarvalue const &x2) const

{

if (b_homogeneous) {

if (is_enabled(f_cv_homogeneous)) {

return (cvcs[0])->dist2(x1, x2);

} else {

return x1.dist2(x2);

Line 1339

Line 1364

colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,

colvarvalue const &x2) const

{

if (b_homogeneous) {

if (is_enabled(f_cv_homogeneous)) {

return (cvcs[0])->dist2_lgrad(x1, x2);

} else {

return x1.dist2_grad(x2);

Line 1349

Line 1374

colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,

colvarvalue const &x2) const

{

if (b_homogeneous) {

if (is_enabled(f_cv_homogeneous)) {

return (cvcs[0])->dist2_rgrad(x1, x2);

} else {

return x2.dist2_grad(x1);

Line 1358

Line 1383

void colvar::wrap(colvarvalue &x) const

{

if (b_homogeneous) {

if (is_enabled(f_cv_homogeneous)) {

(cvcs[0])->wrap(x);

}

return;

Line 1404

Line 1429

cvm::to_str(x)+"\n");

}

if (tasks[colvar::task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

if ( !(get_keyval(conf, "extended_x", xr,

colvarvalue(x.type()), colvarparse::parse_silent)) &&

Line 1414

Line 1439

"\"extended_x\" or \"extended_v\" for the colvar \""+

name+"\", but you requested \"extendedLagrangian\".\n");

}

if (tasks[task_extended_lagrangian]) {

x_reported = xr;

} else {

x_reported = x;

}

if (tasks[task_output_velocity]) {

if (is_enabled(f_cv_output_velocity)) {

if ( !(get_keyval(conf, "v", v_fdiff,

colvarvalue(x.type()), colvarparse::parse_silent)) ) {

Line 1431

Line 1453

name+"\", but you requested \"outputVelocity\".\n");

}

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

v_reported = vr;

} else {

v_reported = v_fdiff;

Line 1446

Line 1468

{

size_t const start_pos = is.tellg();

if (tasks[task_output_value]) {

if (is_enabled(f_cv_output_value)) {

if (!(is >> x)) {

cvm::log("Error: in reading the value of colvar \""+

Line 1457

Line 1479

return is;

}

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

is >> xr;

x_reported = xr;

} else {

Line 1465

Line 1487

}

if (tasks[task_output_velocity]) {

if (is_enabled(f_cv_output_velocity)) {

is >> v_fdiff;

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

is >> vr;

v_reported = vr;

} else {

Line 1477

Line 1499

}

if (tasks[task_output_system_force]) {

if (is_enabled(f_cv_output_total_force)) {

is >> ft;

ft_reported = ft;

}

if (tasks[task_output_applied_force]) {

if (is_enabled(f_cv_output_applied_force)) {

is >> f;

}

Line 1501

Line 1523

<< std::setw(cvm::cv_width)

<< x << "\n";

if (tasks[task_output_velocity]) {

if (is_enabled(f_cv_output_velocity)) {

os << " v "

<< std::setprecision(cvm::cv_prec)

<< std::setw(cvm::cv_width)

<< v_reported << "\n";

}

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

os << " extended_x "

<< std::setprecision(cvm::cv_prec)

<< std::setw(cvm::cv_width)

Line 1531

Line 1553

os << " ";

if (tasks[task_output_value]) {

if (is_enabled(f_cv_output_value)) {

os << " "

<< cvm::wrap_string(this->name, this_cv_width);

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

// extended DOF

os << " r_"

<< cvm::wrap_string(this->name, this_cv_width-2);

}

if (tasks[task_output_velocity]) {

if (is_enabled(f_cv_output_velocity)) {

os << " v_"

<< cvm::wrap_string(this->name, this_cv_width-2);

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

// extended DOF

os << " vr_"

<< cvm::wrap_string(this->name, this_cv_width-3);

}

if (tasks[task_output_energy]) {

if (is_enabled(f_cv_output_energy)) {

os << " Ep_"

<< cvm::wrap_string(this->name, this_cv_width-3)

<< " Ek_"

<< cvm::wrap_string(this->name, this_cv_width-3);

}

if (tasks[task_output_system_force]) {

if (is_enabled(f_cv_output_total_force)) {

os << " fs_"

os << " ft_"

<< cvm::wrap_string(this->name, this_cv_width-3);

}

if (tasks[task_output_applied_force]) {

if (is_enabled(f_cv_output_applied_force)) {

os << " fa_"

<< cvm::wrap_string(this->name, this_cv_width-3);

}

Line 1580

Line 1602

{

os << " ";

if (tasks[task_output_value]) {

if (is_enabled(f_cv_output_value)) {

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< x;

Line 1593

Line 1615

<< x_reported;

}

if (tasks[task_output_velocity]) {

if (is_enabled(f_cv_output_velocity)) {

if (tasks[task_extended_lagrangian]) {

if (is_enabled(f_cv_extended_Lagrangian)) {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< v_fdiff;

Line 1606

Line 1628

<< v_reported;

}

if (tasks[task_output_energy]) {

if (is_enabled(f_cv_output_energy)) {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< potential_energy

Line 1615

Line 1637

}

if (tasks[task_output_system_force]) {

if (is_enabled(f_cv_output_total_force)) {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< ft_reported;

}

if (tasks[task_output_applied_force]) {

if (is_enabled(f_cv_output_applied_force)) {

if (tasks[task_extended_lagrangian]) {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< fr;

} else {

os << " "

<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)

<< f;

<< applied_force();

}

return os;

Line 1643

Line 1659

if (acf.size()) {

cvm::log("Writing acf to file \""+acf_outfile+"\".\n");

cvm::backup_file(acf_outfile.c_str());

cvm::ofstream acf_os(acf_outfile.c_str());

if (! acf_os.good()) {

if (! acf_os.is_open()) {

cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR);

}

write_acf(acf_os);

acf_os.close();

}

if (runave_os.good()) {

if (runave_os.is_open()) {

runave_os.close();

}

Line 1664

Line 1681

void colvar::analyze()

{

if (tasks[task_runave]) {

if (is_enabled(f_cv_runave)) {

calc_runave();

}

if (tasks[task_corrfunc]) {

if (is_enabled(f_cv_corrfunc)) {

calc_acf();

}

Line 1753

Line 1770

case acf_vel:

if (tasks[task_fdiff_velocity]) {

if (is_enabled(f_cv_fdiff_velocity)) {

// calc() should do this already, but this only happens in a

// simulation; better do it again in case a trajectory is

// being read

Line 1786

Line 1803

}

if (tasks[task_fdiff_velocity]) {

if (is_enabled(f_cv_fdiff_velocity)) {

// set it for the next step

x_old = x;

}

Line 1944

Line 1961

}

// Static members

std::vector<colvardeps::feature *> colvar::cv_features;

Legend:

Removed in v.1.25
changed lines
	Added in v.1.26

Made by using version 1.53 of cvs2html