colvar.C

00001 // -*- c++ -*-
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <list>
00011 #include <vector>
00012 #include <algorithm>
00013 #include <iostream>
00014 #include <iomanip>
00015 
00016 #include "colvarmodule.h"
00017 #include "colvarvalue.h"
00018 #include "colvarparse.h"
00019 #include "colvar.h"
00020 #include "colvarcomp.h"
00021 #include "colvarscript.h"
00022 
00023 #if (__cplusplus >= 201103L)
00024 std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>> colvar::global_cvc_map = std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>>();
00025 #endif
00026 
00027 colvar::colvar()
00028 {
00029 prev_timestep = -1L;
00030 after_restart = false;
00031 kinetic_energy = 0.0;
00032 potential_energy = 0.0;
00033 
00034 #ifdef LEPTON
00035 dev_null = 0.0;
00036 #endif
00037 
00038 expand_boundaries = false;
00039 
00040 description = "uninitialized colvar";
00041 colvar::init_dependencies();
00042 }
00043 
00044 
00047 bool colvar::compare_cvc(const colvar::cvc* const i, const colvar::cvc* const j)
00048 {
00049 return i->name < j->name;
00050 }
00051 
00052 
00053 int colvar::init(std::string const &conf)
00054 {
00055 cvm::log("Initializing a new collective variable.\n");
00056 colvarparse::set_string(conf);
00057 
00058 int error_code = COLVARS_OK;
00059 
00060 colvarmodule *cv = cvm::main();
00061 
00062 get_keyval(conf, "name", this->name,
00063 (std::string("colvar")+cvm::to_str(cv->variables()->size())));
00064 
00065 if ((cvm::colvar_by_name(this->name) != NULL) &&
00066 (cvm::colvar_by_name(this->name) != this)) {
00067 cvm::error("Error: this colvar cannot have the same name, \""+this->name+
00068 "\", as another colvar.\n",
00069 COLVARS_INPUT_ERROR);
00070 return COLVARS_INPUT_ERROR;
00071 }
00072 
00073 // Initialize dependency members
00074 // Could be a function defined in a different source file, for space?
00075 
00076 this->description = "colvar " + this->name;
00077 
00078 error_code |= init_components(conf);
00079 if (error_code != COLVARS_OK) {
00080 return cvm::get_error();
00081 }
00082 
00083 size_t i;
00084 
00085 #ifdef LEPTON
00086 error_code |= init_custom_function(conf);
00087 if (error_code != COLVARS_OK) {
00088 return cvm::get_error();
00089 }
00090 #endif
00091 
00092 // Setup colvar as scripted function of components
00093 if (get_keyval(conf, "scriptedFunction", scripted_function,
00094 "", colvarparse::parse_silent)) {
00095 
00096 enable(f_cv_scripted);
00097 cvm::log("This colvar uses scripted function \"" + scripted_function + "\".\n");
00098 cvm::main()->cite_feature("Scripted functions (Tcl)");
00099 
00100 std::string type_str;
00101 get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
00102 
00103 x.type(colvarvalue::type_notset);
00104 int t;
00105 for (t = 0; t < colvarvalue::type_all; t++) {
00106 if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
00107 x.type(colvarvalue::Type(t));
00108 break;
00109 }
00110 }
00111 if (x.type() == colvarvalue::type_notset) {
00112 cvm::error("Could not parse scripted colvar type.", COLVARS_INPUT_ERROR);
00113 return COLVARS_INPUT_ERROR;
00114 }
00115 
00116 cvm::log(std::string("Expecting colvar value of type ")
00117 + colvarvalue::type_desc(x.type()));
00118 
00119 if (x.type() == colvarvalue::type_vector) {
00120 int size;
00121 if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
00122 cvm::error("Error: no size specified for vector scripted function.",
00123 COLVARS_INPUT_ERROR);
00124 return COLVARS_INPUT_ERROR;
00125 }
00126 x.vector1d_value.resize(size);
00127 }
00128 
00129 x_reported.type(x);
00130 
00131 // Sort array of cvcs based on their names
00132 // Note: default CVC names are in input order for same type of CVC
00133 std::sort(cvcs.begin(), cvcs.end(), colvar::compare_cvc);
00134 
00135 if(cvcs.size() > 1) {
00136 cvm::log("Sorted list of components for this scripted colvar:\n");
00137 for (i = 0; i < cvcs.size(); i++) {
00138 cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name);
00139 }
00140 }
00141 
00142 // Build ordered list of component values that will be
00143 // passed to the script
00144 for (i = 0; i < cvcs.size(); i++) {
00145 sorted_cvc_values.push_back(&(cvcs[i]->value()));
00146 }
00147 }
00148 
00149 if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function))) {
00150 colvarvalue const &cvc_value = (cvcs[0])->value();
00151 if (cvm::debug())
00152 cvm::log ("This collective variable is a "+
00153 colvarvalue::type_desc(cvc_value.type())+
00154 ((cvc_value.size() > 1) ? " with "+
00155 cvm::to_str(cvc_value.size())+" individual components.\n" :
00156 ".\n"));
00157 x.type(cvc_value);
00158 x_reported.type(cvc_value);
00159 }
00160 
00161 set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
00162 
00163 // If using scripted biases, any colvar may receive bias forces
00164 // and will need its gradient
00165 if (cvm::scripted_forces()) {
00166 enable(f_cv_gradient);
00167 }
00168 
00169 // check for linear combinations
00170 {
00171 bool lin = !(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function));
00172 for (i = 0; i < cvcs.size(); i++) {
00173 
00174 // FIXME this is a reverse dependency, ie. cv feature depends on cvc flag
00175 // need to clarify this case
00176 // if ((cvcs[i])->b_debug_gradients)
00177 // enable(task_gradients);
00178 
00179 if ((cvcs[i])->sup_np != 1) {
00180 if (cvm::debug() && lin)
00181 cvm::log("Warning: You are using a non-linear polynomial "
00182 "combination to define this collective variable, "
00183 "some biasing methods may be unavailable.\n");
00184 lin = false;
00185 
00186 if ((cvcs[i])->sup_np < 0) {
00187 cvm::log("Warning: you chose a negative exponent in the combination; "
00188 "if you apply forces, the simulation may become unstable "
00189 "when the component \""+
00190 (cvcs[i])->function_type+"\" approaches zero.\n");
00191 }
00192 }
00193 }
00194 set_enabled(f_cv_linear, lin);
00195 }
00196 
00197 // Colvar is homogeneous if:
00198 // - it is linear (hence not scripted)
00199 // - all cvcs have coefficient 1 or -1
00200 // i.e. sum or difference of cvcs
00201 {
00202 bool homogeneous = is_enabled(f_cv_linear);
00203 for (i = 0; i < cvcs.size(); i++) {
00204 if (cvm::fabs(cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {
00205 homogeneous = false;
00206 }
00207 }
00208 set_enabled(f_cv_homogeneous, homogeneous);
00209 }
00210 
00211 // A single-component variable almost concides with its CVC object
00212 if ((cvcs.size() == 1) && is_enabled(f_cv_homogeneous)) {
00213 if ( !is_enabled(f_cv_scripted) && !is_enabled(f_cv_custom_function) &&
00214 (cvm::fabs(cvcs[0]->sup_coeff - 1.0) < 1.0e-10) &&
00215 (cvcs[0]->sup_np == 1) ) {
00216 enable(f_cv_single_cvc);
00217 }
00218 }
00219 
00220 // Colvar is deemed periodic if:
00221 // - it is homogeneous
00222 // - all cvcs are periodic
00223 // - all cvcs have the same period
00224 if (is_enabled(f_cv_homogeneous) && cvcs[0]->is_enabled(f_cvc_periodic)) {
00225 bool b_periodic = true;
00226 period = cvcs[0]->period;
00227 wrap_center = cvcs[0]->wrap_center;
00228 for (i = 1; i < cvcs.size(); i++) {
00229 if (!cvcs[i]->is_enabled(f_cvc_periodic) || cvcs[i]->period != period) {
00230 b_periodic = false;
00231 period = 0.0;
00232 cvm::log("Warning: although one component is periodic, this colvar will "
00233 "not be treated as periodic, either because the exponent is not "
00234 "1, or because components of different periodicity are defined. "
00235 "Make sure that you know what you are doing!");
00236 }
00237 }
00238 set_enabled(f_cv_periodic, b_periodic);
00239 }
00240 
00241 // Allow scripted/custom functions to be defined as periodic
00242 if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) {
00243 if (get_keyval(conf, "period", period, 0.)) {
00244 enable(f_cv_periodic);
00245 get_keyval(conf, "wrapAround", wrap_center, 0.);
00246 }
00247 }
00248 
00249 // check that cvcs are compatible
00250 
00251 for (i = 0; i < cvcs.size(); i++) {
00252 
00253 // components may have different types only for scripted functions
00254 if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && (colvarvalue::check_types(cvcs[i]->value(),
00255 cvcs[0]->value())) ) {
00256 cvm::error("ERROR: you are defining this collective variable "
00257 "by using components of different types. "
00258 "You must use the same type in order to "
00259 "sum them together.\n", COLVARS_INPUT_ERROR);
00260 return COLVARS_INPUT_ERROR;
00261 }
00262 }
00263 
00264 active_cvc_square_norm = 0.;
00265 for (i = 0; i < cvcs.size(); i++) {
00266 active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
00267 }
00268 
00269 // at this point, the colvar's type is defined
00270 f.type(value());
00271 
00272 x_old.type(value());
00273 v_fdiff.type(value());
00274 v_reported.type(value());
00275 fj.type(value());
00276 ft.type(value());
00277 ft_reported.type(value());
00278 f_old.type(value());
00279 f_old.reset();
00280 
00281 x_restart.type(value());
00282 
00283 reset_bias_force();
00284 
00285 get_keyval(conf, "timeStepFactor", time_step_factor, 1);
00286 if (time_step_factor < 0) {
00287 cvm::error("Error: timeStepFactor must be positive.\n");
00288 return COLVARS_ERROR;
00289 }
00290 if (time_step_factor != 1) {
00291 enable(f_cv_multiple_ts);
00292 }
00293 
00294 error_code |= init_grid_parameters(conf);
00295 
00296 // Detect if we have a single component that is an alchemical lambda
00297 if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alchLambda") {
00298 enable(f_cv_external);
00299 }
00300 
00301 error_code |= init_extended_Lagrangian(conf);
00302 error_code |= init_output_flags(conf);
00303 
00304 // Now that the children are defined we can solve dependencies
00305 enable(f_cv_active);
00306 
00307 error_code |= parse_analysis(conf);
00308 
00309 if (cvm::debug())
00310 cvm::log("Done initializing collective variable \""+this->name+"\".\n");
00311 
00312 return error_code;
00313 }
00314 
00315 
00316 #ifdef LEPTON
00317 int colvar::init_custom_function(std::string const &conf)
00318 {
00319 std::string expr, expr_in; // expr_in is a buffer to remember expr after unsuccessful parsing
00320 std::vector<Lepton::ParsedExpression> pexprs;
00321 Lepton::ParsedExpression pexpr;
00322 size_t pos = 0; // current position in config string
00323 double *ref;
00324 
00325 if (!key_lookup(conf, "customFunction", &expr_in, &pos)) {
00326 return COLVARS_OK;
00327 }
00328 
00329 cvm::main()->cite_feature("Custom functions (Lepton)");
00330 
00331 enable(f_cv_custom_function);
00332 cvm::log("This colvar uses a custom function.\n");
00333 
00334 do {
00335 expr = expr_in;
00336 if (cvm::debug())
00337 cvm::log("Parsing expression \"" + expr + "\".\n");
00338 try {
00339 pexpr = Lepton::Parser::parse(expr);
00340 pexprs.push_back(pexpr);
00341 }
00342 catch (...) {
00343 cvm::error("Error parsing expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR);
00344 return COLVARS_INPUT_ERROR;
00345 }
00346 
00347 try {
00348 value_evaluators.push_back(
00349 new Lepton::CompiledExpression(pexpr.createCompiledExpression()));
00350 // Define variables for cvc values
00351 // Stored in order: expr1, cvc1, cvc2, expr2, cvc1...
00352 for (size_t i = 0; i < cvcs.size(); i++) {
00353 for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
00354 std::string vn = cvcs[i]->name +
00355 (cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
00356 try {
00357 ref =&value_evaluators.back()->getVariableReference(vn);
00358 }
00359 catch (...) { // Variable is absent from expression
00360 // To keep the same workflow, we use a pointer to a double here
00361 // that will receive CVC values - even though none was allocated by Lepton
00362 ref = &dev_null;
00363 cvm::log("Warning: Variable " + vn + " is absent from expression \"" + expr + "\".\n");
00364 }
00365 value_eval_var_refs.push_back(ref);
00366 }
00367 }
00368 }
00369 catch (...) {
00370 cvm::error("Error compiling expression \"" + expr + "\".\n", COLVARS_INPUT_ERROR);
00371 return COLVARS_INPUT_ERROR;
00372 }
00373 } while (key_lookup(conf, "customFunction", &expr_in, &pos));
00374 
00375 
00376 // Now define derivative with respect to each scalar sub-component
00377 for (size_t i = 0; i < cvcs.size(); i++) {
00378 for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
00379 std::string vn = cvcs[i]->name +
00380 (cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
00381 // Element ordering: we want the
00382 // gradient vector of derivatives of all elements of the colvar
00383 // wrt to a given element of a cvc ([i][j])
00384 for (size_t c = 0; c < pexprs.size(); c++) {
00385 gradient_evaluators.push_back(
00386 new Lepton::CompiledExpression(pexprs[c].differentiate(vn).createCompiledExpression()));
00387 // and record the refs to each variable in those expressions
00388 for (size_t k = 0; k < cvcs.size(); k++) {
00389 for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
00390 std::string vvn = cvcs[k]->name +
00391 (cvcs[k]->value().size() > 1 ? cvm::to_str(l+1) : "");
00392 try {
00393 ref = &gradient_evaluators.back()->getVariableReference(vvn);
00394 }
00395 catch (...) { // Variable is absent from derivative
00396 // To keep the same workflow, we use a pointer to a double here
00397 // that will receive CVC values - even though none was allocated by Lepton
00398 if (cvm::debug()) {
00399 cvm::log("Warning: Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n");
00400 }
00401 ref = &dev_null;
00402 }
00403 grad_eval_var_refs.push_back(ref);
00404 }
00405 }
00406 }
00407 }
00408 }
00409 
00410 
00411 if (value_evaluators.size() == 0) {
00412 cvm::error("Error: no custom function defined.\n", COLVARS_INPUT_ERROR);
00413 return COLVARS_INPUT_ERROR;
00414 }
00415 
00416 std::string type_str;
00417 bool b_type_specified = get_keyval(conf, "customFunctionType",
00418 type_str, "scalar", parse_silent);
00419 x.type(colvarvalue::type_notset);
00420 int t;
00421 for (t = 0; t < colvarvalue::type_all; t++) {
00422 if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
00423 x.type(colvarvalue::Type(t));
00424 break;
00425 }
00426 }
00427 if (x.type() == colvarvalue::type_notset) {
00428 cvm::error("Could not parse custom colvar type.", COLVARS_INPUT_ERROR);
00429 return COLVARS_INPUT_ERROR;
00430 }
00431 
00432 // Guess type based on number of expressions
00433 if (!b_type_specified) {
00434 if (value_evaluators.size() == 1) {
00435 x.type(colvarvalue::type_scalar);
00436 } else {
00437 x.type(colvarvalue::type_vector);
00438 }
00439 }
00440 
00441 if (x.type() == colvarvalue::type_vector) {
00442 x.vector1d_value.resize(value_evaluators.size());
00443 }
00444 
00445 x_reported.type(x);
00446 cvm::log(std::string("Expecting colvar value of type ")
00447 + colvarvalue::type_desc(x.type())
00448 + (x.type()==colvarvalue::type_vector ? " of size " + cvm::to_str(x.size()) : "")
00449 + ".\n");
00450 
00451 if (x.size() != value_evaluators.size()) {
00452 cvm::error("Error: based on custom function type, expected "
00453 + cvm::to_str(x.size()) + " scalar expressions, but "
00454 + cvm::to_str(value_evaluators.size()) + " were found.\n");
00455 return COLVARS_INPUT_ERROR;
00456 }
00457 
00458 return COLVARS_OK;
00459 }
00460 
00461 #else
00462 
00463 int colvar::init_custom_function(std::string const &conf)
00464 {
00465 
00466 std::string expr;
00467 size_t pos = 0;
00468 if (key_lookup(conf, "customFunction", &expr, &pos)) {
00469 std::string msg("Error: customFunction requires the Lepton library.");
00470 #if (__cplusplus < 201103L)
00471 // NOTE: this is not ideal; testing for the Lepton library's version would
00472 // be more accurate, but also less portable
00473 msg +=
00474 std::string(" Note also that recent versions of Lepton require C++11: "
00475 "please see https://colvars.github.io/README-c++11.html.");
00476 #endif
00477 return cvm::error(msg, COLVARS_NOT_IMPLEMENTED);
00478 }
00479 
00480 return COLVARS_OK;
00481 }
00482 
00483 #endif // #ifdef LEPTON
00484 
00485 
00486 int colvar::init_grid_parameters(std::string const &conf)
00487 {
00488 int error_code = COLVARS_OK;
00489 
00490 colvarmodule *cv = cvm::main();
00491 
00492 cvm::real default_width = width;
00493 
00494 if (!key_already_set("width")) {
00495 // The first time, check if the CVC has a width to provide
00496 default_width = 1.0;
00497 if (is_enabled(f_cv_single_cvc) && cvcs[0]->is_enabled(f_cvc_width)) {
00498 cvm::real const cvc_width = cvcs[0]->get_param("width");
00499 default_width = cvc_width;
00500 }
00501 }
00502 
00503 get_keyval(conf, "width", width, default_width);
00504 
00505 if (width <= 0.0) {
00506 cvm::error("Error: \"width\" must be positive.\n", COLVARS_INPUT_ERROR);
00507 return COLVARS_INPUT_ERROR;
00508 }
00509 
00510 lower_boundary.type(value());
00511 upper_boundary.type(value());
00512 lower_boundary.real_value = 0.0;
00513 upper_boundary.real_value = width; // Default to 1-wide grids
00514 
00515 if (is_enabled(f_cv_scalar)) {
00516 
00517 if (is_enabled(f_cv_single_cvc)) {
00518 // Get the default boundaries from the component
00519 if (cvcs[0]->is_enabled(f_cvc_lower_boundary)) {
00520 enable(f_cv_lower_boundary);
00521 enable(f_cv_hard_lower_boundary);
00522 lower_boundary =
00523 *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("lowerBoundary")));
00524 }
00525 if (cvcs[0]->is_enabled(f_cvc_upper_boundary)) {
00526 enable(f_cv_upper_boundary);
00527 enable(f_cv_hard_upper_boundary);
00528 upper_boundary =
00529 *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("upperBoundary")));
00530 }
00531 }
00532 
00533 if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
00534 enable(f_cv_lower_boundary);
00535 // Because this is the user's choice, we cannot assume it is a true
00536 // physical boundary
00537 disable(f_cv_hard_lower_boundary);
00538 }
00539 
00540 if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
00541 enable(f_cv_upper_boundary);
00542 disable(f_cv_hard_upper_boundary);
00543 }
00544 
00545 // Parse legacy wall options and set up a harmonicWalls bias if needed
00546 cvm::real lower_wall_k = 0.0, upper_wall_k = 0.0;
00547 cvm::real lower_wall = 0.0, upper_wall = 0.0;
00548 std::string lw_conf, uw_conf;
00549 
00550 if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
00551 parse_silent)) {
00552 cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
00553 "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
00554 if (!get_keyval(conf, "lowerWall", lower_wall)) {
00555 error_code |= cvm::error("Error: the value of lowerWall must be set "
00556 "explicitly.\n", COLVARS_INPUT_ERROR);
00557 }
00558 lw_conf = std::string("\n\
00559 lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
00560 lowerWalls "+cvm::to_str(lower_wall)+"\n");
00561 }
00562 
00563 if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
00564 parse_silent)) {
00565 cvm::log("Reading legacy options upperWall and upperWallConstant: "
00566 "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
00567 if (!get_keyval(conf, "upperWall", upper_wall)) {
00568 error_code |= cvm::error("Error: the value of upperWall must be set "
00569 "explicitly.\n", COLVARS_INPUT_ERROR);
00570 }
00571 uw_conf = std::string("\n\
00572 upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
00573 upperWalls "+cvm::to_str(upper_wall)+"\n");
00574 }
00575 
00576 if (lw_conf.size() && uw_conf.size()) {
00577 if (lower_wall >= upper_wall) {
00578 error_code |= cvm::error("Error: the upper wall, "+
00579 cvm::to_str(upper_wall)+
00580 ", is not higher than the lower wall, "+
00581 cvm::to_str(lower_wall)+".\n",
00582 COLVARS_INPUT_ERROR);
00583 }
00584 }
00585 
00586 if (lw_conf.size() || uw_conf.size()) {
00587 cvm::log("Generating a new harmonicWalls bias for compatibility purposes.\n");
00588 std::string const walls_conf("\n\
00589 harmonicWalls {\n\
00590 name "+this->name+"w\n\
00591 colvars "+this->name+"\n"+lw_conf+uw_conf+"\
00592 timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+
00593 "}\n");
00594 error_code |= cv->append_new_config(walls_conf);
00595 }
00596 }
00597 
00598 get_keyval_feature(this, conf, "hardLowerBoundary", f_cv_hard_lower_boundary,
00599 is_enabled(f_cv_hard_lower_boundary));
00600 
00601 get_keyval_feature(this, conf, "hardUpperBoundary", f_cv_hard_upper_boundary,
00602 is_enabled(f_cv_hard_upper_boundary));
00603 
00604 // consistency checks for boundaries and walls
00605 if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) {
00606 if (lower_boundary >= upper_boundary) {
00607 error_code |= cvm::error("Error: the upper boundary, "+
00608 cvm::to_str(upper_boundary)+
00609 ", is not higher than the lower boundary, "+
00610 cvm::to_str(lower_boundary)+".\n",
00611 COLVARS_INPUT_ERROR);
00612 }
00613 }
00614 
00615 get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries);
00616 if (expand_boundaries && periodic_boundaries()) {
00617 error_code |= cvm::error("Error: trying to expand boundaries that already "
00618 "cover a whole period of a periodic colvar.\n",
00619 COLVARS_INPUT_ERROR);
00620 }
00621 
00622 if (expand_boundaries && is_enabled(f_cv_hard_lower_boundary) &&
00623 is_enabled(f_cv_hard_upper_boundary)) {
00624 error_code |= cvm::error("Error: inconsistent configuration "
00625 "(trying to expand boundaries, but both "
00626 "hardLowerBoundary and hardUpperBoundary "
00627 "are enabled).\n", COLVARS_INPUT_ERROR);
00628 }
00629 
00630 return error_code;
00631 }
00632 
00633 
00634 int colvar::init_extended_Lagrangian(std::string const &conf)
00635 {
00636 colvarproxy *proxy = cvm::main()->proxy;
00637 get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false);
00638 
00639 if (is_enabled(f_cv_extended_Lagrangian)) {
00640 cvm::real temp, tolerance, extended_period;
00641 
00642 cvm::log("Enabling the extended Lagrangian term for colvar \""+
00643 this->name+"\".\n");
00644 
00645 // Mark x_ext as uninitialized so we can initialize it to the colvar value when updating
00646 x_ext.type(colvarvalue::type_notset);
00647 v_ext.type(value());
00648 fr.type(value());
00649 const bool temp_provided = get_keyval(conf, "extendedTemp", temp,
00650 proxy->target_temperature());
00651 if (is_enabled(f_cv_external)) {
00652 // In the case of an "external" coordinate, there is no coupling potential:
00653 // only the fictitious mass is meaningful
00654 get_keyval(conf, "extendedMass", ext_mass);
00655 // Ensure that the computed restraint energy term is zero
00656 ext_force_k = 0.0;
00657 } else {
00658 // Standard case of coupling to a geometric colvar
00659 if (temp <= 0.0) { // Then a finite temperature is required
00660 if (temp_provided)
00661 cvm::error("Error: \"extendedTemp\" must be positive.\n", COLVARS_INPUT_ERROR);
00662 else
00663 cvm::error("Error: a positive temperature must be provided, either "
00664 "by enabling a thermostat, or through \"extendedTemp\".\n",
00665 COLVARS_INPUT_ERROR);
00666 return COLVARS_INPUT_ERROR;
00667 }
00668 get_keyval(conf, "extendedFluctuation", tolerance);
00669 if (tolerance <= 0.0) {
00670 cvm::error("Error: \"extendedFluctuation\" must be positive.\n", COLVARS_INPUT_ERROR);
00671 return COLVARS_INPUT_ERROR;
00672 }
00673 ext_force_k = proxy->boltzmann() * temp / (tolerance * tolerance);
00674 cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2\n");
00675 
00676 get_keyval(conf, "extendedTimeConstant", extended_period, 200.0);
00677 if (extended_period <= 0.0) {
00678 cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", COLVARS_INPUT_ERROR);
00679 }
00680 ext_mass = (proxy->boltzmann() * temp * extended_period * extended_period)
00681 / (4.0 * PI * PI * tolerance * tolerance);
00682 cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)\n");
00683 }
00684 {
00685 bool b_output_energy;
00686 get_keyval(conf, "outputEnergy", b_output_energy, false);
00687 if (b_output_energy) {
00688 enable(f_cv_output_energy);
00689 }
00690 }
00691 
00692 get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
00693 if (ext_gamma < 0.0) {
00694 cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", COLVARS_INPUT_ERROR);
00695 return COLVARS_INPUT_ERROR;
00696 }
00697 if (ext_gamma != 0.0) {
00698 enable(f_cv_Langevin);
00699 ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs
00700 // Adjust Langevin sigma for slow time step if time_step_factor != 1
00701 ext_sigma = cvm::sqrt(2.0 * proxy->boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
00702 }
00703 
00704 get_keyval_feature(this, conf, "reflectingLowerBoundary", f_cv_reflecting_lower_boundary, false);
00705 get_keyval_feature(this, conf, "reflectingUpperBoundary", f_cv_reflecting_upper_boundary, false);
00706 }
00707 
00708 return COLVARS_OK;
00709 }
00710 
00711 
00712 int colvar::init_output_flags(std::string const &conf)
00713 {
00714 {
00715 bool b_output_value;
00716 get_keyval(conf, "outputValue", b_output_value, true);
00717 if (b_output_value) {
00718 enable(f_cv_output_value);
00719 }
00720 }
00721 
00722 {
00723 bool b_output_velocity;
00724 get_keyval(conf, "outputVelocity", b_output_velocity, false);
00725 if (b_output_velocity) {
00726 enable(f_cv_output_velocity);
00727 }
00728 }
00729 
00730 {
00731 bool temp;
00732 if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
00733 cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
00734 "The two are NOT identical: see https://colvars.github.io/totalforce.html.\n", COLVARS_INPUT_ERROR);
00735 return COLVARS_INPUT_ERROR;
00736 }
00737 }
00738 
00739 get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
00740 get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
00741 get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
00742 
00743 return COLVARS_OK;
00744 }
00745 
00746 #if (__cplusplus >= 201103L)
00747 // C++11
00748 template<typename def_class_name> int colvar::init_components_type(std::string const &,
00749 char const * /* def_desc */,
00750 char const *def_config_key) {
00751 // global_cvc_map is only supported in the C++11 case
00752 global_cvc_map[def_config_key] = [](const std::string& cvc_conf){return new def_class_name(cvc_conf);};
00753 // TODO: maybe it is better to do more check to avoid duplication in the map?
00754 return COLVARS_OK;
00755 }
00756 
00757 int colvar::init_components_type_from_global_map(const std::string& conf,
00758 const char* def_config_key) {
00759 #else
00760 template<typename def_class_name> int colvar::init_components_type(std::string const & conf,
00761 char const * /* def_desc */,
00762 char const *def_config_key) {
00763 #endif
00764 size_t def_count = 0;
00765 std::string def_conf = "";
00766 size_t pos = 0;
00767 while ( this->key_lookup(conf,
00768 def_config_key,
00769 &def_conf,
00770 &pos) ) {
00771 if (!def_conf.size()) continue;
00772 cvm::log("Initializing "
00773 "a new \""+std::string(def_config_key)+"\" component"+
00774 (cvm::debug() ? ", with configuration:\n"+def_conf
00775 : ".\n"));
00776 cvm::increase_depth();
00777 // only the following line is different from init_components_type
00778 // in the non-C++11 case
00779 #if (__cplusplus >= 201103L)
00780 cvc *cvcp = global_cvc_map.at(def_config_key)(def_conf);
00781 #else
00782 cvc *cvcp = new def_class_name(def_conf);
00783 #endif
00784 if (cvcp != NULL) {
00785 cvcs.push_back(cvcp);
00786 cvcp->check_keywords(def_conf, def_config_key);
00787 cvcp->set_function_type(def_config_key);
00788 if (cvm::get_error()) {
00789 cvm::error("Error: in setting up component \""+
00790 std::string(def_config_key)+"\".\n", COLVARS_INPUT_ERROR);
00791 return COLVARS_INPUT_ERROR;
00792 }
00793 cvm::decrease_depth();
00794 } else {
00795 cvm::decrease_depth();
00796 cvm::error("Error: in allocating component \""+
00797 std::string(def_config_key)+"\".\n",
00798 COLVARS_MEMORY_ERROR);
00799 return COLVARS_MEMORY_ERROR;
00800 }
00801 
00802 if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) {
00803 if (! cvcp->is_enabled(f_cvc_periodic)) {
00804 cvm::error("Error: invalid use of period and/or "
00805 "wrapAround in a \""+
00806 std::string(def_config_key)+
00807 "\" component.\n"+
00808 "Period: "+cvm::to_str(cvcp->period) +
00809 " wrapAround: "+cvm::to_str(cvcp->wrap_center),
00810 COLVARS_INPUT_ERROR);
00811 return COLVARS_INPUT_ERROR;
00812 }
00813 }
00814 
00815 if ( ! cvcs.back()->name.size()) {
00816 std::ostringstream s;
00817 s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;
00818 cvcs.back()->name = s.str();
00819 /* pad cvc number for correct ordering when sorting by name */
00820 }
00821 
00822 cvcs.back()->setup();
00823 if (cvm::debug()) {
00824 cvm::log("Done initializing a \""+
00825 std::string(def_config_key)+
00826 "\" component"+
00827 (cvm::debug() ?
00828 ", named \""+cvcs.back()->name+"\""
00829 : "")+".\n");
00830 }
00831 def_conf = "";
00832 if (cvm::debug()) {
00833 cvm::log("Parsed "+cvm::to_str(cvcs.size())+
00834 " components at this time.\n");
00835 }
00836 }
00837 
00838 return COLVARS_OK;
00839 }
00840 
00841 int colvar::init_components(std::string const &conf)
00842 {
00843 int error_code = COLVARS_OK;
00844 size_t i = 0, j = 0;
00845 
00846 // in the non-C++11 case, the components are initialized directly by init_components_type;
00847 // in the C++11 case, the components are stored in the global_cvc_map at first
00848 // by init_components_type, and then the map is iterated to initialize all components.
00849 error_code |= init_components_type<distance>(conf, "distance", "distance");
00850 error_code |= init_components_type<distance_vec>(conf, "distance vector", "distanceVec");
00851 error_code |= init_components_type<cartesian>(conf, "Cartesian coordinates", "cartesian");
00852 error_code |= init_components_type<distance_dir>(conf, "distance vector "
00853 "direction", "distanceDir");
00854 error_code |= init_components_type<distance_z>(conf, "distance projection "
00855 "on an axis", "distanceZ");
00856 error_code |= init_components_type<distance_xy>(conf, "distance projection "
00857 "on a plane", "distanceXY");
00858 error_code |= init_components_type<polar_theta>(conf, "spherical polar angle theta",
00859 "polarTheta");
00860 error_code |= init_components_type<polar_phi>(conf, "spherical azimuthal angle phi",
00861 "polarPhi");
00862 error_code |= init_components_type<distance_inv>(conf, "average distance "
00863 "weighted by inverse power", "distanceInv");
00864 error_code |= init_components_type<distance_pairs>(conf, "N1xN2-long vector "
00865 "of pairwise distances", "distancePairs");
00866 error_code |= init_components_type<dipole_magnitude>(conf, "dipole magnitude",
00867 "dipoleMagnitude");
00868 error_code |= init_components_type<coordnum>(conf, "coordination "
00869 "number", "coordNum");
00870 error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
00871 "number", "selfCoordNum");
00872 error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
00873 "number", "groupCoord");
00874 error_code |= init_components_type<angle>(conf, "angle", "angle");
00875 error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
00876 error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
00877 error_code |= init_components_type<h_bond>(conf, "hydrogen bond", "hBond");
00878 error_code |= init_components_type<alpha_angles>(conf, "alpha helix", "alpha");
00879 error_code |= init_components_type<dihedPC>(conf, "dihedral "
00880 "principal component", "dihedralPC");
00881 error_code |= init_components_type<orientation>(conf, "orientation", "orientation");
00882 error_code |= init_components_type<orientation_angle>(conf, "orientation "
00883 "angle", "orientationAngle");
00884 error_code |= init_components_type<orientation_proj>(conf, "orientation "
00885 "projection", "orientationProj");
00886 error_code |= init_components_type<tilt>(conf, "tilt", "tilt");
00887 error_code |= init_components_type<spin_angle>(conf, "spin angle", "spinAngle");
00888 error_code |= init_components_type<rmsd>(conf, "RMSD", "rmsd");
00889 error_code |= init_components_type<gyration>(conf, "radius of "
00890 "gyration", "gyration");
00891 error_code |= init_components_type<inertia>(conf, "moment of "
00892 "inertia", "inertia");
00893 error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");
00894 error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
00895 error_code |= init_components_type<alch_lambda>(conf, "alchemical coupling parameter", "alchLambda");
00896 error_code |= init_components_type<alch_Flambda>(conf, "force on alchemical coupling parameter", "alchFLambda");
00897 error_code |= init_components_type<gspath>(conf, "geometrical path collective variables (s)", "gspath");
00898 error_code |= init_components_type<gzpath>(conf, "geometrical path collective variables (z)", "gzpath");
00899 error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "linearCombination");
00900 error_code |= init_components_type<gspathCV>(conf, "geometrical path collective variables (s) for other CVs", "gspathCV");
00901 error_code |= init_components_type<gzpathCV>(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV");
00902 error_code |= init_components_type<aspathCV>(conf, "arithmetic path collective variables (s) for other CVs", "aspathCV");
00903 error_code |= init_components_type<azpathCV>(conf, "arithmetic path collective variables (s) for other CVs", "azpathCV");
00904 error_code |= init_components_type<euler_phi>(conf, "euler phi angle of the optimal orientation", "eulerPhi");
00905 error_code |= init_components_type<euler_psi>(conf, "euler psi angle of the optimal orientation", "eulerPsi");
00906 error_code |= init_components_type<euler_theta>(conf, "euler theta angle of the optimal orientation", "eulerTheta");
00907 #ifdef LEPTON
00908 error_code |= init_components_type<customColvar>(conf, "CV with support of the lepton custom function", "customColvar");
00909 #endif
00910 error_code |= init_components_type<neuralNetwork>(conf, "neural network CV for other CVs", "NeuralNetwork");
00911 
00912 error_code |= init_components_type<map_total>(conf, "total value of atomic map", "mapTotal");
00913 #if (__cplusplus >= 201103L)
00914 // iterate over all available CVC in the map
00915 for (auto it = global_cvc_map.begin(); it != global_cvc_map.end(); ++it) {
00916 error_code |= init_components_type_from_global_map(conf, it->first.c_str());
00917 // TODO: is it better to check the error code here?
00918 if (error_code != COLVARS_OK) {
00919 cvm::log("Failed to initialize " + it->first + " with the following configuration:\n");
00920 cvm::log(conf);
00921 // TODO: should it stop here?
00922 }
00923 }
00924 #endif
00925 if (!cvcs.size() || (error_code != COLVARS_OK)) {
00926 cvm::error("Error: no valid components were provided "
00927 "for this collective variable.\n",
00928 COLVARS_INPUT_ERROR);
00929 return COLVARS_INPUT_ERROR;
00930 }
00931 
00932 // Check for uniqueness of CVC names (esp. if user-provided)
00933 for (i = 0; i < cvcs.size(); i++) {
00934 for (j = i+1; j < cvcs.size(); j++) {
00935 if (cvcs[i]->name == cvcs[j]->name) {
00936 cvm::error("Components " + cvm::to_str(i) + " and " + cvm::to_str(j) +\
00937 " cannot have the same name \"" + cvcs[i]->name+ "\".\n", COLVARS_INPUT_ERROR);
00938 return COLVARS_INPUT_ERROR;
00939 }
00940 }
00941 }
00942 
00943 n_active_cvcs = cvcs.size();
00944 
00945 // Store list of children cvcs for dependency checking purposes
00946 for (i = 0; i < cvcs.size(); i++) {
00947 add_child(cvcs[i]);
00948 }
00949 
00950 cvm::log("All components initialized.\n");
00951 
00952 return COLVARS_OK;
00953 }
00954 
00955 
00956 void colvar::do_feature_side_effects(int id)
00957 {
00958 switch (id) {
00959 case f_cv_total_force_calc:
00960 cvm::request_total_force();
00961 break;
00962 case f_cv_collect_atom_ids:
00963 // Needed for getting gradients vias collect_gradients
00964 // or via atomic forces e.g. in Colvars Dashboard in VMD
00965 if (atom_ids.size() == 0) {
00966 build_atom_list();
00967 }
00968 break;
00969 }
00970 }
00971 
00972 
00973 void colvar::build_atom_list(void)
00974 {
00975 // If atomic gradients are requested, build full list of atom ids from all cvcs
00976 std::list<int> temp_id_list;
00977 
00978 for (size_t i = 0; i < cvcs.size(); i++) {
00979 for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
00980 cvm::atom_group const &ag = *(cvcs[i]->atom_groups[j]);
00981 for (size_t k = 0; k < ag.size(); k++) {
00982 temp_id_list.push_back(ag[k].id);
00983 }
00984 if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) {
00985 cvm::atom_group const &fg = *(ag.fitting_group);
00986 for (size_t k = 0; k < fg.size(); k++) {
00987 temp_id_list.push_back(fg[k].id);
00988 }
00989 }
00990 }
00991 }
00992 
00993 temp_id_list.sort();
00994 temp_id_list.unique();
00995 
00996 std::list<int>::iterator li;
00997 for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
00998 atom_ids.push_back(*li);
00999 }
01000 
01001 temp_id_list.clear();
01002 
01003 atomic_gradients.resize(atom_ids.size());
01004 if (atom_ids.size()) {
01005 if (cvm::debug())
01006 cvm::log("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
01007 } else {
01008 cvm::log("Warning: colvar components communicated no atom IDs.\n");
01009 }
01010 }
01011 
01012 
01013 int colvar::parse_analysis(std::string const &conf)
01014 {
01015 
01016 // if (cvm::debug())
01017 // cvm::log ("Parsing analysis flags for collective variable \""+
01018 // this->name+"\".\n");
01019 
01020 runave_length = 0;
01021 bool b_runave = false;
01022 if (get_keyval(conf, "runAve", b_runave) && b_runave) {
01023 
01024 enable(f_cv_runave);
01025 
01026 get_keyval(conf, "runAveLength", runave_length, 1000);
01027 get_keyval(conf, "runAveStride", runave_stride, 1);
01028 
01029 if ((cvm::restart_out_freq % runave_stride) != 0) {
01030 cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", COLVARS_INPUT_ERROR);
01031 }
01032 
01033 get_keyval(conf, "runAveOutputFile", runave_outfile, runave_outfile);
01034 }
01035 
01036 acf_length = 0;
01037 bool b_acf = false;
01038 if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {
01039 
01040 enable(f_cv_corrfunc);
01041 
01042 get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
01043 if (acf_colvar_name == this->name) {
01044 cvm::log("Calculating auto-correlation function.\n");
01045 } else {
01046 cvm::log("Calculating correlation function with \""+
01047 this->name+"\".\n");
01048 }
01049 
01050 std::string acf_type_str;
01051 get_keyval(conf, "corrFuncType", acf_type_str, to_lower_cppstr(std::string("velocity")));
01052 if (acf_type_str == to_lower_cppstr(std::string("coordinate"))) {
01053 acf_type = acf_coor;
01054 } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
01055 acf_type = acf_vel;
01056 enable(f_cv_fdiff_velocity);
01057 colvar *cv2 = cvm::colvar_by_name(acf_colvar_name);
01058 if (cv2 == NULL) {
01059 return cvm::error("Error: collective variable \""+acf_colvar_name+
01060 "\" is not defined at this time.\n", COLVARS_INPUT_ERROR);
01061 }
01062 cv2->enable(f_cv_fdiff_velocity); // Manual dependency to object of same type
01063 } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
01064 acf_type = acf_p2coor;
01065 } else {
01066 cvm::log("Unknown type of correlation function, \""+
01067 acf_type_str+"\".\n");
01068 cvm::set_error_bits(COLVARS_INPUT_ERROR);
01069 }
01070 
01071 get_keyval(conf, "corrFuncOffset", acf_offset, 0);
01072 get_keyval(conf, "corrFuncLength", acf_length, 1000);
01073 get_keyval(conf, "corrFuncStride", acf_stride, 1);
01074 
01075 if ((cvm::restart_out_freq % acf_stride) != 0) {
01076 cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", COLVARS_INPUT_ERROR);
01077 }
01078 
01079 get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
01080 get_keyval(conf, "corrFuncOutputFile", acf_outfile, acf_outfile);
01081 }
01082 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
01083 }
01084 
01085 
01086 int colvar::init_dependencies() {
01087 size_t i;
01088 if (features().size() == 0) {
01089 for (i = 0; i < f_cv_ntot; i++) {
01090 modify_features().push_back(new feature);
01091 }
01092 
01093 init_feature(f_cv_active, "active", f_type_dynamic);
01094 // Do not require f_cvc_active in children, as some components may be disabled
01095 // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
01096 require_feature_alt(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
01097 
01098 init_feature(f_cv_awake, "awake", f_type_static);
01099 require_feature_self(f_cv_awake, f_cv_active);
01100 
01101 init_feature(f_cv_gradient, "gradient", f_type_dynamic);
01102 require_feature_children(f_cv_gradient, f_cvc_gradient);
01103 
01104 init_feature(f_cv_collect_gradient, "collect_gradient", f_type_dynamic);
01105 require_feature_self(f_cv_collect_gradient, f_cv_gradient);
01106 require_feature_self(f_cv_collect_gradient, f_cv_scalar);
01107 require_feature_self(f_cv_collect_gradient, f_cv_collect_atom_ids);
01108 // The following exclusions could be lifted by implementing the feature
01109 exclude_feature_self(f_cv_collect_gradient, f_cv_scripted);
01110 exclude_feature_self(f_cv_collect_gradient, f_cv_custom_function);
01111 require_feature_children(f_cv_collect_gradient, f_cvc_explicit_gradient);
01112 
01113 init_feature(f_cv_collect_atom_ids, "collect_atom_ids", f_type_dynamic);
01114 require_feature_children(f_cv_collect_atom_ids, f_cvc_collect_atom_ids);
01115 
01116 init_feature(f_cv_fdiff_velocity, "velocity_from_finite_differences", f_type_dynamic);
01117 
01118 // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
01119 init_feature(f_cv_total_force, "total_force", f_type_dynamic);
01120 require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
01121 
01122 // Deps for explicit total force calculation
01123 init_feature(f_cv_total_force_calc, "total_force_calculation", f_type_dynamic);
01124 require_feature_self(f_cv_total_force_calc, f_cv_scalar);
01125 require_feature_self(f_cv_total_force_calc, f_cv_linear);
01126 require_feature_children(f_cv_total_force_calc, f_cvc_inv_gradient);
01127 require_feature_self(f_cv_total_force_calc, f_cv_Jacobian);
01128 
01129 init_feature(f_cv_Jacobian, "Jacobian_derivative", f_type_dynamic);
01130 require_feature_self(f_cv_Jacobian, f_cv_scalar);
01131 require_feature_self(f_cv_Jacobian, f_cv_linear);
01132 require_feature_children(f_cv_Jacobian, f_cvc_Jacobian);
01133 
01134 init_feature(f_cv_hide_Jacobian, "hide_Jacobian_force", f_type_user);
01135 require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
01136 exclude_feature_self(f_cv_hide_Jacobian, f_cv_extended_Lagrangian);
01137 
01138 init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user);
01139 require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
01140 require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient);
01141 
01142 init_feature(f_cv_Langevin, "Langevin_dynamics", f_type_user);
01143 require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
01144 
01145 init_feature(f_cv_external, "external", f_type_user);
01146 require_feature_self(f_cv_external, f_cv_single_cvc);
01147 
01148 init_feature(f_cv_single_cvc, "single_component", f_type_static);
01149 
01150 init_feature(f_cv_linear, "linear", f_type_static);
01151 
01152 init_feature(f_cv_scalar, "scalar", f_type_static);
01153 
01154 init_feature(f_cv_output_energy, "output_energy", f_type_user);
01155 
01156 init_feature(f_cv_output_value, "output_value", f_type_user);
01157 
01158 init_feature(f_cv_output_velocity, "output_velocity", f_type_user);
01159 require_feature_self(f_cv_output_velocity, f_cv_fdiff_velocity);
01160 
01161 init_feature(f_cv_output_applied_force, "output_applied_force", f_type_user);
01162 
01163 init_feature(f_cv_output_total_force, "output_total_force", f_type_user);
01164 require_feature_self(f_cv_output_total_force, f_cv_total_force);
01165 
01166 init_feature(f_cv_subtract_applied_force, "subtract_applied_force_from_total_force", f_type_user);
01167 require_feature_self(f_cv_subtract_applied_force, f_cv_total_force);
01168 
01169 init_feature(f_cv_lower_boundary, "lower_boundary", f_type_user);
01170 require_feature_self(f_cv_lower_boundary, f_cv_scalar);
01171 
01172 init_feature(f_cv_upper_boundary, "upper_boundary", f_type_user);
01173 require_feature_self(f_cv_upper_boundary, f_cv_scalar);
01174 
01175 init_feature(f_cv_hard_lower_boundary, "hard_lower_boundary", f_type_user);
01176 require_feature_self(f_cv_hard_lower_boundary, f_cv_lower_boundary);
01177 
01178 init_feature(f_cv_hard_upper_boundary, "hard_upper_boundary", f_type_user);
01179 require_feature_self(f_cv_hard_upper_boundary, f_cv_upper_boundary);
01180 
01181 init_feature(f_cv_reflecting_lower_boundary, "reflecting_lower_boundary", f_type_user);
01182 require_feature_self(f_cv_reflecting_lower_boundary, f_cv_lower_boundary);
01183 require_feature_self(f_cv_reflecting_lower_boundary, f_cv_extended_Lagrangian);
01184 
01185 init_feature(f_cv_reflecting_upper_boundary, "reflecting_upper_boundary", f_type_user);
01186 require_feature_self(f_cv_reflecting_upper_boundary, f_cv_upper_boundary);
01187 require_feature_self(f_cv_reflecting_upper_boundary, f_cv_extended_Lagrangian);
01188 
01189 init_feature(f_cv_grid, "grid", f_type_dynamic);
01190 require_feature_self(f_cv_grid, f_cv_scalar);
01191 
01192 init_feature(f_cv_runave, "running_average", f_type_user);
01193 
01194 init_feature(f_cv_corrfunc, "correlation_function", f_type_user);
01195 
01196 init_feature(f_cv_scripted, "scripted", f_type_user);
01197 
01198 init_feature(f_cv_custom_function, "custom_function", f_type_user);
01199 exclude_feature_self(f_cv_custom_function, f_cv_scripted);
01200 
01201 init_feature(f_cv_periodic, "periodic", f_type_static);
01202 require_feature_self(f_cv_periodic, f_cv_scalar);
01203 init_feature(f_cv_scalar, "scalar", f_type_static);
01204 init_feature(f_cv_linear, "linear", f_type_static);
01205 init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
01206 
01207 // because total forces are obtained from the previous time step,
01208 // we cannot (currently) have colvar values and total forces for the same timestep
01209 init_feature(f_cv_multiple_ts, "multiple_timestep", f_type_static);
01210 exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
01211 
01212 // check that everything is initialized
01213 for (i = 0; i < colvardeps::f_cv_ntot; i++) {
01214 if (is_not_set(i)) {
01215 cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
01216 }
01217 }
01218 }
01219 
01220 // Initialize feature_states for each instance
01221 feature_states.reserve(f_cv_ntot);
01222 for (i = 0; i < f_cv_ntot; i++) {
01223 feature_states.push_back(feature_state(true, false));
01224 // Most features are available, so we set them so
01225 // and list exceptions below
01226 }
01227 
01228 feature_states[f_cv_fdiff_velocity].available =
01229 cvm::main()->proxy->simulation_running();
01230 
01231 return COLVARS_OK;
01232 }
01233 
01234 
01235 void colvar::setup()
01236 {
01237 // loop over all components to update masses and charges of all groups
01238 for (size_t i = 0; i < cvcs.size(); i++) {
01239 for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
01240 cvm::atom_group *atoms = cvcs[i]->atom_groups[ig];
01241 atoms->setup();
01242 atoms->print_properties(name, i, ig);
01243 atoms->read_positions();
01244 }
01245 }
01246 }
01247 
01248 
01249 std::vector<std::vector<int> > colvar::get_atom_lists()
01250 {
01251 std::vector<std::vector<int> > lists;
01252 for (size_t i = 0; i < cvcs.size(); i++) {
01253 std::vector<std::vector<int> > li = cvcs[i]->get_atom_lists();
01254 lists.insert(lists.end(), li.begin(), li.end());
01255 }
01256 return lists;
01257 }
01258 
01259 
01260 std::vector<int> const &colvar::get_volmap_ids()
01261 {
01262 volmap_ids_.resize(cvcs.size());
01263 for (size_t i = 0; i < cvcs.size(); i++) {
01264 if (cvcs[i]->param_exists("mapID") == COLVARS_OK) {
01265 volmap_ids_[i] =
01266 *(reinterpret_cast<int const *>(cvcs[i]->get_param_ptr("mapID")));
01267 } else {
01268 volmap_ids_[i] = -1;
01269 }
01270 }
01271 return volmap_ids_;
01272 }
01273 
01274 
01275 colvar::~colvar()
01276 {
01277 // There is no need to call free_children_deps() here
01278 // because the children are cvcs and will be deleted
01279 // just below
01280 
01281 // Clear references to this colvar's cvcs as children
01282 // for dependency purposes
01283 remove_all_children();
01284 
01285 for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();
01286 ci != cvcs.rend();
01287 ++ci) {
01288 // clear all children of this cvc (i.e. its atom groups)
01289 // because the cvc base class destructor can't do it early enough
01290 // and we don't want to have each cvc derived class do it separately
01291 (*ci)->remove_all_children();
01292 delete *ci;
01293 }
01294 cvcs.clear();
01295 
01296 while (biases.size() > 0) {
01297 size_t const i = biases.size()-1;
01298 cvm::log("Warning: before deleting colvar " + name
01299 + ", deleting related bias " + biases[i]->name);
01300 delete biases[i];
01301 }
01302 biases.clear();
01303 
01304 // remove reference to this colvar from the module
01305 colvarmodule *cv = cvm::main();
01306 for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
01307 cvi != cv->variables()->end();
01308 ++cvi) {
01309 if ( *cvi == this) {
01310 cv->variables()->erase(cvi);
01311 break;
01312 }
01313 }
01314 
01315 cv->config_changed();
01316 
01317 #ifdef LEPTON
01318 for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
01319 cei != value_evaluators.end();
01320 ++cei) {
01321 if (*cei != NULL) delete (*cei);
01322 }
01323 value_evaluators.clear();
01324 
01325 for (std::vector<Lepton::CompiledExpression *>::iterator gei = gradient_evaluators.begin();
01326 gei != gradient_evaluators.end();
01327 ++gei) {
01328 if (*gei != NULL) delete (*gei);
01329 }
01330 gradient_evaluators.clear();
01331 #endif
01332 }
01333 
01334 
01335 
01336 // ******************** CALC FUNCTIONS ********************
01337 
01338 
01339 // Default schedule (everything is serialized)
01340 int colvar::calc()
01341 {
01342 // Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
01343 int error_code = COLVARS_OK;
01344 if (is_enabled(f_cv_active)) {
01345 error_code |= update_cvc_flags();
01346 if (error_code != COLVARS_OK) return error_code;
01347 error_code |= calc_cvcs();
01348 if (error_code != COLVARS_OK) return error_code;
01349 error_code |= collect_cvc_data();
01350 }
01351 return error_code;
01352 }
01353 
01354 
01355 int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
01356 {
01357 if (cvm::debug())
01358 cvm::log("Calculating colvar \""+this->name+"\", components "+
01359 cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
01360 
01361 colvarproxy *proxy = cvm::main()->proxy;
01362 int error_code = COLVARS_OK;
01363 
01364 error_code |= check_cvc_range(first_cvc, num_cvcs);
01365 if (error_code != COLVARS_OK) {
01366 return error_code;
01367 }
01368 
01369 if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
01370 // Use Jacobian derivative from previous timestep
01371 error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
01372 }
01373 // atom coordinates are updated by the next line
01374 error_code |= calc_cvc_values(first_cvc, num_cvcs);
01375 error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
01376 error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
01377 if (proxy->total_forces_same_step()){
01378 // Use Jacobian derivative from this timestep
01379 error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
01380 }
01381 
01382 if (cvm::debug())
01383 cvm::log("Done calculating colvar \""+this->name+"\".\n");
01384 
01385 return error_code;
01386 }
01387 
01388 
01389 int colvar::collect_cvc_data()
01390 {
01391 if (cvm::debug())
01392 cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
01393 
01394 colvarproxy *proxy = cvm::main()->proxy;
01395 int error_code = COLVARS_OK;
01396 
01397 if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
01398 // Total force depends on Jacobian derivative from previous timestep
01399 // collect_cvc_total_forces() uses the previous value of jd
01400 error_code |= collect_cvc_total_forces();
01401 }
01402 error_code |= collect_cvc_values();
01403 error_code |= collect_cvc_gradients();
01404 error_code |= collect_cvc_Jacobians();
01405 if (proxy->total_forces_same_step()){
01406 // Use Jacobian derivative from this timestep
01407 error_code |= collect_cvc_total_forces();
01408 }
01409 error_code |= calc_colvar_properties();
01410 
01411 if (cvm::debug())
01412 cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
01413 
01414 return error_code;
01415 }
01416 
01417 
01418 int colvar::check_cvc_range(int first_cvc, size_t /* num_cvcs */)
01419 {
01420 if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) {
01421 cvm::error("Error: trying to address a component outside the "
01422 "range defined for colvar \""+name+"\".\n", COLVARS_BUG_ERROR);
01423 return COLVARS_BUG_ERROR;
01424 }
01425 return COLVARS_OK;
01426 }
01427 
01428 
01429 int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)
01430 {
01431 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
01432 size_t i, cvc_count;
01433 
01434 // calculate the value of the colvar
01435 
01436 if (cvm::debug())
01437 cvm::log("Calculating colvar components.\n");
01438 
01439 // First, calculate component values
01440 cvm::increase_depth();
01441 for (i = first_cvc, cvc_count = 0;
01442 (i < cvcs.size()) && (cvc_count < cvc_max_count);
01443 i++) {
01444 if (!cvcs[i]->is_enabled()) continue;
01445 cvc_count++;
01446 (cvcs[i])->read_data();
01447 (cvcs[i])->calc_value();
01448 if (cvm::debug())
01449 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
01450 " within colvar \""+this->name+"\" has value "+
01451 cvm::to_str((cvcs[i])->value(),
01452 cvm::cv_width, cvm::cv_prec)+".\n");
01453 }
01454 cvm::decrease_depth();
01455 
01456 return COLVARS_OK;
01457 }
01458 
01459 
01460 int colvar::collect_cvc_values()
01461 {
01462 x.reset();
01463 
01464 // combine them appropriately, using either a scripted function or a polynomial
01465 if (is_enabled(f_cv_scripted)) {
01466 // cvcs combined by user script
01467 int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);
01468 if (res == COLVARS_NOT_IMPLEMENTED) {
01469 cvm::error("Scripted colvars are not implemented.");
01470 return COLVARS_NOT_IMPLEMENTED;
01471 }
01472 if (res != COLVARS_OK) {
01473 cvm::error("Error running scripted colvar");
01474 return COLVARS_OK;
01475 }
01476 
01477 #ifdef LEPTON
01478 } else if (is_enabled(f_cv_custom_function)) {
01479 
01480 size_t l = 0; // index in the vector of variable references
01481 
01482 for (size_t i = 0; i < x.size(); i++) {
01483 // Fill Lepton evaluator variables with CVC values, serialized into scalars
01484 for (size_t j = 0; j < cvcs.size(); j++) {
01485 for (size_t k = 0; k < cvcs[j]->value().size(); k++) {
01486 *(value_eval_var_refs[l++]) = cvcs[j]->value()[k];
01487 }
01488 }
01489 x[i] = value_evaluators[i]->evaluate();
01490 }
01491 #endif
01492 
01493 } else if (x.type() == colvarvalue::type_scalar) {
01494 // polynomial combination allowed
01495 for (size_t i = 0; i < cvcs.size(); i++) {
01496 if (!cvcs[i]->is_enabled()) continue;
01497 x += (cvcs[i])->sup_coeff *
01498 ( ((cvcs[i])->sup_np != 1) ?
01499 cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
01500 (cvcs[i])->value().real_value );
01501 }
01502 } else {
01503 for (size_t i = 0; i < cvcs.size(); i++) {
01504 if (!cvcs[i]->is_enabled()) continue;
01505 x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
01506 }
01507 }
01508 
01509 if (cvm::debug())
01510 cvm::log("Colvar \""+this->name+"\" has value "+
01511 cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
01512 
01513 if (after_restart) {
01514 if (cvm::proxy->simulation_running()) {
01515 cvm::real const jump2 = dist2(x, x_restart) / (width*width);
01516 if (jump2 > 0.25) {
01517 cvm::error("Error: the calculated value of colvar \""+name+
01518 "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
01519 "last read from the state file:\n"+cvm::to_str(x_restart)+
01520 "\nPossible causes are changes in configuration, "
01521 "wrong state file, or how PBC wrapping is handled.\n",
01522 COLVARS_INPUT_ERROR);
01523 return COLVARS_INPUT_ERROR;
01524 }
01525 }
01526 }
01527 
01528 return COLVARS_OK;
01529 }
01530 
01531 
01532 int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
01533 {
01534 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
01535 size_t i, cvc_count;
01536 
01537 if (cvm::debug())
01538 cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
01539 
01540 // calculate the gradients of each component
01541 cvm::increase_depth();
01542 for (i = first_cvc, cvc_count = 0;
01543 (i < cvcs.size()) && (cvc_count < cvc_max_count);
01544 i++) {
01545 if (!cvcs[i]->is_enabled()) continue;
01546 cvc_count++;
01547 
01548 if ((cvcs[i])->is_enabled(f_cvc_gradient)) {
01549 (cvcs[i])->calc_gradients();
01550 // if requested, propagate (via chain rule) the gradients above
01551 // to the atoms used to define the roto-translation
01552 (cvcs[i])->calc_fit_gradients();
01553 if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
01554 (cvcs[i])->debug_gradients();
01555 }
01556 
01557 cvm::decrease_depth();
01558 
01559 if (cvm::debug())
01560 cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
01561 }
01562 
01563 return COLVARS_OK;
01564 }
01565 
01566 
01567 int colvar::collect_cvc_gradients()
01568 {
01569 size_t i;
01570 if (is_enabled(f_cv_collect_gradient)) {
01571 // Collect the atomic gradients inside colvar object
01572 for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
01573 atomic_gradients[a].reset();
01574 }
01575 for (i = 0; i < cvcs.size(); i++) {
01576 if (!cvcs[i]->is_enabled()) continue;
01577 cvcs[i]->collect_gradients(atom_ids, atomic_gradients);
01578 }
01579 }
01580 return COLVARS_OK;
01581 }
01582 
01583 
01584 int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
01585 {
01586 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
01587 size_t i, cvc_count;
01588 
01589 if (is_enabled(f_cv_total_force_calc)) {
01590 if (cvm::debug())
01591 cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
01592 
01593 cvm::increase_depth();
01594 
01595 for (i = first_cvc, cvc_count = 0;
01596 (i < cvcs.size()) && (cvc_count < cvc_max_count);
01597 i++) {
01598 if (!cvcs[i]->is_enabled()) continue;
01599 cvc_count++;
01600 (cvcs[i])->calc_force_invgrads();
01601 }
01602 cvm::decrease_depth();
01603 
01604 
01605 if (cvm::debug())
01606 cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
01607 }
01608 
01609 return COLVARS_OK;
01610 }
01611 
01612 
01613 int colvar::collect_cvc_total_forces()
01614 {
01615 if (is_enabled(f_cv_total_force_calc)) {
01616 ft.reset();
01617 
01618 if (cvm::step_relative() > 0) {
01619 // get from the cvcs the total forces from the PREVIOUS step
01620 for (size_t i = 0; i < cvcs.size(); i++) {
01621 if (!cvcs[i]->is_enabled()) continue;
01622 if (cvm::debug())
01623 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
01624 " within colvar \""+this->name+"\" has total force "+
01625 cvm::to_str((cvcs[i])->total_force(),
01626 cvm::cv_width, cvm::cv_prec)+".\n");
01627 // linear combination is assumed
01628 ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
01629 }
01630 }
01631 
01632 if (!(is_enabled(f_cv_hide_Jacobian) && is_enabled(f_cv_subtract_applied_force))) {
01633 // add the Jacobian force to the total force, and don't apply any silent
01634 // correction internally: biases such as colvarbias_abf will handle it
01635 // If f_cv_hide_Jacobian is enabled, a force of -fj is present in ft due to the
01636 // Jacobian-compensating force
01637 ft += fj;
01638 }
01639 }
01640 
01641 return COLVARS_OK;
01642 }
01643 
01644 
01645 int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs)
01646 {
01647 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
01648 
01649 if (is_enabled(f_cv_Jacobian)) {
01650 cvm::increase_depth();
01651 size_t i, cvc_count;
01652 for (i = first_cvc, cvc_count = 0;
01653 (i < cvcs.size()) && (cvc_count < cvc_max_count);
01654 i++) {
01655 if (!cvcs[i]->is_enabled()) continue;
01656 cvc_count++;
01657 (cvcs[i])->calc_Jacobian_derivative();
01658 }
01659 cvm::decrease_depth();
01660 }
01661 
01662 return COLVARS_OK;
01663 }
01664 
01665 
01666 int colvar::collect_cvc_Jacobians()
01667 {
01668 colvarproxy *proxy = cvm::main()->proxy;
01669 if (is_enabled(f_cv_Jacobian)) {
01670 fj.reset();
01671 for (size_t i = 0; i < cvcs.size(); i++) {
01672 if (!cvcs[i]->is_enabled()) continue;
01673 if (cvm::debug())
01674 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
01675 " within colvar \""+this->name+"\" has Jacobian derivative"+
01676 cvm::to_str((cvcs[i])->Jacobian_derivative(),
01677 cvm::cv_width, cvm::cv_prec)+".\n");
01678 // linear combination is assumed
01679 fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
01680 }
01681 fj *= proxy->boltzmann() * proxy->target_temperature();
01682 }
01683 
01684 return COLVARS_OK;
01685 }
01686 
01687 
01688 int colvar::calc_colvar_properties()
01689 {
01690 if (is_enabled(f_cv_fdiff_velocity)) {
01691 // calculate the velocity by finite differences
01692 if (cvm::step_relative() == 0) {
01693 x_old = x;
01694 v_fdiff.reset(); // Do not pretend we know anything about the actual velocity
01695 // eg. upon restarting. That would require saving v_fdiff or x_old to the state file
01696 } else {
01697 v_fdiff = fdiff_velocity(x_old, x);
01698 v_reported = v_fdiff;
01699 }
01700 }
01701 
01702 if (is_enabled(f_cv_extended_Lagrangian)) {
01703 // initialize the restraint center in the first step to the value
01704 // just calculated from the cvcs
01705 // Do the same if no simulation is running (eg. VMD postprocessing)
01706 if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset || !cvm::proxy->simulation_running()) {
01707 x_ext = x;
01708 if (is_enabled(f_cv_reflecting_lower_boundary) && x_ext < lower_boundary) {
01709 cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.");
01710 x_ext = lower_boundary;
01711 }
01712 if (is_enabled(f_cv_reflecting_upper_boundary) && x_ext > upper_boundary) {
01713 cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.");
01714 x_ext = upper_boundary;
01715 }
01716 
01717 v_ext.reset(); // (already 0; added for clarity)
01718 }
01719 
01720 // Special case of a repeated timestep (eg. multiple NAMD "run" statements)
01721 // revert values of the extended coordinate and velocity prior to latest integration
01722 if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) {
01723 x_ext = prev_x_ext;
01724 v_ext = prev_v_ext;
01725 }
01726 // report the restraint center as "value"
01727 // These position and velocities come from integration at the _previous timestep_ in update_forces_energy()
01728 // But we report values at the beginning of the timestep (value at t=0 on the first timestep)
01729 x_reported = x_ext;
01730 v_reported = v_ext;
01731 // the "total force" with the extended Lagrangian is
01732 // calculated in update_forces_energy() below
01733 
01734 } else {
01735 
01736 if (is_enabled(f_cv_subtract_applied_force)) {
01737 // correct the total force only if it has been measured
01738 // TODO add a specific test instead of relying on sq norm
01739 if (ft.norm2() > 0.0) {
01740 ft -= f_old;
01741 }
01742 }
01743 
01744 x_reported = x;
01745 ft_reported = ft;
01746 }
01747 
01748 // At the end of the first update after a restart, we can reset the flag
01749 after_restart = false;
01750 return COLVARS_OK;
01751 }
01752 
01753 
01754 cvm::real colvar::update_forces_energy()
01755 {
01756 if (cvm::debug())
01757 cvm::log("Updating colvar \""+this->name+"\".\n");
01758 
01759 // set to zero the applied force
01760 f.type(value());
01761 f.reset();
01762 fr.reset();
01763 
01764 // If we are not active at this timestep, that's all we have to do
01765 // return with energy == zero
01766 if (!is_enabled(f_cv_active)) return 0.;
01767 
01768 // add the biases' force, which at this point should already have
01769 // been summed over each bias using this colvar
01770 // fb is already multiplied by the relevant time step factor for each bias
01771 f += fb;
01772 
01773 if (is_enabled(f_cv_Jacobian)) {
01774 // the instantaneous Jacobian force was not included in the reported total force;
01775 // instead, it is subtracted from the applied force (silent Jacobian correction)
01776 // This requires the Jacobian term for the *current* timestep
01777 // Need to scale it for impulse MTS
01778 if (is_enabled(f_cv_hide_Jacobian))
01779 f -= fj * cvm::real(time_step_factor);
01780 }
01781 
01782 // At this point f is the force f from external biases that will be applied to the
01783 // extended variable if there is one
01784 if (is_enabled(f_cv_extended_Lagrangian) && cvm::proxy->simulation_running()) {
01785 update_extended_Lagrangian();
01786 }
01787 
01788 if (!is_enabled(f_cv_external)) {
01789 // Now adding the force on the actual colvar (for those biases that
01790 // bypass the extended Lagrangian mass)
01791 f += fb_actual;
01792 }
01793 
01794 if (cvm::debug())
01795 cvm::log("Done updating colvar \""+this->name+"\".\n");
01796 return (potential_energy + kinetic_energy);
01797 }
01798 
01799 
01800 void colvar::update_extended_Lagrangian()
01801 {
01802 if (cvm::debug()) {
01803 cvm::log("Updating extended-Lagrangian degree of freedom.\n");
01804 }
01805 
01806 if (prev_timestep > -1L) {
01807 // Keep track of slow timestep to integrate MTS colvars
01808 // the colvar checks the interval after waking up twice
01809 cvm::step_number n_timesteps = cvm::step_relative() - prev_timestep;
01810 if (n_timesteps != 0 && n_timesteps != time_step_factor) {
01811 cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
01812 cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
01813 " steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " +
01814 cvm::to_str(cvm::step_relative()) + ").\n" +
01815 "Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n");
01816 return;
01817 }
01818 }
01819 
01820 // Integrate with slow timestep (if time_step_factor != 1)
01821 cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
01822 
01823 colvarvalue f_ext(fr.type()); // force acting on the extended variable
01824 f_ext.reset();
01825 
01826 if (is_enabled(f_cv_external)) {
01827 // There are no forces on the "actual colvar" bc there is no gradient wrt atomic coordinates
01828 // So we apply this to the extended DOF
01829 f += fb_actual;
01830 }
01831 
01832 fr = f;
01833 // External force has been scaled for a 1-timestep impulse, scale it back because we will
01834 // integrate it with the colvar's own timestep factor
01835 f_ext = f / cvm::real(time_step_factor);
01836 
01837 colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
01838 
01839 if (is_enabled(f_cv_external)) {
01840 // Add "alchemical" force from external variable
01841 f_system = cvcs[0]->total_force();
01842 // f is now irrelevant because we are not applying atomic forces in the simulation
01843 // just driving the external variable lambda
01844 } else {
01845 // the total force is applied to the fictitious mass, while the
01846 // atoms only feel the harmonic force + wall force
01847 // fr: bias force on extended variable (without harmonic spring), for output in trajectory
01848 // f_ext: total force on extended variable (including harmonic spring)
01849 // f: - initially, external biasing force
01850 // - after this code block, colvar force to be applied to atomic coordinates
01851 // ie. spring force (fb_actual will be added just below)
01852 f_system = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
01853 f = -1.0 * f_system;
01854 // Coupling force is a slow force, to be applied to atomic coords impulse-style
01855 // over a single MD timestep
01856 f *= cvm::real(time_step_factor);
01857 }
01858 f_ext += f_system;
01859 
01860 if (is_enabled(f_cv_subtract_applied_force)) {
01861 // Report a "system" force without the biases on this colvar
01862 // that is, just the spring force (or alchemical force)
01863 ft_reported = f_system;
01864 } else {
01865 // The total force acting on the extended variable is f_ext
01866 // This will be used in the next timestep
01867 ft_reported = f_ext;
01868 }
01869 
01870 // backup in case we need to revert this integration timestep
01871 // if the same MD timestep is re-run
01872 prev_x_ext = x_ext;
01873 prev_v_ext = v_ext;
01874 
01875 // leapfrog: starting from x_i, f_i, v_(i-1/2)
01876 v_ext += (0.5 * dt) * f_ext / ext_mass;
01877 // Because of leapfrog, kinetic energy at time i is approximate
01878 kinetic_energy = 0.5 * ext_mass * v_ext * v_ext;
01879 potential_energy = 0.5 * ext_force_k * this->dist2(x_ext, x);
01880 // leap to v_(i+1/2)
01881 if (is_enabled(f_cv_Langevin)) {
01882 v_ext -= dt * ext_gamma * v_ext;
01883 colvarvalue rnd(x);
01884 rnd.set_random();
01885 v_ext += dt * ext_sigma * rnd / ext_mass;
01886 }
01887 v_ext += (0.5 * dt) * f_ext / ext_mass;
01888 x_ext += dt * v_ext;
01889 
01890 cvm::real delta = 0; // Length of overshoot past either reflecting boundary
01891 if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
01892 (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
01893 x_ext -= 2.0 * delta;
01894 v_ext *= -1.0;
01895 if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
01896 (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
01897 cvm::error("Error: extended coordinate value " + cvm::to_str(x_ext) + " is still outside boundaries after reflection.\n");
01898 }
01899 }
01900 
01901 x_ext.apply_constraints();
01902 this->wrap(x_ext);
01903 if (is_enabled(f_cv_external)) {
01904 // Colvar value is constrained to the extended value
01905 x = x_ext;
01906 cvcs[0]->set_value(x_ext);
01907 }
01908 }
01909 
01910 
01911 int colvar::end_of_step()
01912 {
01913 if (cvm::debug())
01914 cvm::log("End of step for colvar \""+this->name+"\".\n");
01915 
01916 if (is_enabled(f_cv_fdiff_velocity)) {
01917 x_old = x;
01918 }
01919 
01920 if (is_enabled(f_cv_subtract_applied_force)) {
01921 f_old = f;
01922 }
01923 
01924 prev_timestep = cvm::step_relative();
01925 
01926 return COLVARS_OK;
01927 }
01928 
01929 
01930 void colvar::communicate_forces()
01931 {
01932 size_t i;
01933 if (cvm::debug()) {
01934 cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
01935 cvm::log("Force to be applied: " + cvm::to_str(f) + "\n");
01936 }
01937 
01938 if (is_enabled(f_cv_scripted)) {
01939 std::vector<cvm::matrix2d<cvm::real> > func_grads;
01940 func_grads.reserve(cvcs.size());
01941 for (i = 0; i < cvcs.size(); i++) {
01942 if (!cvcs[i]->is_enabled()) continue;
01943 func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
01944 cvcs[i]->value().size()));
01945 }
01946 int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
01947 
01948 if (res != COLVARS_OK) {
01949 if (res == COLVARS_NOT_IMPLEMENTED) {
01950 cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);
01951 } else {
01952 cvm::error("Error running colvar gradient script");
01953 }
01954 return;
01955 }
01956 
01957 int grad_index = 0; // index in the scripted gradients, to account for some components being disabled
01958 for (i = 0; i < cvcs.size(); i++) {
01959 if (!cvcs[i]->is_enabled()) continue;
01960 // cvc force is colvar force times colvar/cvc Jacobian
01961 // (vector-matrix product)
01962 (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
01963 cvcs[i]->value().type()));
01964 }
01965 
01966 #ifdef LEPTON
01967 } else if (is_enabled(f_cv_custom_function)) {
01968 
01969 size_t r = 0; // index in the vector of variable references
01970 size_t e = 0; // index of the gradient evaluator
01971 
01972 for (i = 0; i < cvcs.size(); i++) { // gradient with respect to cvc i
01973 cvm::matrix2d<cvm::real> jacobian (x.size(), cvcs[i]->value().size());
01974 for (size_t j = 0; j < cvcs[i]->value().size(); j++) { // j-th element
01975 for (size_t c = 0; c < x.size(); c++) { // derivative of scalar element c of the colvarvalue
01976 
01977 // Feed cvc values to the evaluator
01978 for (size_t k = 0; k < cvcs.size(); k++) { //
01979 for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
01980 *(grad_eval_var_refs[r++]) = cvcs[k]->value()[l];
01981 }
01982 }
01983 jacobian[c][j] = gradient_evaluators[e++]->evaluate();
01984 }
01985 }
01986 // cvc force is colvar force times colvar/cvc Jacobian
01987 // (vector-matrix product)
01988 (cvcs[i])->apply_force(colvarvalue(f.as_vector() * jacobian,
01989 cvcs[i]->value().type()));
01990 }
01991 #endif
01992 
01993 } else if (x.type() == colvarvalue::type_scalar) {
01994 
01995 for (i = 0; i < cvcs.size(); i++) {
01996 if (!cvcs[i]->is_enabled()) continue;
01997 (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
01998 cvm::real((cvcs[i])->sup_np) *
01999 (cvm::integer_power((cvcs[i])->value().real_value,
02000 (cvcs[i])->sup_np-1)) );
02001 }
02002 
02003 } else {
02004 
02005 for (i = 0; i < cvcs.size(); i++) {
02006 if (!cvcs[i]->is_enabled()) continue;
02007 (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
02008 }
02009 }
02010 
02011 if (cvm::debug())
02012 cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
02013 }
02014 
02015 
02016 int colvar::set_cvc_flags(std::vector<bool> const &flags)
02017 {
02018 if (flags.size() != cvcs.size()) {
02019 cvm::error("ERROR: Wrong number of CVC flags provided.");
02020 return COLVARS_ERROR;
02021 }
02022 // We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
02023 // so we store the flags that will be enforced at the next call to calc()
02024 cvc_flags = flags;
02025 return COLVARS_OK;
02026 }
02027 
02028 
02029 void colvar::update_active_cvc_square_norm()
02030 {
02031 active_cvc_square_norm = 0.0;
02032 for (size_t i = 0; i < cvcs.size(); i++) {
02033 if (cvcs[i]->is_enabled()) {
02034 active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
02035 }
02036 }
02037 }
02038 
02039 
02040 int colvar::update_cvc_flags()
02041 {
02042 // Update the enabled/disabled status of cvcs if necessary
02043 if (cvc_flags.size()) {
02044 n_active_cvcs = 0;
02045 for (size_t i = 0; i < cvcs.size(); i++) {
02046 cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
02047 if (cvcs[i]->is_enabled()) {
02048 n_active_cvcs++;
02049 }
02050 }
02051 if (!n_active_cvcs) {
02052 cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
02053 return COLVARS_ERROR;
02054 }
02055 cvc_flags.clear();
02056 
02057 update_active_cvc_square_norm();
02058 }
02059 
02060 return COLVARS_OK;
02061 }
02062 
02063 
02064 int colvar::update_cvc_config(std::vector<std::string> const &confs)
02065 {
02066 cvm::log("Updating configuration for colvar \""+name+"\"\n");
02067 
02068 if (confs.size() != cvcs.size()) {
02069 return cvm::error("Error: Wrong number of CVC config strings. "
02070 "For those CVCs that are not being changed, try passing "
02071 "an empty string.", COLVARS_INPUT_ERROR);
02072 }
02073 
02074 int error_code = COLVARS_OK;
02075 int num_changes = 0;
02076 for (size_t i = 0; i < cvcs.size(); i++) {
02077 if (confs[i].size()) {
02078 std::string conf(confs[i]);
02079 cvm::increase_depth();
02080 error_code |= cvcs[i]->colvar::cvc::init(conf);
02081 error_code |= cvcs[i]->check_keywords(conf,
02082 cvcs[i]->config_key.c_str());
02083 cvm::decrease_depth();
02084 num_changes++;
02085 }
02086 }
02087 
02088 if (num_changes == 0) {
02089 cvm::log("Warning: no changes were applied through modifycvcs; "
02090 "please check that its argument is a list of strings.\n");
02091 }
02092 
02093 update_active_cvc_square_norm();
02094 
02095 return error_code;
02096 }
02097 
02098 
02099 int colvar::cvc_param_exists(std::string const &param_name)
02100 {
02101 if (is_enabled(f_cv_single_cvc)) {
02102 return cvcs[0]->param_exists(param_name);
02103 }
02104 return cvm::error("Error: calling colvar::cvc_param_exists() for a variable "
02105 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
02106 }
02107 
02108 
02109 cvm::real colvar::get_cvc_param(std::string const &param_name)
02110 {
02111 if (is_enabled(f_cv_single_cvc)) {
02112 return cvcs[0]->get_param(param_name);
02113 }
02114 cvm::error("Error: calling colvar::get_cvc_param() for a variable "
02115 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
02116 return 0.0;
02117 }
02118 
02119 
02120 void const *colvar::get_cvc_param_ptr(std::string const &param_name)
02121 {
02122 if (is_enabled(f_cv_single_cvc)) {
02123 return cvcs[0]->get_param_ptr(param_name);
02124 }
02125 cvm::error("Error: calling colvar::get_cvc_param() for a variable "
02126 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
02127 return NULL;
02128 }
02129 
02130 
02131 colvarvalue const *colvar::get_cvc_param_grad(std::string const &param_name)
02132 {
02133 if (is_enabled(f_cv_single_cvc)) {
02134 return cvcs[0]->get_param_grad(param_name);
02135 }
02136 cvm::error("Error: calling colvar::get_cvc_param_grad() for a variable "
02137 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
02138 return NULL;
02139 }
02140 
02141 
02142 int colvar::set_cvc_param(std::string const &param_name, void const *new_value)
02143 {
02144 if (is_enabled(f_cv_single_cvc)) {
02145 return cvcs[0]->set_param(param_name, new_value);
02146 }
02147 return cvm::error("Error: calling colvar::set_cvc_param() for a variable "
02148 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
02149 }
02150 
02151 
02152 // ******************** METRIC FUNCTIONS ********************
02153 // Use the metrics defined by \link colvar::cvc \endlink objects
02154 
02155 
02156 bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
02157 {
02158 if (period > 0.0) {
02159 if ( ((cvm::sqrt(this->dist2(lb, ub))) / this->width)
02160 < 1.0E-10 ) {
02161 return true;
02162 }
02163 }
02164 
02165 return false;
02166 }
02167 
02168 bool colvar::periodic_boundaries() const
02169 {
02170 if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
02171 // Return false if answer is unknown at this time
02172 return false;
02173 }
02174 
02175 return periodic_boundaries(lower_boundary, upper_boundary);
02176 }
02177 
02178 
02179 cvm::real colvar::dist2(colvarvalue const &x1,
02180 colvarvalue const &x2) const
02181 {
02182 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
02183 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
02184 cvm::real diff = x1.real_value - x2.real_value;
02185 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
02186 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
02187 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
02188 return diff * diff;
02189 }
02190 }
02191 if (is_enabled(f_cv_homogeneous)) {
02192 return (cvcs[0])->dist2(x1, x2);
02193 } else {
02194 return x1.dist2(x2);
02195 }
02196 }
02197 
02198 colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
02199 colvarvalue const &x2) const
02200 {
02201 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
02202 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
02203 cvm::real diff = x1.real_value - x2.real_value;
02204 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
02205 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
02206 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
02207 return 2.0 * diff;
02208 }
02209 }
02210 if (is_enabled(f_cv_homogeneous)) {
02211 return (cvcs[0])->dist2_lgrad(x1, x2);
02212 } else {
02213 return x1.dist2_grad(x2);
02214 }
02215 }
02216 
02217 colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
02218 colvarvalue const &x2) const
02219 {
02220 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
02221 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
02222 cvm::real diff = x1.real_value - x2.real_value;
02223 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
02224 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
02225 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
02226 return (-2.0) * diff;
02227 }
02228 }
02229 if (is_enabled(f_cv_homogeneous)) {
02230 return (cvcs[0])->dist2_rgrad(x1, x2);
02231 } else {
02232 return x2.dist2_grad(x1);
02233 }
02234 }
02235 
02236 
02237 void colvar::wrap(colvarvalue &x_unwrapped) const
02238 {
02239 if (!is_enabled(f_cv_periodic)) {
02240 return;
02241 }
02242 
02243 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
02244 // Scripted functions do their own wrapping, as cvcs might not be periodic
02245 cvm::real shift = cvm::floor((x_unwrapped.real_value - wrap_center) /
02246 period + 0.5);
02247 x_unwrapped.real_value -= shift * period;
02248 } else {
02249 cvcs[0]->wrap(x_unwrapped);
02250 }
02251 }
02252 
02253 
02254 // ******************** INPUT FUNCTIONS ********************
02255 
02256 std::istream & colvar::read_state(std::istream &is)
02257 {
02258 std::streampos const start_pos = is.tellg();
02259 
02260 std::string conf;
02261 if ( !(is >> colvarparse::read_block("colvar", &conf)) ) {
02262 // this is not a colvar block
02263 is.clear();
02264 is.seekg(start_pos, std::ios::beg);
02265 is.setstate(std::ios::failbit);
02266 return is;
02267 }
02268 
02269 {
02270 std::string check_name = "";
02271 get_keyval(conf, "name", check_name,
02272 std::string(""), colvarparse::parse_silent);
02273 if (check_name.size() == 0) {
02274 cvm::error("Error: Collective variable in the "
02275 "restart file without any identifier.\n", COLVARS_INPUT_ERROR);
02276 is.clear();
02277 is.seekg(start_pos, std::ios::beg);
02278 is.setstate(std::ios::failbit);
02279 return is;
02280 }
02281 
02282 if (check_name != name) {
02283 if (cvm::debug()) {
02284 cvm::log("Ignoring state of colvar \""+check_name+
02285 "\": this colvar is named \""+name+"\".\n");
02286 }
02287 is.seekg(start_pos, std::ios::beg);
02288 return is;
02289 }
02290 }
02291 
02292 if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
02293 cvm::log("Error: restart file does not contain "
02294 "the value of the colvar \""+
02295 name+"\" .\n");
02296 } else {
02297 cvm::log("Restarting collective variable \""+name+"\" from value: "+
02298 cvm::to_str(x)+"\n");
02299 x_restart = x;
02300 after_restart = true;
02301 }
02302 
02303 if (is_enabled(f_cv_extended_Lagrangian)) {
02304 if ( !(get_keyval(conf, "extended_x", x_ext,
02305 colvarvalue(x.type()), colvarparse::parse_silent)) ||
02306 !(get_keyval(conf, "extended_v", v_ext,
02307 colvarvalue(x.type()), colvarparse::parse_silent)) ) {
02308 cvm::log("Error: restart file does not contain "
02309 "\"extended_x\" or \"extended_v\" for the colvar \""+
02310 name+"\", but you requested \"extendedLagrangian\".\n");
02311 }
02312 x_reported = x_ext;
02313 } else {
02314 x_reported = x;
02315 }
02316 
02317 if (is_enabled(f_cv_output_velocity)) {
02318 
02319 if ( !(get_keyval(conf, "v", v_fdiff,
02320 colvarvalue(x.type()), colvarparse::parse_silent)) ) {
02321 cvm::log("Error: restart file does not contain "
02322 "the velocity for the colvar \""+
02323 name+"\", but you requested \"outputVelocity\".\n");
02324 }
02325 
02326 if (is_enabled(f_cv_extended_Lagrangian)) {
02327 v_reported = v_ext;
02328 } else {
02329 v_reported = v_fdiff;
02330 }
02331 }
02332 
02333 return is;
02334 }
02335 
02336 
02337 std::istream & colvar::read_traj(std::istream &is)
02338 {
02339 std::streampos const start_pos = is.tellg();
02340 
02341 if (is_enabled(f_cv_output_value)) {
02342 
02343 if (!(is >> x)) {
02344 cvm::log("Error: in reading the value of colvar \""+
02345 this->name+"\" from trajectory.\n");
02346 is.clear();
02347 is.seekg(start_pos, std::ios::beg);
02348 is.setstate(std::ios::failbit);
02349 return is;
02350 }
02351 
02352 if (is_enabled(f_cv_extended_Lagrangian)) {
02353 is >> x_ext;
02354 x_reported = x_ext;
02355 } else {
02356 x_reported = x;
02357 }
02358 }
02359 
02360 if (is_enabled(f_cv_output_velocity)) {
02361 
02362 is >> v_fdiff;
02363 
02364 if (is_enabled(f_cv_extended_Lagrangian)) {
02365 is >> v_ext;
02366 v_reported = v_ext;
02367 } else {
02368 v_reported = v_fdiff;
02369 }
02370 }
02371 
02372 if (is_enabled(f_cv_output_total_force)) {
02373 is >> ft;
02374 ft_reported = ft;
02375 }
02376 
02377 if (is_enabled(f_cv_output_applied_force)) {
02378 is >> f;
02379 }
02380 
02381 return is;
02382 }
02383 
02384 
02385 // ******************** OUTPUT FUNCTIONS ********************
02386 
02387 std::ostream & colvar::write_state(std::ostream &os) {
02388 
02389 os << "colvar {\n"
02390 << " name " << name << "\n"
02391 << " x "
02392 << std::setprecision(cvm::cv_prec)
02393 << std::setw(cvm::cv_width)
02394 << x << "\n";
02395 
02396 if (is_enabled(f_cv_output_velocity)) {
02397 os << " v "
02398 << std::setprecision(cvm::cv_prec)
02399 << std::setw(cvm::cv_width)
02400 << v_reported << "\n";
02401 }
02402 
02403 if (is_enabled(f_cv_extended_Lagrangian)) {
02404 os << " extended_x "
02405 << std::setprecision(cvm::cv_prec)
02406 << std::setw(cvm::cv_width)
02407 << x_reported << "\n"
02408 << " extended_v "
02409 << std::setprecision(cvm::cv_prec)
02410 << std::setw(cvm::cv_width)
02411 << v_reported << "\n";
02412 }
02413 
02414 os << "}\n\n";
02415 
02416 if (runave_outfile.size() > 0) {
02417 cvm::main()->proxy->flush_output_stream(runave_outfile);
02418 }
02419 
02420 return os;
02421 }
02422 
02423 
02424 std::ostream & colvar::write_traj_label(std::ostream & os)
02425 {
02426 size_t const this_cv_width = x.output_width(cvm::cv_width);
02427 
02428 os << " ";
02429 
02430 if (is_enabled(f_cv_output_value)) {
02431 
02432 os << " "
02433 << cvm::wrap_string(this->name, this_cv_width);
02434 
02435 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
02436 // extended DOF
02437 os << " r_"
02438 << cvm::wrap_string(this->name, this_cv_width-2);
02439 }
02440 }
02441 
02442 if (is_enabled(f_cv_output_velocity)) {
02443 
02444 os << " v_"
02445 << cvm::wrap_string(this->name, this_cv_width-2);
02446 
02447 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
02448 // extended DOF
02449 os << " vr_"
02450 << cvm::wrap_string(this->name, this_cv_width-3);
02451 }
02452 }
02453 
02454 if (is_enabled(f_cv_output_energy)) {
02455 os << " Ep_"
02456 << cvm::wrap_string(this->name, this_cv_width-3)
02457 << " Ek_"
02458 << cvm::wrap_string(this->name, this_cv_width-3);
02459 }
02460 
02461 if (is_enabled(f_cv_output_total_force)) {
02462 os << " ft_"
02463 << cvm::wrap_string(this->name, this_cv_width-3);
02464 }
02465 
02466 if (is_enabled(f_cv_output_applied_force)) {
02467 os << " fa_"
02468 << cvm::wrap_string(this->name, this_cv_width-3);
02469 }
02470 
02471 return os;
02472 }
02473 
02474 
02475 std::ostream & colvar::write_traj(std::ostream &os)
02476 {
02477 os << " ";
02478 if (is_enabled(f_cv_output_value)) {
02479 
02480 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
02481 os << " "
02482 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02483 << x;
02484 }
02485 
02486 os << " "
02487 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02488 << x_reported;
02489 }
02490 
02491 if (is_enabled(f_cv_output_velocity)) {
02492 
02493 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
02494 os << " "
02495 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02496 << v_fdiff;
02497 }
02498 
02499 os << " "
02500 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02501 << v_reported;
02502 }
02503 
02504 if (is_enabled(f_cv_output_energy)) {
02505 os << " "
02506 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02507 << potential_energy
02508 << " "
02509 << kinetic_energy;
02510 }
02511 
02512 
02513 if (is_enabled(f_cv_output_total_force)) {
02514 os << " "
02515 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02516 << ft_reported;
02517 }
02518 
02519 if (is_enabled(f_cv_output_applied_force)) {
02520 os << " "
02521 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
02522 << applied_force();
02523 }
02524 
02525 return os;
02526 }
02527 
02528 
02529 int colvar::write_output_files()
02530 {
02531 int error_code = COLVARS_OK;
02532 
02533 if (is_enabled(f_cv_corrfunc)) {
02534 if (acf.size()) {
02535 if (acf_outfile.size() == 0) {
02536 acf_outfile = std::string(cvm::output_prefix()+"."+this->name+
02537 ".corrfunc.dat");
02538 }
02539 cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n");
02540 cvm::backup_file(acf_outfile.c_str());
02541 std::ostream &acf_os = cvm::proxy->output_stream(acf_outfile,
02542 "colvar ACF file");
02543 if (!acf_os) {
02544 error_code |= COLVARS_FILE_ERROR;
02545 } else {
02546 error_code |= write_acf(acf_os);
02547 }
02548 cvm::proxy->close_output_stream(acf_outfile);
02549 }
02550 }
02551 
02552 return error_code;
02553 }
02554 
02555 
02556 
02557 // ******************** ANALYSIS FUNCTIONS ********************
02558 
02559 int colvar::analyze()
02560 {
02561 int error_code = COLVARS_OK;
02562 
02563 if (is_enabled(f_cv_runave)) {
02564 error_code |= calc_runave();
02565 }
02566 
02567 if (is_enabled(f_cv_corrfunc)) {
02568 error_code |= calc_acf();
02569 }
02570 
02571 return error_code;
02572 }
02573 
02574 
02575 inline void history_add_value(size_t const &history_length,
02576 std::list<colvarvalue> &history,
02577 colvarvalue const &new_value)
02578 {
02579 history.push_front(new_value);
02580 if (history.size() > history_length)
02581 history.pop_back();
02582 }
02583 
02584 
02585 inline void history_incr(std::list< std::list<colvarvalue> > &history,
02586 std::list< std::list<colvarvalue> >::iterator &history_p)
02587 {
02588 if ((++history_p) == history.end())
02589 history_p = history.begin();
02590 }
02591 
02592 
02593 int colvar::calc_acf()
02594 {
02595 // using here an acf_stride-long list of vectors for either
02596 // coordinates (acf_x_history) or velocities (acf_v_history); each vector can
02597 // contain up to acf_length values, which are contiguous in memory
02598 // representation but separated by acf_stride in the time series;
02599 // the pointer to each vector is changed at every step
02600 
02601 colvar const *cfcv = cvm::colvar_by_name(acf_colvar_name);
02602 if (cfcv == NULL) {
02603 return cvm::error("Error: collective variable \""+acf_colvar_name+
02604 "\" is not defined at this time.\n", COLVARS_INPUT_ERROR);
02605 }
02606 
02607 if (acf_x_history.empty() && acf_v_history.empty()) {
02608 
02609 // first-step operations
02610 
02611 if (colvarvalue::check_types(cfcv->value(), value())) {
02612 cvm::error("Error: correlation function between \""+cfcv->name+
02613 "\" and \""+this->name+"\" cannot be calculated, "
02614 "because their value types are different.\n",
02615 COLVARS_INPUT_ERROR);
02616 }
02617 acf_nframes = 0;
02618 
02619 cvm::log("Colvar \""+this->name+"\": initializing correlation function "
02620 "calculation.\n");
02621 
02622 if (acf.size() < acf_length+1)
02623 acf.resize(acf_length+1, 0.0);
02624 
02625 size_t i;
02626 switch (acf_type) {
02627 
02628 case acf_vel:
02629 // allocate space for the velocities history
02630 for (i = 0; i < acf_stride; i++) {
02631 acf_v_history.push_back(std::list<colvarvalue>());
02632 }
02633 acf_v_history_p = acf_v_history.begin();
02634 break;
02635 
02636 case acf_coor:
02637 case acf_p2coor:
02638 // allocate space for the coordinates history
02639 for (i = 0; i < acf_stride; i++) {
02640 acf_x_history.push_back(std::list<colvarvalue>());
02641 }
02642 acf_x_history_p = acf_x_history.begin();
02643 break;
02644 
02645 case acf_notset:
02646 default:
02647 break;
02648 }
02649 
02650 } else if (cvm::step_relative() > prev_timestep) {
02651 
02652 switch (acf_type) {
02653 
02654 case acf_vel:
02655 
02656 calc_vel_acf((*acf_v_history_p), cfcv->velocity());
02657 history_add_value(acf_length+acf_offset, *acf_v_history_p,
02658 cfcv->velocity());
02659 history_incr(acf_v_history, acf_v_history_p);
02660 break;
02661 
02662 case acf_coor:
02663 
02664 calc_coor_acf((*acf_x_history_p), cfcv->value());
02665 history_add_value(acf_length+acf_offset, *acf_x_history_p,
02666 cfcv->value());
02667 history_incr(acf_x_history, acf_x_history_p);
02668 break;
02669 
02670 case acf_p2coor:
02671 
02672 calc_p2coor_acf((*acf_x_history_p), cfcv->value());
02673 history_add_value(acf_length+acf_offset, *acf_x_history_p,
02674 cfcv->value());
02675 history_incr(acf_x_history, acf_x_history_p);
02676 break;
02677 
02678 case acf_notset:
02679 default:
02680 break;
02681 }
02682 }
02683 
02684 return COLVARS_OK;
02685 }
02686 
02687 
02688 void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
02689 colvarvalue const &v)
02690 {
02691 // loop over stored velocities and add to the ACF, but only if the
02692 // length is sufficient to hold an entire row of ACF values
02693 if (v_list.size() >= acf_length+acf_offset) {
02694 std::list<colvarvalue>::iterator vs_i = v_list.begin();
02695 std::vector<cvm::real>::iterator acf_i = acf.begin();
02696 
02697 for (size_t i = 0; i < acf_offset; i++)
02698 ++vs_i;
02699 
02700 // current vel with itself
02701 *(acf_i) += v.norm2();
02702 ++acf_i;
02703 
02704 // inner products of previous velocities with current (acf_i and
02705 // vs_i are updated)
02706 colvarvalue::inner_opt(v, vs_i, v_list.end(), acf_i);
02707 
02708 acf_nframes++;
02709 }
02710 }
02711 
02712 
02713 void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
02714 colvarvalue const &x_now)
02715 {
02716 // same as above but for coordinates
02717 if (x_list.size() >= acf_length+acf_offset) {
02718 std::list<colvarvalue>::iterator xs_i = x_list.begin();
02719 std::vector<cvm::real>::iterator acf_i = acf.begin();
02720 
02721 for (size_t i = 0; i < acf_offset; i++)
02722 ++xs_i;
02723 
02724 *(acf_i++) += x.norm2();
02725 
02726 colvarvalue::inner_opt(x_now, xs_i, x_list.end(), acf_i);
02727 
02728 acf_nframes++;
02729 }
02730 }
02731 
02732 
02733 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
02734 colvarvalue const &x_now)
02735 {
02736 // same as above but with second order Legendre polynomial instead
02737 // of just the scalar product
02738 if (x_list.size() >= acf_length+acf_offset) {
02739 std::list<colvarvalue>::iterator xs_i = x_list.begin();
02740 std::vector<cvm::real>::iterator acf_i = acf.begin();
02741 
02742 for (size_t i = 0; i < acf_offset; i++)
02743 ++xs_i;
02744 
02745 // value of P2(0) = 1
02746 *(acf_i++) += 1.0;
02747 
02748 colvarvalue::p2leg_opt(x_now, xs_i, x_list.end(), acf_i);
02749 
02750 acf_nframes++;
02751 }
02752 }
02753 
02754 
02755 int colvar::write_acf(std::ostream &os)
02756 {
02757 if (!acf_nframes) {
02758 return COLVARS_OK;
02759 }
02760 
02761 os.setf(std::ios::scientific, std::ios::floatfield);
02762 os << "# ";
02763 switch (acf_type) {
02764 case acf_vel:
02765 os << "Velocity";
02766 break;
02767 case acf_coor:
02768 os << "Coordinate";
02769 break;
02770 case acf_p2coor:
02771 os << "Coordinate (2nd Legendre poly)";
02772 break;
02773 case acf_notset:
02774 default:
02775 break;
02776 }
02777 
02778 if (acf_colvar_name == name) {
02779 os << " autocorrelation function for variable \""
02780 << this->name << "\"\n";
02781 } else {
02782 os << " correlation function between variables \"" //
02783 << this->name << "\" and \"" << acf_colvar_name << "\"\n";
02784 }
02785 
02786 os << "# Number of samples = ";
02787 if (acf_normalize) {
02788 os << (acf_nframes-1) << " (one DoF is used for normalization)\n";
02789 } else {
02790 os << acf_nframes << "\n";
02791 }
02792 
02793 os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "
02794 << cvm::wrap_string("corrfunc(step)", cvm::cv_width) << "\n";
02795 
02796 cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
02797 
02798 std::vector<cvm::real>::iterator acf_i;
02799 size_t it = acf_offset;
02800 for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
02801 os << std::setw(cvm::it_width) << acf_stride * (it++) << " "
02802 << std::setprecision(cvm::cv_prec)
02803 << std::setw(cvm::cv_width)
02804 << ( acf_normalize ?
02805 (*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
02806 (*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
02807 }
02808 
02809 return os.good() ? COLVARS_OK : COLVARS_FILE_ERROR;
02810 }
02811 
02812 
02813 int colvar::calc_runave()
02814 {
02815 int error_code = COLVARS_OK;
02816 colvarproxy *proxy = cvm::main()->proxy;
02817 
02818 if (x_history.empty()) {
02819 
02820 runave.type(value().type());
02821 runave.reset();
02822 
02823 // first-step operationsf
02824 
02825 if (cvm::debug())
02826 cvm::log("Colvar \""+this->name+
02827 "\": initializing running average calculation.\n");
02828 
02829 acf_nframes = 0;
02830 
02831 x_history.push_back(std::list<colvarvalue>());
02832 x_history_p = x_history.begin();
02833 
02834 } else {
02835 
02836 if ( (cvm::step_relative() % runave_stride) == 0 &&
02837 (cvm::step_relative() > prev_timestep) ) {
02838 
02839 if ((*x_history_p).size() >= runave_length-1) {
02840 
02841 if (runave_outfile.size() == 0) {
02842 runave_outfile = std::string(cvm::output_prefix()+"."+
02843 this->name+".runave.traj");
02844 }
02845 
02846 if (! proxy->output_stream_exists(runave_outfile)) {
02847 size_t const this_cv_width = x.output_width(cvm::cv_width);
02848 std::ostream &runave_os = proxy->output_stream(runave_outfile,
02849 "colvar running average");
02850 runave_os.setf(std::ios::scientific, std::ios::floatfield);
02851 runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
02852 << " "
02853 << cvm::wrap_string("running average", this_cv_width)
02854 << " "
02855 << cvm::wrap_string("running stddev", this_cv_width)
02856 << "\n";
02857 }
02858 
02859 runave = x;
02860 std::list<colvarvalue>::iterator xs_i;
02861 for (xs_i = (*x_history_p).begin();
02862 xs_i != (*x_history_p).end(); ++xs_i) {
02863 runave += (*xs_i);
02864 }
02865 runave *= 1.0 / cvm::real(runave_length);
02866 runave.apply_constraints();
02867 
02868 runave_variance = 0.0;
02869 runave_variance += this->dist2(x, runave);
02870 for (xs_i = (*x_history_p).begin();
02871 xs_i != (*x_history_p).end(); ++xs_i) {
02872 runave_variance += this->dist2(x, (*xs_i));
02873 }
02874 runave_variance *= 1.0 / cvm::real(runave_length-1);
02875 
02876 if (runave_outfile.size() > 0) {
02877 std::ostream &runave_os = proxy->output_stream(runave_outfile);
02878 runave_os << std::setw(cvm::it_width) << cvm::step_relative()
02879 << " "
02880 << std::setprecision(cvm::cv_prec)
02881 << std::setw(cvm::cv_width)
02882 << runave << " "
02883 << std::setprecision(cvm::cv_prec)
02884 << std::setw(cvm::cv_width)
02885 << cvm::sqrt(runave_variance) << "\n";
02886 }
02887 }
02888 
02889 history_add_value(runave_length, *x_history_p, x);
02890 }
02891 }
02892 
02893 return error_code;
02894 }
02895 
02896 // Static members
02897 
02898 std::vector<colvardeps::feature *> colvar::cv_features;