f3cf407a21
This commit includes several fixes to moving restraints; also added is support for runtime integration of 2D and 3D PMFs from ABF. Mostly changes to existing member functions, with few additions in classes not directly accessible by LAMMPS. Also removed are calls to std::pow(), replaced by a copy of MathSpecial::powint(). Relevant commits in Colvars repository: 7307b5c 2017-12-14 Doc improvements [Giacomo Fiorin] 7f86f37 2017-12-14 Allow K-changing restraints computing accumulated work; fix staged-k TI estimator [Giacomo Fiorin] 7c1c175 2017-12-14 Fix 1D ABF trying to do pABF [Jérôme Hénin] b94aa7e 2017-11-16 Unify PMF output for 1D, 2D and 3D in ABF [Jérôme Hénin] 771a88f 2017-11-15 Poisson integration for all BC in 2d and 3d [Jérôme Hénin] 6af4d60 2017-12-01 Print message when issuing cv delete in VMD [Giacomo Fiorin] 4413972 2017-11-30 Check for homogeneous colvar to set it periodic [Jérôme Hénin] 95fe4b2 2017-11-06 Allow abf_integrate to start in bin with 1 sample [Jérôme Hénin] 06eea27 2017-10-23 Shorten a few constructs by using the power function [Giacomo Fiorin] 3165dfb 2017-10-20 Move includes of colvarproxy.h from headers to files [Giacomo Fiorin] 32a867b 2017-10-20 Add optimized powint function from LAMMPS headers [Giacomo Fiorin] 3ad070a 2017-10-20 Remove some unused includes, isolate calls to std::pow() [Giacomo Fiorin] 0aaf540 2017-10-20 Replace all calls to std::pow() where the exponent is not an integer [Giacomo Fiorin]
2298 lines
69 KiB
C++
2298 lines
69 KiB
C++
// -*- c++ -*-
|
|
|
|
// This file is part of the Collective Variables module (Colvars).
|
|
// The original version of Colvars and its updates are located at:
|
|
// https://github.com/colvars/colvars
|
|
// Please update all Colvars source files before making any changes.
|
|
// If you wish to distribute your changes, please submit them to the
|
|
// Colvars repository at GitHub.
|
|
|
|
#include "colvarmodule.h"
|
|
#include "colvarvalue.h"
|
|
#include "colvarparse.h"
|
|
#include "colvar.h"
|
|
#include "colvarcomp.h"
|
|
#include "colvarscript.h"
|
|
|
|
// used in build_atom_list()
|
|
#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()
|
|
: prev_timestep(-1)
|
|
{
|
|
// Initialize static array once and for all
|
|
runave_os = NULL;
|
|
init_cv_requires();
|
|
}
|
|
|
|
|
|
int colvar::init(std::string const &conf)
|
|
{
|
|
cvm::log("Initializing a new collective variable.\n");
|
|
colvarparse::init(conf);
|
|
|
|
int error_code = COLVARS_OK;
|
|
|
|
colvarmodule *cv = cvm::main();
|
|
|
|
get_keyval(conf, "name", this->name,
|
|
(std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));
|
|
|
|
if ((cvm::colvar_by_name(this->name) != NULL) &&
|
|
(cvm::colvar_by_name(this->name) != this)) {
|
|
cvm::error("Error: this colvar cannot have the same name, \""+this->name+
|
|
"\", as another colvar.\n",
|
|
INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
// Initialize dependency members
|
|
// Could be a function defined in a different source file, for space?
|
|
|
|
this->description = "colvar " + this->name;
|
|
|
|
kinetic_energy = 0.0;
|
|
potential_energy = 0.0;
|
|
|
|
error_code |= init_components(conf);
|
|
if (error_code != COLVARS_OK) {
|
|
return cvm::get_error();
|
|
}
|
|
|
|
size_t i;
|
|
|
|
#ifdef LEPTON
|
|
error_code |= init_custom_function(conf);
|
|
if (error_code != COLVARS_OK) {
|
|
return cvm::get_error();
|
|
}
|
|
#endif
|
|
|
|
// Setup colvar as scripted function of components
|
|
if (get_keyval(conf, "scriptedFunction", scripted_function,
|
|
"", colvarparse::parse_silent)) {
|
|
|
|
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.", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
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.",
|
|
INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
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) || is_enabled(f_cv_custom_function))) {
|
|
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) || is_enabled(f_cv_custom_function));
|
|
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");
|
|
}
|
|
}
|
|
}
|
|
set_enabled(f_cv_linear, lin);
|
|
}
|
|
|
|
// Colvar is homogeneous if:
|
|
// - 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;
|
|
}
|
|
}
|
|
set_enabled(f_cv_homogeneous, homogeneous);
|
|
}
|
|
|
|
// Colvar is deemed periodic if:
|
|
// - it is homogeneous
|
|
// - all cvcs are periodic
|
|
// - all cvcs have the same period
|
|
if (is_enabled(f_cv_homogeneous) && cvcs[0]->b_periodic) { // TODO make this a CVC feature
|
|
bool b_periodic = true;
|
|
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;
|
|
cvm::log("Warning: although one component is periodic, this colvar will "
|
|
"not be treated as periodic, either because the exponent is not "
|
|
"1, or because components of different periodicity are defined. "
|
|
"Make sure that you know what you are doing!");
|
|
}
|
|
}
|
|
set_enabled(f_cv_periodic, b_periodic);
|
|
}
|
|
|
|
// check that cvcs are compatible
|
|
|
|
for (i = 0; i < cvcs.size(); i++) {
|
|
|
|
// components may have different types only for scripted functions
|
|
if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && (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 INPUT_ERROR;
|
|
}
|
|
}
|
|
|
|
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());
|
|
|
|
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();
|
|
|
|
x_restart.type(value());
|
|
after_restart = false;
|
|
|
|
reset_bias_force();
|
|
|
|
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
|
|
if (time_step_factor < 0) {
|
|
cvm::error("Error: timeStepFactor must be positive.\n");
|
|
return COLVARS_ERROR;
|
|
}
|
|
if (time_step_factor != 1) {
|
|
enable(f_cv_multiple_ts);
|
|
}
|
|
|
|
// TODO use here information from the CVCs' own natural boundaries
|
|
error_code |= init_grid_parameters(conf);
|
|
|
|
error_code |= init_extended_Lagrangian(conf);
|
|
error_code |= init_output_flags(conf);
|
|
|
|
// Now that the children are defined we can solve dependencies
|
|
enable(f_cv_active);
|
|
|
|
if (cvm::b_analysis)
|
|
parse_analysis(conf);
|
|
|
|
if (cvm::debug())
|
|
cvm::log("Done initializing collective variable \""+this->name+"\".\n");
|
|
|
|
return error_code;
|
|
}
|
|
|
|
|
|
#ifdef LEPTON
|
|
int colvar::init_custom_function(std::string const &conf)
|
|
{
|
|
std::string expr;
|
|
std::vector<Lepton::ParsedExpression> pexprs;
|
|
Lepton::ParsedExpression pexpr;
|
|
size_t pos = 0; // current position in config string
|
|
double *ref;
|
|
|
|
if (!key_lookup(conf, "customFunction", &expr, &pos)) {
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
enable(f_cv_custom_function);
|
|
cvm::log("This colvar uses a custom function.\n");
|
|
|
|
do {
|
|
if (cvm::debug())
|
|
cvm::log("Parsing expression \"" + expr + "\".\n");
|
|
try {
|
|
pexpr = Lepton::Parser::parse(expr);
|
|
pexprs.push_back(pexpr);
|
|
}
|
|
catch (...) {
|
|
cvm::error("Error parsing expression \"" + expr + "\".\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
try {
|
|
value_evaluators.push_back(
|
|
new Lepton::CompiledExpression(pexpr.createCompiledExpression()));
|
|
// Define variables for cvc values
|
|
// Stored in order: expr1, cvc1, cvc2, expr2, cvc1...
|
|
for (size_t i = 0; i < cvcs.size(); i++) {
|
|
for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
|
|
std::string vn = cvcs[i]->name +
|
|
(cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
|
|
try {
|
|
ref =&value_evaluators.back()->getVariableReference(vn);
|
|
}
|
|
catch (...) { // Variable is absent from expression
|
|
// To keep the same workflow, we use a pointer to a double here
|
|
// that will receive CVC values - even though none was allocated by Lepton
|
|
ref = &dev_null;
|
|
if (cvm::debug())
|
|
cvm::log("Variable " + vn + " is absent from expression \"" + expr + "\".\n");
|
|
}
|
|
value_eval_var_refs.push_back(ref);
|
|
}
|
|
}
|
|
}
|
|
catch (...) {
|
|
cvm::error("Error compiling expression \"" + expr + "\".\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
} while (key_lookup(conf, "customFunction", &expr, &pos));
|
|
|
|
|
|
// Now define derivative with respect to each scalar sub-component
|
|
for (size_t i = 0; i < cvcs.size(); i++) {
|
|
for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
|
|
std::string vn = cvcs[i]->name +
|
|
(cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
|
|
// Element ordering: we want the
|
|
// gradient vector of derivatives of all elements of the colvar
|
|
// wrt to a given element of a cvc ([i][j])
|
|
for (size_t c = 0; c < pexprs.size(); c++) {
|
|
gradient_evaluators.push_back(
|
|
new Lepton::CompiledExpression(pexprs[c].differentiate(vn).createCompiledExpression()));
|
|
// and record the refs to each variable in those expressions
|
|
for (size_t k = 0; k < cvcs.size(); k++) {
|
|
for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
|
|
std::string vvn = cvcs[k]->name +
|
|
(cvcs[k]->value().size() > 1 ? cvm::to_str(l+1) : "");
|
|
try {
|
|
ref = &gradient_evaluators.back()->getVariableReference(vvn);
|
|
}
|
|
catch (...) { // Variable is absent from derivative
|
|
// To keep the same workflow, we use a pointer to a double here
|
|
// that will receive CVC values - even though none was allocated by Lepton
|
|
if (cvm::debug())
|
|
cvm::log("Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n");
|
|
ref = &dev_null;
|
|
}
|
|
grad_eval_var_refs.push_back(ref);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (value_evaluators.size() == 0) {
|
|
cvm::error("Error: no custom function defined.\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
std::string type_str;
|
|
bool b_type_specified = get_keyval(conf, "customFunctionType",
|
|
type_str, "scalar", parse_silent);
|
|
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 custom colvar type.", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
// Guess type based on number of expressions
|
|
if (!b_type_specified) {
|
|
if (value_evaluators.size() == 1) {
|
|
x.type(colvarvalue::type_scalar);
|
|
} else {
|
|
x.type(colvarvalue::type_vector);
|
|
}
|
|
}
|
|
|
|
if (x.type() == colvarvalue::type_vector) {
|
|
x.vector1d_value.resize(value_evaluators.size());
|
|
}
|
|
|
|
x_reported.type(x);
|
|
cvm::log(std::string("Expecting colvar value of type ")
|
|
+ colvarvalue::type_desc(x.type())
|
|
+ (x.type()==colvarvalue::type_vector ? " of size " + cvm::to_str(x.size()) : "")
|
|
+ ".\n");
|
|
|
|
if (x.size() != value_evaluators.size()) {
|
|
cvm::error("Error: based on custom function type, expected "
|
|
+ cvm::to_str(x.size()) + " scalar expressions, but "
|
|
+ cvm::to_str(value_evaluators.size() + " were found.\n"));
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
#else
|
|
|
|
int colvar::init_custom_function(std::string const &conf)
|
|
{
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
#endif // #ifdef LEPTON
|
|
|
|
|
|
int colvar::init_grid_parameters(std::string const &conf)
|
|
{
|
|
colvarmodule *cv = cvm::main();
|
|
|
|
get_keyval(conf, "width", width, 1.0);
|
|
if (width <= 0.0) {
|
|
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
lower_boundary.type(value());
|
|
|
|
upper_boundary.type(value());
|
|
upper_wall.type(value());
|
|
|
|
set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
|
|
|
|
if (is_enabled(f_cv_scalar)) {
|
|
|
|
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
|
|
enable(f_cv_lower_boundary);
|
|
}
|
|
std::string lw_conf, uw_conf;
|
|
|
|
if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) {
|
|
cvm::log("Warning: lowerWallConstant and lowerWall are deprecated, "
|
|
"please define a harmonicWalls bias instead.\n");
|
|
lower_wall.type(value());
|
|
get_keyval(conf, "lowerWall", lower_wall, lower_boundary);
|
|
lw_conf = std::string("\n\
|
|
lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
|
|
lowerWalls "+cvm::to_str(lower_wall)+"\n");
|
|
}
|
|
|
|
if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
|
|
enable(f_cv_upper_boundary);
|
|
}
|
|
|
|
if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) {
|
|
cvm::log("Warning: upperWallConstant and upperWall are deprecated, "
|
|
"please define a harmonicWalls bias instead.\n");
|
|
upper_wall.type(value());
|
|
get_keyval(conf, "upperWall", upper_wall, upper_boundary);
|
|
uw_conf = std::string("\n\
|
|
upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
|
|
upperWalls "+cvm::to_str(upper_wall)+"\n");
|
|
}
|
|
|
|
if (lw_conf.size() && uw_conf.size()) {
|
|
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);
|
|
return INPUT_ERROR;
|
|
}
|
|
}
|
|
|
|
if (lw_conf.size() || uw_conf.size()) {
|
|
cvm::log("Generating a new harmonicWalls bias for compatibility purposes.\n");
|
|
std::string const walls_conf("\n\
|
|
harmonicWalls {\n\
|
|
name "+this->name+"w\n\
|
|
colvars "+this->name+"\n"+lw_conf+uw_conf+"\
|
|
timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+
|
|
"}\n");
|
|
cv->append_new_config(walls_conf);
|
|
}
|
|
}
|
|
|
|
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);
|
|
return 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);
|
|
return 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);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
|
|
int colvar::init_extended_Lagrangian(std::string const &conf)
|
|
{
|
|
get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false);
|
|
|
|
if (is_enabled(f_cv_extended_Lagrangian)) {
|
|
cvm::real temp, tolerance, period;
|
|
|
|
cvm::log("Enabling the extended Lagrangian term for colvar \""+
|
|
this->name+"\".\n");
|
|
|
|
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);
|
|
return INPUT_ERROR;
|
|
}
|
|
|
|
get_keyval(conf, "extendedFluctuation", tolerance);
|
|
if (tolerance <= 0.0) {
|
|
cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
|
|
cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/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) + " [E]/(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);
|
|
return INPUT_ERROR;
|
|
}
|
|
if (ext_gamma != 0.0) {
|
|
enable(f_cv_Langevin);
|
|
ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs
|
|
// Adjust Langevin sigma for slow time step if time_step_factor != 1
|
|
ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
|
|
}
|
|
}
|
|
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
|
|
int colvar::init_output_flags(std::string const &conf)
|
|
{
|
|
{
|
|
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 temp;
|
|
if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
|
|
cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
|
|
"The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
}
|
|
|
|
get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
|
|
get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
|
|
get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
|
|
|
|
return COLVARS_OK;
|
|
}
|
|
|
|
|
|
|
|
|
|
// 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("polar_phi")) &&
|
|
(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<polar_theta>(conf, "spherical polar angle theta",
|
|
"polarTheta");
|
|
error_code |= init_components_type<polar_phi>(conf, "spherical azimuthal angle phi",
|
|
"polarPhi");
|
|
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;
|
|
}
|
|
|
|
|
|
void colvar::do_feature_side_effects(int id)
|
|
{
|
|
switch (id) {
|
|
case f_cv_total_force_calc:
|
|
cvm::request_total_force();
|
|
break;
|
|
case f_cv_collect_gradient:
|
|
if (atom_ids.size() == 0) {
|
|
build_atom_list();
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
|
|
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);
|
|
}
|
|
|
|
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::proxy->backup_file(runave_outfile);
|
|
runave_os = cvm::proxy->output_stream(runave_outfile);
|
|
*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_bits(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()
|
|
{
|
|
// There is no need to call free_children_deps() here
|
|
// because the children are cvcs and will be deleted
|
|
// just below
|
|
|
|
// 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
|
|
colvarmodule *cv = cvm::main();
|
|
for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
|
|
cvi != cv->variables()->end();
|
|
++cvi) {
|
|
if ( *cvi == this) {
|
|
cv->variables()->erase(cvi);
|
|
break;
|
|
}
|
|
}
|
|
|
|
#ifdef LEPTON
|
|
for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
|
|
cei != value_evaluators.end();
|
|
++cei) {
|
|
if (*cei != NULL) delete (*cei);
|
|
}
|
|
value_evaluators.clear();
|
|
|
|
for (std::vector<Lepton::CompiledExpression *>::iterator gei = gradient_evaluators.begin();
|
|
gei != gradient_evaluators.end();
|
|
++gei) {
|
|
if (*gei != NULL) delete (*gei);
|
|
}
|
|
gradient_evaluators.clear();
|
|
#endif
|
|
}
|
|
|
|
|
|
|
|
// ******************** 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)) {
|
|
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;
|
|
}
|
|
|
|
|
|
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
|
|
{
|
|
colvarproxy *proxy = cvm::main()->proxy;
|
|
|
|
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");
|
|
|
|
error_code |= check_cvc_range(first_cvc, num_cvcs);
|
|
if (error_code != COLVARS_OK) {
|
|
return error_code;
|
|
}
|
|
|
|
if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
|
|
// Use Jacobian derivative from previous timestep
|
|
error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
|
|
}
|
|
// atom coordinates are updated by the next line
|
|
error_code |= calc_cvc_values(first_cvc, num_cvcs);
|
|
error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
|
|
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
|
|
if (proxy->total_forces_same_step()){
|
|
// Use Jacobian derivative from this timestep
|
|
error_code |= calc_cvc_total_force(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;
|
|
|
|
if (cvm::step_relative() > 0) {
|
|
// Total force depends on Jacobian derivative from previous timestep
|
|
// collect_cvc_total_forces() uses the previous value of jd
|
|
error_code |= collect_cvc_total_forces();
|
|
}
|
|
error_code |= collect_cvc_values();
|
|
error_code |= collect_cvc_gradients();
|
|
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");
|
|
|
|
// 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();
|
|
|
|
// 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;
|
|
}
|
|
|
|
#ifdef LEPTON
|
|
} else if (is_enabled(f_cv_custom_function)) {
|
|
|
|
size_t l = 0; // index in the vector of variable references
|
|
|
|
for (size_t i = 0; i < x.size(); i++) {
|
|
// Fill Lepton evaluator variables with CVC values, serialized into scalars
|
|
for (size_t j = 0; j < cvcs.size(); j++) {
|
|
for (size_t k = 0; k < cvcs[j]->value().size(); k++) {
|
|
*(value_eval_var_refs[l++]) = cvcs[j]->value()[k];
|
|
}
|
|
}
|
|
x[i] = value_evaluators[i]->evaluate();
|
|
}
|
|
#endif
|
|
|
|
} else if (x.type() == colvarvalue::type_scalar) {
|
|
// polynomial combination allowed
|
|
for (size_t i = 0; i < cvcs.size(); i++) {
|
|
if (!cvcs[i]->is_enabled()) continue;
|
|
x += (cvcs[i])->sup_coeff *
|
|
( ((cvcs[i])->sup_np != 1) ?
|
|
cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
|
|
(cvcs[i])->value().real_value );
|
|
}
|
|
} else {
|
|
for (size_t 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");
|
|
|
|
if (after_restart) {
|
|
if (cvm::proxy->simulation_running()) {
|
|
cvm::real const jump2 = dist2(x, x_restart) / (width*width);
|
|
if (jump2 > 0.25) {
|
|
cvm::error("Error: the calculated value of colvar \""+name+
|
|
"\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
|
|
"last read from the state file:\n"+cvm::to_str(x_restart)+
|
|
"\nPossible causes are changes in configuration, "
|
|
"wrong state file, or how PBC wrapping is handled.\n",
|
|
INPUT_ERROR);
|
|
return INPUT_ERROR;
|
|
}
|
|
}
|
|
}
|
|
|
|
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");
|
|
|
|
// 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++;
|
|
|
|
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
|
|
(cvcs[i])->calc_fit_gradients();
|
|
if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
|
|
(cvcs[i])->debug_gradients();
|
|
}
|
|
|
|
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)) {
|
|
// 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) *
|
|
cvm::integer_power((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_total_force(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_total_force_calc)) {
|
|
if (cvm::debug())
|
|
cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
|
|
|
|
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_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;
|
|
if (cvm::debug())
|
|
cvm::log("Colvar component no. "+cvm::to_str(i+1)+
|
|
" within colvar \""+this->name+"\" has total force "+
|
|
cvm::to_str((cvcs[i])->total_force(),
|
|
cvm::cv_width, cvm::cv_prec)+".\n");
|
|
// linear combination is assumed
|
|
ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
|
|
}
|
|
}
|
|
|
|
if (!is_enabled(f_cv_hide_Jacobian)) {
|
|
// 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;
|
|
}
|
|
}
|
|
|
|
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;
|
|
if (cvm::debug())
|
|
cvm::log("Colvar component no. "+cvm::to_str(i+1)+
|
|
" within colvar \""+this->name+"\" has Jacobian derivative"+
|
|
cvm::to_str((cvcs[i])->Jacobian_derivative(),
|
|
cvm::cv_width, cvm::cv_prec)+".\n");
|
|
// 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
|
|
if (cvm::step_relative() == 0 && !after_restart) {
|
|
xr = x;
|
|
vr.reset(); // (already 0; added for clarity)
|
|
}
|
|
|
|
// report the restraint center as "value"
|
|
x_reported = xr;
|
|
v_reported = vr;
|
|
// the "total force" with the extended Lagrangian is
|
|
// calculated in update_forces_energy() below
|
|
|
|
} 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;
|
|
}
|
|
|
|
// At the end of the first update after a restart, we can reset the flag
|
|
after_restart = false;
|
|
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();
|
|
fr.reset();
|
|
|
|
// If we are not active at this timestep, that's all we have to do
|
|
// return with energy == zero
|
|
if (!is_enabled(f_cv_active)) return 0.;
|
|
|
|
// 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 total force;
|
|
// instead, it is subtracted from the applied force (silent Jacobian correction)
|
|
// This requires the Jacobian term for the *current* timestep
|
|
if (is_enabled(f_cv_hide_Jacobian))
|
|
f -= fj;
|
|
}
|
|
|
|
// At this point f is the force f from external biases that will be applied to the
|
|
// extended variable if there is one
|
|
|
|
if (is_enabled(f_cv_extended_Lagrangian)) {
|
|
|
|
if (cvm::debug()) {
|
|
cvm::log("Updating extended-Lagrangian degree of freedom.\n");
|
|
}
|
|
|
|
if (prev_timestep > -1) {
|
|
// Keep track of slow timestep to integrate MTS colvars
|
|
// the colvar checks the interval after waking up twice
|
|
int n_timesteps = cvm::step_relative() - prev_timestep;
|
|
if (n_timesteps != 0 && n_timesteps != time_step_factor) {
|
|
cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
|
|
cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
|
|
" steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " +
|
|
cvm::to_str(cvm::step_relative()) + ").\n" +
|
|
"Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n");
|
|
return 0.;
|
|
}
|
|
}
|
|
prev_timestep = cvm::step_relative();
|
|
|
|
// Integrate with slow timestep (if time_step_factor != 1)
|
|
cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
|
|
|
|
colvarvalue f_ext(fr.type()); // force acting on the extended variable
|
|
f_ext.reset();
|
|
|
|
// the total force is applied to the fictitious mass, while the
|
|
// atoms only feel the harmonic force + wall 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 (fb_actual will be added just below)
|
|
fr = f;
|
|
// External force has been scaled for a 1-timestep impulse, scale it back because we will
|
|
// integrate it with the colvar's own timestep factor
|
|
f_ext = f / cvm::real(time_step_factor);
|
|
f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
|
|
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
|
|
// Coupling force is a slow force, to be applied to atomic coords impulse-style
|
|
f *= cvm::real(time_step_factor);
|
|
|
|
if (is_enabled(f_cv_subtract_applied_force)) {
|
|
// Report a "system" force without the biases on this colvar
|
|
// that is, just the spring force
|
|
ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
|
|
} else {
|
|
// 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;
|
|
colvarvalue rnd(x);
|
|
rnd.set_random();
|
|
vr += dt * ext_sigma * rnd / ext_mass;
|
|
}
|
|
vr += (0.5 * dt) * f_ext / ext_mass;
|
|
xr += dt * vr;
|
|
xr.apply_constraints();
|
|
if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
|
|
}
|
|
|
|
// Now adding the force on the actual colvar (for those biases that
|
|
// bypass the extended Lagrangian mass)
|
|
f += fb_actual;
|
|
|
|
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);
|
|
}
|
|
|
|
|
|
void colvar::communicate_forces()
|
|
{
|
|
size_t i;
|
|
if (cvm::debug()) {
|
|
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
|
|
cvm::log("Force to be applied: " + cvm::to_str(f) + "\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.as_vector() * func_grads[grad_index++],
|
|
cvcs[i]->value().type()));
|
|
}
|
|
|
|
#ifdef LEPTON
|
|
} else if (is_enabled(f_cv_custom_function)) {
|
|
|
|
size_t r = 0; // index in the vector of variable references
|
|
size_t e = 0; // index of the gradient evaluator
|
|
|
|
for (size_t i = 0; i < cvcs.size(); i++) { // gradient with respect to cvc i
|
|
cvm::matrix2d<cvm::real> jacobian (x.size(), cvcs[i]->value().size());
|
|
for (size_t j = 0; j < cvcs[i]->value().size(); j++) { // j-th element
|
|
for (size_t c = 0; c < x.size(); c++) { // derivative of scalar element c of the colvarvalue
|
|
|
|
// Feed cvc values to the evaluator
|
|
for (size_t k = 0; k < cvcs.size(); k++) { //
|
|
for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
|
|
*(grad_eval_var_refs[r++]) = cvcs[k]->value()[l];
|
|
}
|
|
}
|
|
jacobian[c][j] = gradient_evaluators[e++]->evaluate();
|
|
}
|
|
}
|
|
// cvc force is colvar force times colvar/cvc Jacobian
|
|
// (vector-matrix product)
|
|
(cvcs[i])->apply_force(colvarvalue(f.as_vector() * jacobian,
|
|
cvcs[i]->value().type()));
|
|
}
|
|
#endif
|
|
|
|
} 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 * (cvcs[i])->sup_coeff *
|
|
cvm::real((cvcs[i])->sup_np) *
|
|
(cvm::integer_power((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 * (cvcs[i])->sup_coeff);
|
|
}
|
|
}
|
|
|
|
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]->set_enabled(f_cvc_active, 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.clear();
|
|
}
|
|
|
|
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_bits(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, x, 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");
|
|
x_restart = x;
|
|
after_restart = true;
|
|
}
|
|
|
|
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_total_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_total_force)) {
|
|
os << " ft_"
|
|
<< 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_total_force)) {
|
|
os << " "
|
|
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
|
|
<< ft_reported;
|
|
}
|
|
|
|
if (is_enabled(f_cv_output_applied_force)) {
|
|
os << " "
|
|
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
|
|
<< applied_force();
|
|
}
|
|
|
|
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());
|
|
std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
|
|
if (!acf_os) return cvm::get_error();
|
|
write_acf(*acf_os);
|
|
cvm::proxy->close_output_stream(acf_outfile);
|
|
}
|
|
|
|
if (runave_os) {
|
|
cvm::proxy->close_output_stream(runave_outfile);
|
|
runave_os = NULL;
|
|
}
|
|
}
|
|
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<colvardeps::feature *> colvar::cv_features;
|