Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76696986
colvar.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 9, 22:14
Size
60 KB
Mime Type
text/x-c
Expires
Sun, Aug 11, 22:14 (2 d)
Engine
blob
Format
Raw Data
Handle
19752440
Attached To
rLAMMPS lammps
colvar.cpp
View Options
// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvarscript.h"
#include <algorithm>
/// Compare two cvcs using their names
/// Used to sort CVC array in scripted coordinates
bool compare(colvar::cvc *i, colvar::cvc *j) {
return i->name < j->name;
}
colvar::colvar(std::string const &conf)
: colvarparse(conf)
{
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)));
if (cvm::colvar_by_name(this->name) != NULL) {
cvm::error("Error: this colvar cannot have the same name, \""+this->name+
"\", as another colvar.\n",
INPUT_ERROR);
return;
}
// Initialize dependency members
// Could be a function defined in a different source file, for space?
this->description = "colvar " + this->name;
// Initialize static array once and for all
init_cv_requires();
kinetic_energy = 0.0;
potential_energy = 0.0;
cvm::combine_errors(error_code, init_components(conf));
size_t i;
// Setup colvar as scripted function of components
if (get_keyval(conf, "scriptedFunction", scripted_function,
"", colvarparse::parse_silent)) {
// 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;
get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
x.type(colvarvalue::type_notset);
int t;
for (t = 0; t < colvarvalue::type_all; t++) {
if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
x.type(colvarvalue::Type(t));
break;
}
}
if (x.type() == colvarvalue::type_notset) {
cvm::error("Could not parse scripted colvar type.");
return;
}
cvm::log(std::string("Expecting colvar value of type ")
+ colvarvalue::type_desc(x.type()));
if (x.type() == colvarvalue::type_vector) {
int size;
if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
cvm::error("Error: no size specified for vector scripted function.");
return;
}
x.vector1d_value.resize(size);
}
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);
if(cvcs.size() > 1) {
cvm::log("Sorted list of components for this scripted colvar:");
for (i = 0; i < cvcs.size(); i++) {
cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name);
}
}
// Build ordered list of component values that will be
// passed to the script
for (i = 0; i < cvcs.size(); i++) {
sorted_cvc_values.push_back(&(cvcs[i]->value()));
}
}
if (!is_enabled(f_cv_scripted)) {
colvarvalue const &cvc_value = (cvcs[0])->value();
if (cvm::debug())
cvm::log ("This collective variable is a "+
colvarvalue::type_desc(cvc_value.type())+
((cvc_value.size() > 1) ? " with "+
cvm::to_str(cvc_value.size())+" individual components.\n" :
".\n"));
x.type(cvc_value);
x_reported.type(cvc_value);
}
// If using scripted biases, any colvar may receive bias forces
// and will need its gradient
if (cvm::scripted_forces()) {
enable(f_cv_gradient);
}
// check for linear combinations
{
bool lin = !is_enabled(f_cv_scripted);
for (i = 0; i < cvcs.size(); i++) {
// 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() && lin)
cvm::log("Warning: You are using a non-linear polynomial "
"combination to define this collective variable, "
"some biasing methods may be unavailable.\n");
lin = false;
if ((cvcs[i])->sup_np < 0) {
cvm::log("Warning: you chose a negative exponent in the combination; "
"if you apply forces, the simulation may become unstable "
"when the component \""+
(cvcs[i])->function_type+"\" approaches zero.\n");
}
}
}
feature_states[f_cv_linear]->enabled = lin;
}
// Colvar is homogeneous iff:
// - it is linear (hence not scripted)
// - all cvcs have coefficient 1 or -1
// i.e. sum or difference of cvcs
{
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) {
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 && 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) {
b_periodic = false;
period = 0.0;
}
}
feature_states[f_cv_periodic]->enabled = b_periodic;
// check that cvcs are compatible
for (i = 0; i < cvcs.size(); i++) {
if ((cvcs[i])->b_periodic && !b_periodic) {
cvm::log("Warning: although this component is periodic, the colvar will "
"not be treated as periodic, either because the exponent is not "
"1, or because multiple components are present. Make sure that "
"you know what you are doing!");
}
// components may have different types only for scripted functions
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. "
"You must use the same type in order to "
" sum them together.\n", INPUT_ERROR);
return;
}
}
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);
if (width <= 0.0) {
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());
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)) {
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(f_cv_lower_wall);
}
if (get_keyval(conf, "upperBoundary", upper_boundary, 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(f_cv_upper_wall);
}
}
if (is_enabled(f_cv_lower_boundary)) {
get_keyval(conf, "hardLowerBoundary", hard_lower_boundary, false);
}
if (is_enabled(f_cv_upper_boundary)) {
get_keyval(conf, "hardUpperBoundary", hard_upper_boundary, false);
}
// consistency checks for boundaries and walls
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)+
", is not higher than the lower boundary, "+
cvm::to_str(lower_boundary)+".\n",
INPUT_ERROR);
}
}
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)+
", is not higher than the lower wall, "+
cvm::to_str(lower_wall)+".\n",
INPUT_ERROR);
}
}
get_keyval(conf, "expandBoundaries", expand_boundaries, false);
if (expand_boundaries && periodic_boundaries()) {
cvm::error("Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.\n",
INPUT_ERROR);
}
if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
cvm::error("Error: inconsistent configuration "
"(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).\n",
INPUT_ERROR);
}
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
{
bool b_extended_Lagrangian;
get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
if (b_extended_Lagrangian) {
cvm::real temp, tolerance, period;
cvm::log("Enabling the extended Lagrangian term for colvar \""+
this->name+"\".\n");
// 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());
fr.type(value());
const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
if (temp <= 0.0) {
if (found)
cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
else
cvm::error("Error: a positive temperature must be provided, either "
"by enabling a thermostat, or through \"extendedTemp\".\n",
INPUT_ERROR);
}
get_keyval(conf, "extendedFluctuation", tolerance);
if (tolerance <= 0.0) {
cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
}
ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
get_keyval(conf, "extendedTimeConstant", period, 200.0);
if (period <= 0.0) {
cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
}
ext_mass = (cvm::boltzmann() * temp * period * period)
/ (4.0 * PI * PI * tolerance * tolerance);
cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");
{
bool b_output_energy;
get_keyval(conf, "outputEnergy", b_output_energy, false);
if (b_output_energy) {
enable(f_cv_output_energy);
}
}
get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
if (ext_gamma < 0.0) {
cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
}
if (ext_gamma != 0.0) {
enable(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());
}
}
}
{
bool b_output_value;
get_keyval(conf, "outputValue", b_output_value, true);
if (b_output_value) {
enable(f_cv_output_value);
}
}
{
bool b_output_velocity;
get_keyval(conf, "outputVelocity", b_output_velocity, false);
if (b_output_velocity) {
enable(f_cv_output_velocity);
}
}
{
bool b_output_system_force;
get_keyval(conf, "outputSystemForce", b_output_system_force, false);
if (b_output_system_force) {
enable(f_cv_output_system_force);
}
}
{
bool b_output_applied_force;
get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);
if (b_output_applied_force) {
enable(f_cv_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());
if (cvm::b_analysis)
parse_analysis(conf);
if (cvm::debug())
cvm::log("Done initializing collective variable \""+this->name+"\".\n");
}
// 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;
cvm::combine_errors(error_code,
init_components_type<distance>(conf,
"distance", "distance"));
cvm::combine_errors(error_code,
init_components_type<distance_vec>(conf,
"distance vector", "distanceVec"));
cvm::combine_errors(error_code,
init_components_type<cartesian>(conf,
"Cartesian coordinates", "cartesian"));
cvm::combine_errors(error_code,
init_components_type<distance_dir>(conf,
"distance vector "
"direction", "distanceDir"));
cvm::combine_errors(error_code,
init_components_type<distance_z>(conf,
"distance projection "
"on an axis", "distanceZ"));
cvm::combine_errors(error_code,
init_components_type<distance_xy>(conf,
"distance projection "
"on a plane", "distanceXY"));
cvm::combine_errors(error_code,
init_components_type<distance_inv>(conf,
"average distance "
"weighted by inverse power",
"distanceInv"));
cvm::combine_errors(error_code,
init_components_type<distance_pairs>(conf,
"N1xN2-long vector "
"of pairwise distances",
"distancePairs"));
cvm::combine_errors(error_code,
init_components_type<coordnum>(conf,
"coordination "
"number", "coordNum"));
cvm::combine_errors(error_code,
init_components_type<selfcoordnum>(conf,
"self-coordination "
"number", "selfCoordNum"));
cvm::combine_errors(error_code,
init_components_type<angle>(conf,
"angle", "angle"));
cvm::combine_errors(error_code,
init_components_type<dipole_angle>(conf,
"dipole angle", "dipoleAngle"));
cvm::combine_errors(error_code,
init_components_type<dihedral>(conf,
"dihedral", "dihedral"));
cvm::combine_errors(error_code,
init_components_type<h_bond>(conf,
"hydrogen bond", "hBond"));
// cvm::combine_errors(error_code, init_components_type<alpha_dihedrals>(conf, "alpha helix", "alphaDihedrals"));
cvm::combine_errors(error_code,
init_components_type<alpha_angles>(conf,
"alpha helix", "alpha"));
cvm::combine_errors(error_code,
init_components_type<dihedPC>(conf,
"dihedral "
"principal component", "dihedralPC"));
cvm::combine_errors(error_code,
init_components_type<orientation>(conf,
"orientation", "orientation"));
cvm::combine_errors(error_code,
init_components_type<orientation_angle>(conf,
"orientation "
"angle", "orientationAngle"));
cvm::combine_errors(error_code,
init_components_type<orientation_proj>(conf,
"orientation "
"projection", "orientationProj"));
cvm::combine_errors(error_code,
init_components_type<tilt>(conf,
"tilt", "tilt"));
cvm::combine_errors(error_code,
init_components_type<spin_angle>(conf,
"spin angle", "spinAngle"));
cvm::combine_errors(error_code,
init_components_type<rmsd>(conf,
"RMSD", "rmsd"));
// cvm::combine_errors(error_code, init_components_type <logmsd>(conf,"logarithm of MSD", "logmsd"));
cvm::combine_errors(error_code,
init_components_type<gyration>(conf,
"radius of "
"gyration", "gyration"));
cvm::combine_errors(error_code,
init_components_type<inertia>(conf,
"moment of "
"inertia", "inertia"));
cvm::combine_errors(error_code,
init_components_type<inertia_z>(conf,
"moment of inertia around an axis",
"inertiaZ"));
cvm::combine_errors(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_system_force_calc)) {
cvm::request_system_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
std::list<int> temp_id_list;
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
for (size_t k = 0; k < ag.size(); k++) {
temp_id_list.push_back(ag[k].id);
}
}
}
temp_id_list.sort();
temp_id_list.unique();
// atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
unsigned int id_i = 0;
std::list<int>::iterator li;
for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
atom_ids[id_i] = *li;
id_i++;
}
temp_id_list.clear();
atomic_gradients.resize(atom_ids.size());
if (atom_ids.size()) {
if (cvm::debug())
cvm::log("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
} else {
cvm::log("Warning: colvar components communicated no atom IDs.\n");
}
}
int colvar::parse_analysis(std::string const &conf)
{
// if (cvm::debug())
// cvm::log ("Parsing analysis flags for collective variable \""+
// this->name+"\".\n");
runave_length = 0;
bool b_runave = false;
if (get_keyval(conf, "runAve", b_runave) && b_runave) {
enable(f_cv_runave);
get_keyval(conf, "runAveLength", runave_length, 1000);
get_keyval(conf, "runAveStride", runave_stride, 1);
if ((cvm::restart_out_freq % runave_stride) != 0) {
cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
}
std::string runave_outfile;
get_keyval(conf, "runAveOutputFile", runave_outfile,
std::string(cvm::output_prefix+"."+
this->name+".runave.traj"));
size_t const this_cv_width = x.output_width(cvm::cv_width);
cvm::backup_file(runave_outfile.c_str());
runave_os.open(runave_outfile.c_str());
runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
<< " "
<< cvm::wrap_string("running average", this_cv_width)
<< " "
<< cvm::wrap_string("running stddev", this_cv_width)
<< "\n";
}
acf_length = 0;
bool b_acf = false;
if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {
enable(f_cv_corrfunc);
std::string acf_colvar_name;
get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
if (acf_colvar_name == this->name) {
cvm::log("Calculating auto-correlation function.\n");
} else {
cvm::log("Calculating correlation function with \""+
this->name+"\".\n");
}
std::string acf_type_str;
get_keyval(conf, "corrFuncType", acf_type_str, to_lower_cppstr(std::string("velocity")));
if (acf_type_str == to_lower_cppstr(std::string("coordinate"))) {
acf_type = acf_coor;
} else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
acf_type = acf_vel;
enable(f_cv_fdiff_velocity);
if (acf_colvar_name.size())
(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 {
cvm::log("Unknown type of correlation function, \""+
acf_type_str+"\".\n");
cvm::set_error_bit(INPUT_ERROR);
}
get_keyval(conf, "corrFuncOffset", acf_offset, 0);
get_keyval(conf, "corrFuncLength", acf_length, 1000);
get_keyval(conf, "corrFuncStride", acf_stride, 1);
if ((cvm::restart_out_freq % acf_stride) != 0) {
cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
}
get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
get_keyval(conf, "corrFuncOutputFile", acf_outfile,
std::string(cvm::output_prefix+"."+this->name+
".corrfunc.dat"));
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void colvar::setup() {
// loop over all components to reset masses of all groups
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
atoms.setup();
atoms.reset_mass(name,i,ig);
atoms.read_positions();
}
}
}
colvar::~colvar()
{
// Clear references to this colvar's cvcs as children
// 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
for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
cvi != cvm::colvars.end();
++cvi) {
if ( *cvi == this) {
cvm::colvars.erase(cvi);
break;
}
}
}
// ******************** CALC FUNCTIONS ********************
// Default schedule (everything is serialized)
int colvar::calc()
{
// Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
int error_code = COLVARS_OK;
if (is_enabled(f_cv_active)) {
cvm::combine_errors(error_code, update_cvc_flags());
cvm::combine_errors(error_code, calc_cvcs());
cvm::combine_errors(error_code, collect_cvc_data());
}
return error_code;
}
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\", components "+
cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
cvm::combine_errors(error_code, check_cvc_range(first_cvc, num_cvcs));
if (error_code != COLVARS_OK) {
return error_code;
}
cvm::combine_errors(error_code, calc_cvc_values(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_gradients(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_sys_forces(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_Jacobians(first_cvc, num_cvcs));
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\".\n");
return error_code;
}
int colvar::collect_cvc_data()
{
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
int error_code = COLVARS_OK;
cvm::combine_errors(error_code, collect_cvc_values());
cvm::combine_errors(error_code, collect_cvc_gradients());
cvm::combine_errors(error_code, collect_cvc_sys_forces());
cvm::combine_errors(error_code, collect_cvc_Jacobians());
cvm::combine_errors(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");
// First, calculate component values
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();
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;
// combine them appropriately, using either a scripted function or a polynomial
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 COLVARS_NOT_IMPLEMENTED;
}
if (res != COLVARS_OK) {
cvm::error("Error running scripted colvar");
return COLVARS_OK;
}
} else if (x.type() == colvarvalue::type_scalar) {
// polynomial combination allowed
for (i = 0; i < cvcs.size(); i++) {
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 {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
}
}
if (cvm::debug())
cvm::log("Colvar \""+this->name+"\" has value "+
cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
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 (is_enabled(f_cv_gradient)) {
if (cvm::debug())
cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
// calculate the gradients of each component
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])->calc_gradients();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
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();
if (cvm::debug())
cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
}
return COLVARS_OK;
}
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.", COLVARS_NOT_IMPLEMENTED);
return COLVARS_NOT_IMPLEMENTED;
}
// Collect the atomic gradients inside colvar object
for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
atomic_gradients[a].reset();
}
for (i = 0; i < cvcs.size(); i++) {
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 (ag.b_rotate) {
cvm::rotation const rot_inv = ag.rot.inverse();
for (size_t k = 0; k < ag.size(); k++) {
size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
ag[k].id) - atom_ids.begin();
atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
}
} else {
for (size_t k = 0; k < ag.size(); k++) {
size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
ag[k].id) - atom_ids.begin();
atomic_gradients[a] += coeff * ag[k].grad;
}
}
}
}
}
return COLVARS_OK;
}
int colvar::calc_cvc_sys_forces(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 (is_enabled(f_cv_system_force_calc)) {
if (cvm::debug())
cvm::log("Calculating system force of colvar \""+this->name+"\".\n");
// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
// Disabled check to allow for explicit system force calculation
// even with extended Lagrangian
if (cvm::step_relative() > 0) {
cvm::increase_depth();
// get from the cvcs the system forces from the PREVIOUS step
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->is_enabled()) continue;
cvc_count++;
(cvcs[i])->calc_force_invgrads();
}
cvm::decrease_depth();
}
if (cvm::debug())
cvm::log("Done calculating system force of colvar \""+this->name+"\".\n");
}
return COLVARS_OK;
}
int colvar::collect_cvc_sys_forces()
{
if (is_enabled(f_cv_system_force_calc)) {
ft.reset();
if (cvm::step_relative() > 0) {
// get from the cvcs the system 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])->system_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
}
}
if (!is_enabled(f_cv_hide_Jacobian)) {
// add the Jacobian force to the system force, and don't apply any silent
// correction internally: biases such as colvarbias_abf will handle it
ft += fj;
}
}
return COLVARS_OK;
}
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;
else {
v_fdiff = fdiff_velocity(x_old, x);
v_reported = v_fdiff;
}
}
if (is_enabled(f_cv_extended_Lagrangian)) {
// initialize the restraint center in the first step to the value
// just calculated from the cvcs
// TODO: put it in the restart information
if (cvm::step_relative() == 0) {
xr = x;
vr = 0.0; // (already 0; added for clarity)
}
// report the restraint center as "value"
x_reported = xr;
v_reported = vr;
// the "system force" with the extended Lagrangian is just the
// harmonic term acting on the extended coordinate
// Note: this is the force for current timestep
ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
} else {
x_reported = x;
ft_reported = ft;
}
return COLVARS_OK;
}
cvm::real colvar::update_forces_energy()
{
if (cvm::debug())
cvm::log("Updating colvar \""+this->name+"\".\n");
// set to zero the applied force
f.type(value());
f.reset();
// add the biases' force, which at this point should already have
// been summed over each bias using this colvar
f += fb;
if (is_enabled(f_cv_Jacobian)) {
// the instantaneous Jacobian force was not included in the reported system force;
// instead, it is subtracted from the applied force (silent Jacobian correction)
if (is_enabled(f_cv_hide_Jacobian))
f -= fj;
}
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 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 (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);
}
// Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use
if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
// Wall force
colvarvalue fw(x);
fw.reset();
if (cvm::debug())
cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
if ( (!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);
if (grad < 0.0) {
fw = -0.5 * lower_wall_k * grad;
f += fw;
if (cvm::debug())
cvm::log("Applying a lower wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
}
} else {
cvm::real const grad = this->dist2_lgrad(x, upper_wall);
if (grad > 0.0) {
fw = -0.5 * upper_wall_k * grad;
f += fw;
if (cvm::debug())
cvm::log("Applying an upper wall force("+
cvm::to_str(fw)+") to \""+this->name+"\".\n");
}
}
}
f_accumulated += f;
if (is_enabled(f_cv_fdiff_velocity)) {
// set it for the next step
x_old = x;
}
if (cvm::debug())
cvm::log("Done updating colvar \""+this->name+"\".\n");
return (potential_energy + kinetic_energy);
}
void colvar::communicate_forces()
{
size_t i;
if (cvm::debug())
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
if (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]->is_enabled()) continue;
func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
cvcs[i]->value().size()));
}
int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
if (res != COLVARS_OK) {
if (res == COLVARS_NOT_IMPLEMENTED) {
cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);
} else {
cvm::error("Error running colvar gradient script");
}
return;
}
int grad_index = 0; // index in the scripted gradients, to account for some components being disabled
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
// cvc force is colvar force times colvar/cvc Jacobian
// (vector-matrix product)
(cvcs[i])->apply_force(colvarvalue(f_accumulated.as_vector() * func_grads[grad_index++],
cvcs[i]->value().type()));
}
} else if (x.type() == colvarvalue::type_scalar) {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff *
cvm::real((cvcs[i])->sup_np) *
(std::pow((cvcs[i])->value().real_value,
(cvcs[i])->sup_np-1)) );
}
} else {
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff);
}
}
// 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)
{
if (flags.size() != cvcs.size()) {
cvm::error("ERROR: Wrong number of CVC flags provided.");
return COLVARS_ERROR;
}
// We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
// so we store the flags that will be enforced at the next call to calc()
cvc_flags = flags;
return COLVARS_OK;
}
int colvar::update_cvc_flags()
{
// 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 ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
"requires lower and upper boundaries to be defined.\n");
cvm::set_error_bit(INPUT_ERROR);
}
if (period > 0.0) {
if ( ((std::sqrt(this->dist2(lb, ub))) / this->width)
< 1.0E-10 ) {
return true;
}
}
return false;
}
bool colvar::periodic_boundaries() const
{
if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
"requires lower and upper boundaries to be defined.\n");
}
return periodic_boundaries(lower_boundary, upper_boundary);
}
cvm::real colvar::dist2(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (is_enabled(f_cv_homogeneous)) {
return (cvcs[0])->dist2(x1, x2);
} else {
return x1.dist2(x2);
}
}
colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (is_enabled(f_cv_homogeneous)) {
return (cvcs[0])->dist2_lgrad(x1, x2);
} else {
return x1.dist2_grad(x2);
}
}
colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const
{
if (is_enabled(f_cv_homogeneous)) {
return (cvcs[0])->dist2_rgrad(x1, x2);
} else {
return x2.dist2_grad(x1);
}
}
void colvar::wrap(colvarvalue &x) const
{
if (is_enabled(f_cv_homogeneous)) {
(cvcs[0])->wrap(x);
}
return;
}
// ******************** INPUT FUNCTIONS ********************
std::istream & colvar::read_restart(std::istream &is)
{
size_t const start_pos = is.tellg();
std::string conf;
if ( !(is >> colvarparse::read_block("colvar", conf)) ) {
// this is not a colvar block
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
{
std::string check_name = "";
if ( (get_keyval(conf, "name", check_name,
std::string(""), colvarparse::parse_silent)) &&
(check_name != name) ) {
cvm::error("Error: the state file does not match the "
"configuration file, at colvar \""+name+"\".\n");
}
if (check_name.size() == 0) {
cvm::error("Error: Collective variable in the "
"restart file without any identifier.\n");
}
}
if ( !(get_keyval(conf, "x", x,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"the value of the colvar \""+
name+"\" .\n");
} else {
cvm::log("Restarting collective variable \""+name+"\" from value: "+
cvm::to_str(x)+"\n");
}
if (is_enabled(f_cv_extended_Lagrangian)) {
if ( !(get_keyval(conf, "extended_x", xr,
colvarvalue(x.type()), colvarparse::parse_silent)) &&
!(get_keyval(conf, "extended_v", vr,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"\"extended_x\" or \"extended_v\" for the colvar \""+
name+"\", but you requested \"extendedLagrangian\".\n");
}
x_reported = xr;
} else {
x_reported = x;
}
if (is_enabled(f_cv_output_velocity)) {
if ( !(get_keyval(conf, "v", v_fdiff,
colvarvalue(x.type()), colvarparse::parse_silent)) ) {
cvm::log("Error: restart file does not contain "
"the velocity for the colvar \""+
name+"\", but you requested \"outputVelocity\".\n");
}
if (is_enabled(f_cv_extended_Lagrangian)) {
v_reported = vr;
} else {
v_reported = v_fdiff;
}
}
return is;
}
std::istream & colvar::read_traj(std::istream &is)
{
size_t const start_pos = is.tellg();
if (is_enabled(f_cv_output_value)) {
if (!(is >> x)) {
cvm::log("Error: in reading the value of colvar \""+
this->name+"\" from trajectory.\n");
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
return is;
}
if (is_enabled(f_cv_extended_Lagrangian)) {
is >> xr;
x_reported = xr;
} else {
x_reported = x;
}
}
if (is_enabled(f_cv_output_velocity)) {
is >> v_fdiff;
if (is_enabled(f_cv_extended_Lagrangian)) {
is >> vr;
v_reported = vr;
} else {
v_reported = v_fdiff;
}
}
if (is_enabled(f_cv_output_system_force)) {
is >> ft;
ft_reported = ft;
}
if (is_enabled(f_cv_output_applied_force)) {
is >> f;
}
return is;
}
// ******************** OUTPUT FUNCTIONS ********************
std::ostream & colvar::write_restart(std::ostream &os) {
os << "colvar {\n"
<< " name " << name << "\n"
<< " x "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< x << "\n";
if (is_enabled(f_cv_output_velocity)) {
os << " v "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< v_reported << "\n";
}
if (is_enabled(f_cv_extended_Lagrangian)) {
os << " extended_x "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< xr << "\n"
<< " extended_v "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< vr << "\n";
}
os << "}\n\n";
return os;
}
std::ostream & colvar::write_traj_label(std::ostream & os)
{
size_t const this_cv_width = x.output_width(cvm::cv_width);
os << " ";
if (is_enabled(f_cv_output_value)) {
os << " "
<< cvm::wrap_string(this->name, this_cv_width);
if (is_enabled(f_cv_extended_Lagrangian)) {
// extended DOF
os << " r_"
<< cvm::wrap_string(this->name, this_cv_width-2);
}
}
if (is_enabled(f_cv_output_velocity)) {
os << " v_"
<< cvm::wrap_string(this->name, this_cv_width-2);
if (is_enabled(f_cv_extended_Lagrangian)) {
// extended DOF
os << " vr_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
}
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 (is_enabled(f_cv_output_system_force)) {
os << " fs_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
if (is_enabled(f_cv_output_applied_force)) {
os << " fa_"
<< cvm::wrap_string(this->name, this_cv_width-3);
}
return os;
}
std::ostream & colvar::write_traj(std::ostream &os)
{
os << " ";
if (is_enabled(f_cv_output_value)) {
if (is_enabled(f_cv_extended_Lagrangian)) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< x;
}
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< x_reported;
}
if (is_enabled(f_cv_output_velocity)) {
if (is_enabled(f_cv_extended_Lagrangian)) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< v_fdiff;
}
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< v_reported;
}
if (is_enabled(f_cv_output_energy)) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< potential_energy
<< " "
<< kinetic_energy;
}
if (is_enabled(f_cv_output_system_force)) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< ft_reported;
}
if (is_enabled(f_cv_output_applied_force)) {
if (is_enabled(f_cv_extended_Lagrangian)) {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< fr;
} else {
os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< f;
}
}
return os;
}
int colvar::write_output_files()
{
if (cvm::b_analysis) {
if (acf.size()) {
cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
cvm::backup_file(acf_outfile.c_str());
cvm::ofstream acf_os(acf_outfile.c_str());
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.is_open()) {
runave_os.close();
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
// ******************** ANALYSIS FUNCTIONS ********************
void colvar::analyze()
{
if (is_enabled(f_cv_runave)) {
calc_runave();
}
if (is_enabled(f_cv_corrfunc)) {
calc_acf();
}
}
inline void history_add_value(size_t const &history_length,
std::list<colvarvalue> &history,
colvarvalue const &new_value)
{
history.push_front(new_value);
if (history.size() > history_length)
history.pop_back();
}
inline void history_incr(std::list< std::list<colvarvalue> > &history,
std::list< std::list<colvarvalue> >::iterator &history_p)
{
if ((++history_p) == history.end())
history_p = history.begin();
}
int colvar::calc_acf()
{
// using here an acf_stride-long list of vectors for either
// coordinates(acf_x_history) or velocities (acf_v_history); each vector can
// contain up to acf_length values, which are contiguous in memory
// representation but separated by acf_stride in the time series;
// the pointer to each vector is changed at every step
if (acf_x_history.empty() && acf_v_history.empty()) {
// first-step operations
colvar *cfcv = (acf_colvar_name.size() ?
cvm::colvar_by_name(acf_colvar_name) :
this);
if (colvarvalue::check_types(cfcv->value(), value())) {
cvm::error("Error: correlation function between \""+cfcv->name+
"\" and \""+this->name+"\" cannot be calculated, "
"because their value types are different.\n",
INPUT_ERROR);
}
acf_nframes = 0;
cvm::log("Colvar \""+this->name+"\": initializing ACF calculation.\n");
if (acf.size() < acf_length+1)
acf.resize(acf_length+1, 0.0);
size_t i;
switch (acf_type) {
case acf_vel:
// allocate space for the velocities history
for (i = 0; i < acf_stride; i++) {
acf_v_history.push_back(std::list<colvarvalue>());
}
acf_v_history_p = acf_v_history.begin();
break;
case acf_coor:
case acf_p2coor:
// allocate space for the coordinates history
for (i = 0; i < acf_stride; i++) {
acf_x_history.push_back(std::list<colvarvalue>());
}
acf_x_history_p = acf_x_history.begin();
break;
default:
break;
}
} else {
colvar *cfcv = (acf_colvar_name.size() ?
cvm::colvar_by_name(acf_colvar_name) :
this);
switch (acf_type) {
case acf_vel:
if (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
v_reported = v_fdiff = fdiff_velocity(x_old, cfcv->value());
}
calc_vel_acf((*acf_v_history_p), cfcv->velocity());
// store this value in the history
history_add_value(acf_length+acf_offset, *acf_v_history_p, cfcv->velocity());
// if stride is larger than one, cycle among different histories
history_incr(acf_v_history, acf_v_history_p);
break;
case acf_coor:
calc_coor_acf((*acf_x_history_p), cfcv->value());
history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
history_incr(acf_x_history, acf_x_history_p);
break;
case acf_p2coor:
calc_p2coor_acf((*acf_x_history_p), cfcv->value());
history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
history_incr(acf_x_history, acf_x_history_p);
break;
default:
break;
}
}
if (is_enabled(f_cv_fdiff_velocity)) {
// set it for the next step
x_old = x;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
colvarvalue const &v)
{
// loop over stored velocities and add to the ACF, but only the
// length is sufficient to hold an entire row of ACF values
if (v_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator vs_i = v_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++vs_i;
// current vel with itself
*(acf_i) += v.norm2();
++acf_i;
// inner products of previous velocities with current (acf_i and
// vs_i are updated)
colvarvalue::inner_opt(v, vs_i, v_list.end(), acf_i);
acf_nframes++;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
colvarvalue const &x)
{
// same as above but for coordinates
if (x_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator xs_i = x_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++xs_i;
*(acf_i++) += x.norm2();
colvarvalue::inner_opt(x, xs_i, x_list.end(), acf_i);
acf_nframes++;
}
}
void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
colvarvalue const &x)
{
// same as above but with second order Legendre polynomial instead
// of just the scalar product
if (x_list.size() >= acf_length+acf_offset) {
std::list<colvarvalue>::iterator xs_i = x_list.begin();
std::vector<cvm::real>::iterator acf_i = acf.begin();
for (size_t i = 0; i < acf_offset; i++)
++xs_i;
// value of P2(0) = 1
*(acf_i++) += 1.0;
colvarvalue::p2leg_opt(x, xs_i, x_list.end(), acf_i);
acf_nframes++;
}
}
void colvar::write_acf(std::ostream &os)
{
if (!acf_nframes)
cvm::log("Warning: ACF was not calculated (insufficient frames).\n");
os.setf(std::ios::scientific, std::ios::floatfield);
os << "# Autocorrelation function for collective variable \""
<< this->name << "\"\n";
// one frame is used for normalization, the statistical sample is
// hence decreased
os << "# nframes = " << (acf_normalize ?
acf_nframes - 1 :
acf_nframes) << "\n";
cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
std::vector<cvm::real>::iterator acf_i;
size_t it = acf_offset;
for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
os << std::setw(cvm::it_width) << acf_stride * (it++) << " "
<< std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width)
<< ( acf_normalize ?
(*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
(*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
}
}
void colvar::calc_runave()
{
if (x_history.empty()) {
runave.type(value().type());
runave.reset();
// first-step operations
if (cvm::debug())
cvm::log("Colvar \""+this->name+
"\": initializing running average calculation.\n");
acf_nframes = 0;
x_history.push_back(std::list<colvarvalue>());
x_history_p = x_history.begin();
} else {
if ( (cvm::step_relative() % runave_stride) == 0) {
if ((*x_history_p).size() >= runave_length-1) {
runave = x;
std::list<colvarvalue>::iterator xs_i;
for (xs_i = (*x_history_p).begin();
xs_i != (*x_history_p).end(); ++xs_i) {
runave += (*xs_i);
}
runave *= 1.0 / cvm::real(runave_length);
runave.apply_constraints();
runave_variance = 0.0;
runave_variance += this->dist2(x, runave);
for (xs_i = (*x_history_p).begin();
xs_i != (*x_history_p).end(); ++xs_i) {
runave_variance += this->dist2(x, (*xs_i));
}
runave_variance *= 1.0 / cvm::real(runave_length-1);
runave_os << std::setw(cvm::it_width) << cvm::step_relative()
<< " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< runave << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< std::sqrt(runave_variance) << "\n";
}
history_add_value(runave_length, *x_history_p, x);
}
}
}
// Static members
std::vector<cvm::deps::feature *> colvar::cv_features;
Event Timeline
Log In to Comment