1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/Colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #include <list>
11 #include <vector>
12 #include <algorithm>
13
14 #include "colvarmodule.h"
15 #include "colvarvalue.h"
16 #include "colvarparse.h"
17 #include "colvar.h"
18 #include "colvarcomp.h"
19 #include "colvarscript.h"
20
21 #if (__cplusplus >= 201103L)
22 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)>>();
23 #endif
24
colvar()25 colvar::colvar()
26 {
27 runave_os = NULL;
28
29 prev_timestep = -1L;
30 after_restart = false;
31 kinetic_energy = 0.0;
32 potential_energy = 0.0;
33
34 #ifdef LEPTON
35 dev_null = 0.0;
36 #endif
37
38 expand_boundaries = false;
39
40 description = "uninitialized colvar";
41 init_dependencies();
42 }
43
44
45 namespace {
46 /// Compare two cvcs using their names
47 /// Used to sort CVC array in scripted coordinates
compare(colvar::cvc * i,colvar::cvc * j)48 bool compare(colvar::cvc *i, colvar::cvc *j)
49 {
50 return i->name < j->name;
51 }
52 }
53
54
init(std::string const & conf)55 int colvar::init(std::string const &conf)
56 {
57 cvm::log("Initializing a new collective variable.\n");
58 colvarparse::init(conf);
59
60 int error_code = COLVARS_OK;
61
62 colvarmodule *cv = cvm::main();
63
64 get_keyval(conf, "name", this->name,
65 (std::string("colvar")+cvm::to_str(cv->variables()->size())));
66
67 if ((cvm::colvar_by_name(this->name) != NULL) &&
68 (cvm::colvar_by_name(this->name) != this)) {
69 cvm::error("Error: this colvar cannot have the same name, \""+this->name+
70 "\", as another colvar.\n",
71 INPUT_ERROR);
72 return INPUT_ERROR;
73 }
74
75 // Initialize dependency members
76 // Could be a function defined in a different source file, for space?
77
78 this->description = "colvar " + this->name;
79
80 error_code |= init_components(conf);
81 if (error_code != COLVARS_OK) {
82 return cvm::get_error();
83 }
84
85 size_t i;
86
87 #ifdef LEPTON
88 error_code |= init_custom_function(conf);
89 if (error_code != COLVARS_OK) {
90 return cvm::get_error();
91 }
92 #endif
93
94 // Setup colvar as scripted function of components
95 if (get_keyval(conf, "scriptedFunction", scripted_function,
96 "", colvarparse::parse_silent)) {
97
98 enable(f_cv_scripted);
99 cvm::log("This colvar uses scripted function \"" + scripted_function + "\".\n");
100
101 std::string type_str;
102 get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
103
104 x.type(colvarvalue::type_notset);
105 int t;
106 for (t = 0; t < colvarvalue::type_all; t++) {
107 if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
108 x.type(colvarvalue::Type(t));
109 break;
110 }
111 }
112 if (x.type() == colvarvalue::type_notset) {
113 cvm::error("Could not parse scripted colvar type.", INPUT_ERROR);
114 return INPUT_ERROR;
115 }
116
117 cvm::log(std::string("Expecting colvar value of type ")
118 + colvarvalue::type_desc(x.type()));
119
120 if (x.type() == colvarvalue::type_vector) {
121 int size;
122 if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
123 cvm::error("Error: no size specified for vector scripted function.",
124 INPUT_ERROR);
125 return INPUT_ERROR;
126 }
127 x.vector1d_value.resize(size);
128 }
129
130 x_reported.type(x);
131
132 // Sort array of cvcs based on their names
133 // Note: default CVC names are in input order for same type of CVC
134 std::sort(cvcs.begin(), cvcs.end(), compare);
135
136 if(cvcs.size() > 1) {
137 cvm::log("Sorted list of components for this scripted colvar:\n");
138 for (i = 0; i < cvcs.size(); i++) {
139 cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name);
140 }
141 }
142
143 // Build ordered list of component values that will be
144 // passed to the script
145 for (i = 0; i < cvcs.size(); i++) {
146 sorted_cvc_values.push_back(&(cvcs[i]->value()));
147 }
148 }
149
150 if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function))) {
151 colvarvalue const &cvc_value = (cvcs[0])->value();
152 if (cvm::debug())
153 cvm::log ("This collective variable is a "+
154 colvarvalue::type_desc(cvc_value.type())+
155 ((cvc_value.size() > 1) ? " with "+
156 cvm::to_str(cvc_value.size())+" individual components.\n" :
157 ".\n"));
158 x.type(cvc_value);
159 x_reported.type(cvc_value);
160 }
161
162 set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
163
164 // If using scripted biases, any colvar may receive bias forces
165 // and will need its gradient
166 if (cvm::scripted_forces()) {
167 enable(f_cv_gradient);
168 }
169
170 // check for linear combinations
171 {
172 bool lin = !(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function));
173 for (i = 0; i < cvcs.size(); i++) {
174
175 // FIXME this is a reverse dependency, ie. cv feature depends on cvc flag
176 // need to clarify this case
177 // if ((cvcs[i])->b_debug_gradients)
178 // enable(task_gradients);
179
180 if ((cvcs[i])->sup_np != 1) {
181 if (cvm::debug() && lin)
182 cvm::log("Warning: You are using a non-linear polynomial "
183 "combination to define this collective variable, "
184 "some biasing methods may be unavailable.\n");
185 lin = false;
186
187 if ((cvcs[i])->sup_np < 0) {
188 cvm::log("Warning: you chose a negative exponent in the combination; "
189 "if you apply forces, the simulation may become unstable "
190 "when the component \""+
191 (cvcs[i])->function_type+"\" approaches zero.\n");
192 }
193 }
194 }
195 set_enabled(f_cv_linear, lin);
196 }
197
198 // Colvar is homogeneous if:
199 // - it is linear (hence not scripted)
200 // - all cvcs have coefficient 1 or -1
201 // i.e. sum or difference of cvcs
202 {
203 bool homogeneous = is_enabled(f_cv_linear);
204 for (i = 0; i < cvcs.size(); i++) {
205 if (cvm::fabs(cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {
206 homogeneous = false;
207 }
208 }
209 set_enabled(f_cv_homogeneous, homogeneous);
210 }
211
212 // A single-component variable almost concides with its CVC object
213 if ((cvcs.size() == 1) && is_enabled(f_cv_homogeneous)) {
214 if ( !is_enabled(f_cv_scripted) && !is_enabled(f_cv_custom_function) &&
215 (cvm::fabs(cvcs[0]->sup_coeff - 1.0) < 1.0e-10) &&
216 (cvcs[0]->sup_np == 1) ) {
217 enable(f_cv_single_cvc);
218 }
219 }
220
221 // Colvar is deemed periodic if:
222 // - it is homogeneous
223 // - all cvcs are periodic
224 // - all cvcs have the same period
225 if (is_enabled(f_cv_homogeneous) && cvcs[0]->is_enabled(f_cvc_periodic)) {
226 bool b_periodic = true;
227 period = cvcs[0]->period;
228 wrap_center = cvcs[0]->wrap_center;
229 for (i = 1; i < cvcs.size(); i++) {
230 if (!cvcs[i]->is_enabled(f_cvc_periodic) || cvcs[i]->period != period) {
231 b_periodic = false;
232 period = 0.0;
233 cvm::log("Warning: although one component is periodic, this colvar will "
234 "not be treated as periodic, either because the exponent is not "
235 "1, or because components of different periodicity are defined. "
236 "Make sure that you know what you are doing!");
237 }
238 }
239 set_enabled(f_cv_periodic, b_periodic);
240 }
241
242 // Allow scripted/custom functions to be defined as periodic
243 if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) {
244 if (get_keyval(conf, "period", period, 0.)) {
245 enable(f_cv_periodic);
246 get_keyval(conf, "wrapAround", wrap_center, 0.);
247 }
248 }
249
250 // check that cvcs are compatible
251
252 for (i = 0; i < cvcs.size(); i++) {
253
254 // components may have different types only for scripted functions
255 if (!(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && (colvarvalue::check_types(cvcs[i]->value(),
256 cvcs[0]->value())) ) {
257 cvm::error("ERROR: you are defining this collective variable "
258 "by using components of different types. "
259 "You must use the same type in order to "
260 "sum them together.\n", INPUT_ERROR);
261 return INPUT_ERROR;
262 }
263 }
264
265 active_cvc_square_norm = 0.;
266 for (i = 0; i < cvcs.size(); i++) {
267 active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
268 }
269
270 // at this point, the colvar's type is defined
271 f.type(value());
272
273 x_old.type(value());
274 v_fdiff.type(value());
275 v_reported.type(value());
276 fj.type(value());
277 ft.type(value());
278 ft_reported.type(value());
279 f_old.type(value());
280 f_old.reset();
281
282 x_restart.type(value());
283
284 reset_bias_force();
285
286 get_keyval(conf, "timeStepFactor", time_step_factor, 1);
287 if (time_step_factor < 0) {
288 cvm::error("Error: timeStepFactor must be positive.\n");
289 return COLVARS_ERROR;
290 }
291 if (time_step_factor != 1) {
292 enable(f_cv_multiple_ts);
293 }
294
295 // TODO use here information from the CVCs' own natural boundaries
296 error_code |= init_grid_parameters(conf);
297
298 error_code |= init_extended_Lagrangian(conf);
299 error_code |= init_output_flags(conf);
300
301 // Detect if we have one component that is an alchemical lambda
302 if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alch_lambda") {
303 enable(f_cv_external);
304 }
305
306 // Now that the children are defined we can solve dependencies
307 enable(f_cv_active);
308
309 error_code |= parse_analysis(conf);
310
311 if (cvm::debug())
312 cvm::log("Done initializing collective variable \""+this->name+"\".\n");
313
314 return error_code;
315 }
316
317
318 #ifdef LEPTON
init_custom_function(std::string const & conf)319 int colvar::init_custom_function(std::string const &conf)
320 {
321 std::string expr, expr_in; // expr_in is a buffer to remember expr after unsuccessful parsing
322 std::vector<Lepton::ParsedExpression> pexprs;
323 Lepton::ParsedExpression pexpr;
324 size_t pos = 0; // current position in config string
325 double *ref;
326
327 if (!key_lookup(conf, "customFunction", &expr_in, &pos)) {
328 return COLVARS_OK;
329 }
330
331 enable(f_cv_custom_function);
332 cvm::log("This colvar uses a custom function.\n");
333
334 do {
335 expr = expr_in;
336 if (cvm::debug())
337 cvm::log("Parsing expression \"" + expr + "\".\n");
338 try {
339 pexpr = Lepton::Parser::parse(expr);
340 pexprs.push_back(pexpr);
341 }
342 catch (...) {
343 cvm::error("Error parsing expression \"" + expr + "\".\n", INPUT_ERROR);
344 return INPUT_ERROR;
345 }
346
347 try {
348 value_evaluators.push_back(
349 new Lepton::CompiledExpression(pexpr.createCompiledExpression()));
350 // Define variables for cvc values
351 // Stored in order: expr1, cvc1, cvc2, expr2, cvc1...
352 for (size_t i = 0; i < cvcs.size(); i++) {
353 for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
354 std::string vn = cvcs[i]->name +
355 (cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
356 try {
357 ref =&value_evaluators.back()->getVariableReference(vn);
358 }
359 catch (...) { // Variable is absent from expression
360 // To keep the same workflow, we use a pointer to a double here
361 // that will receive CVC values - even though none was allocated by Lepton
362 ref = &dev_null;
363 cvm::log("Warning: Variable " + vn + " is absent from expression \"" + expr + "\".\n");
364 }
365 value_eval_var_refs.push_back(ref);
366 }
367 }
368 }
369 catch (...) {
370 cvm::error("Error compiling expression \"" + expr + "\".\n", INPUT_ERROR);
371 return INPUT_ERROR;
372 }
373 } while (key_lookup(conf, "customFunction", &expr_in, &pos));
374
375
376 // Now define derivative with respect to each scalar sub-component
377 for (size_t i = 0; i < cvcs.size(); i++) {
378 for (size_t j = 0; j < cvcs[i]->value().size(); j++) {
379 std::string vn = cvcs[i]->name +
380 (cvcs[i]->value().size() > 1 ? cvm::to_str(j+1) : "");
381 // Element ordering: we want the
382 // gradient vector of derivatives of all elements of the colvar
383 // wrt to a given element of a cvc ([i][j])
384 for (size_t c = 0; c < pexprs.size(); c++) {
385 gradient_evaluators.push_back(
386 new Lepton::CompiledExpression(pexprs[c].differentiate(vn).createCompiledExpression()));
387 // and record the refs to each variable in those expressions
388 for (size_t k = 0; k < cvcs.size(); k++) {
389 for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
390 std::string vvn = cvcs[k]->name +
391 (cvcs[k]->value().size() > 1 ? cvm::to_str(l+1) : "");
392 try {
393 ref = &gradient_evaluators.back()->getVariableReference(vvn);
394 }
395 catch (...) { // Variable is absent from derivative
396 // To keep the same workflow, we use a pointer to a double here
397 // that will receive CVC values - even though none was allocated by Lepton
398 cvm::log("Warning: Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n");
399 ref = &dev_null;
400 }
401 grad_eval_var_refs.push_back(ref);
402 }
403 }
404 }
405 }
406 }
407
408
409 if (value_evaluators.size() == 0) {
410 cvm::error("Error: no custom function defined.\n", INPUT_ERROR);
411 return INPUT_ERROR;
412 }
413
414 std::string type_str;
415 bool b_type_specified = get_keyval(conf, "customFunctionType",
416 type_str, "scalar", parse_silent);
417 x.type(colvarvalue::type_notset);
418 int t;
419 for (t = 0; t < colvarvalue::type_all; t++) {
420 if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
421 x.type(colvarvalue::Type(t));
422 break;
423 }
424 }
425 if (x.type() == colvarvalue::type_notset) {
426 cvm::error("Could not parse custom colvar type.", INPUT_ERROR);
427 return INPUT_ERROR;
428 }
429
430 // Guess type based on number of expressions
431 if (!b_type_specified) {
432 if (value_evaluators.size() == 1) {
433 x.type(colvarvalue::type_scalar);
434 } else {
435 x.type(colvarvalue::type_vector);
436 }
437 }
438
439 if (x.type() == colvarvalue::type_vector) {
440 x.vector1d_value.resize(value_evaluators.size());
441 }
442
443 x_reported.type(x);
444 cvm::log(std::string("Expecting colvar value of type ")
445 + colvarvalue::type_desc(x.type())
446 + (x.type()==colvarvalue::type_vector ? " of size " + cvm::to_str(x.size()) : "")
447 + ".\n");
448
449 if (x.size() != value_evaluators.size()) {
450 cvm::error("Error: based on custom function type, expected "
451 + cvm::to_str(x.size()) + " scalar expressions, but "
452 + cvm::to_str(value_evaluators.size()) + " were found.\n");
453 return INPUT_ERROR;
454 }
455
456 return COLVARS_OK;
457 }
458
459 #else
460
init_custom_function(std::string const & conf)461 int colvar::init_custom_function(std::string const &conf)
462 {
463
464 std::string expr;
465 size_t pos = 0;
466 if (key_lookup(conf, "customFunction", &expr, &pos)) {
467 std::string msg("Error: customFunction requires the Lepton library.");
468 #if (__cplusplus < 201103L)
469 // NOTE: this is not ideal; testing for the Lepton library's version would
470 // be more accurate, but also less portable
471 msg +=
472 std::string(" Note also that recent versions of Lepton require C++11: "
473 "please see https://colvars.github.io/README-c++11.html.");
474 #endif
475 return cvm::error(msg, COLVARS_NOT_IMPLEMENTED);
476 }
477
478 return COLVARS_OK;
479 }
480
481 #endif // #ifdef LEPTON
482
483
init_grid_parameters(std::string const & conf)484 int colvar::init_grid_parameters(std::string const &conf)
485 {
486 int error_code = COLVARS_OK;
487
488 colvarmodule *cv = cvm::main();
489
490 cvm::real default_width = width;
491
492 if (!key_already_set("width")) {
493 // The first time, check if the CVC has a width to provide
494 default_width = 1.0;
495 if (is_enabled(f_cv_single_cvc) && cvcs[0]->is_enabled(f_cvc_width)) {
496 cvm::real const cvc_width = cvcs[0]->get_param("width");
497 default_width = cvc_width;
498 }
499 }
500
501 get_keyval(conf, "width", width, default_width);
502
503 if (width <= 0.0) {
504 cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
505 return INPUT_ERROR;
506 }
507
508 lower_boundary.type(value());
509 upper_boundary.type(value());
510
511 if (is_enabled(f_cv_scalar)) {
512
513 if (is_enabled(f_cv_single_cvc)) {
514 // Get the default boundaries from the component
515 if (cvcs[0]->is_enabled(f_cvc_lower_boundary)) {
516 enable(f_cv_lower_boundary);
517 enable(f_cv_hard_lower_boundary);
518 lower_boundary =
519 *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("lowerBoundary")));
520 }
521 if (cvcs[0]->is_enabled(f_cvc_upper_boundary)) {
522 enable(f_cv_upper_boundary);
523 enable(f_cv_hard_upper_boundary);
524 upper_boundary =
525 *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("upperBoundary")));
526 }
527 }
528
529 if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
530 enable(f_cv_lower_boundary);
531 // Because this is the user's choice, we cannot assume it is a true
532 // physical boundary
533 disable(f_cv_hard_lower_boundary);
534 }
535
536 if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
537 enable(f_cv_upper_boundary);
538 disable(f_cv_hard_upper_boundary);
539 }
540
541 // Parse legacy wall options and set up a harmonicWalls bias if needed
542 cvm::real lower_wall_k = 0.0, upper_wall_k = 0.0;
543 cvm::real lower_wall = 0.0, upper_wall = 0.0;
544 std::string lw_conf, uw_conf;
545
546 if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
547 parse_silent)) {
548 cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
549 "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
550 if (!get_keyval(conf, "lowerWall", lower_wall)) {
551 error_code |= cvm::error("Error: the value of lowerWall must be set "
552 "explicitly.\n", INPUT_ERROR);
553 }
554 lw_conf = std::string("\n\
555 lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
556 lowerWalls "+cvm::to_str(lower_wall)+"\n");
557 }
558
559 if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
560 parse_silent)) {
561 cvm::log("Reading legacy options upperWall and upperWallConstant: "
562 "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
563 if (!get_keyval(conf, "upperWall", upper_wall)) {
564 error_code |= cvm::error("Error: the value of upperWall must be set "
565 "explicitly.\n", INPUT_ERROR);
566 }
567 uw_conf = std::string("\n\
568 upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
569 upperWalls "+cvm::to_str(upper_wall)+"\n");
570 }
571
572 if (lw_conf.size() && uw_conf.size()) {
573 if (lower_wall >= upper_wall) {
574 error_code |= cvm::error("Error: the upper wall, "+
575 cvm::to_str(upper_wall)+
576 ", is not higher than the lower wall, "+
577 cvm::to_str(lower_wall)+".\n",
578 INPUT_ERROR);
579 }
580 }
581
582 if (lw_conf.size() || uw_conf.size()) {
583 cvm::log("Generating a new harmonicWalls bias for compatibility purposes.\n");
584 std::string const walls_conf("\n\
585 harmonicWalls {\n\
586 name "+this->name+"w\n\
587 colvars "+this->name+"\n"+lw_conf+uw_conf+"\
588 timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+
589 "}\n");
590 error_code |= cv->append_new_config(walls_conf);
591 }
592 }
593
594 get_keyval_feature(this, conf, "hardLowerBoundary", f_cv_hard_lower_boundary,
595 is_enabled(f_cv_hard_lower_boundary));
596
597 get_keyval_feature(this, conf, "hardUpperBoundary", f_cv_hard_upper_boundary,
598 is_enabled(f_cv_hard_upper_boundary));
599
600 // consistency checks for boundaries and walls
601 if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) {
602 if (lower_boundary >= upper_boundary) {
603 error_code |= cvm::error("Error: the upper boundary, "+
604 cvm::to_str(upper_boundary)+
605 ", is not higher than the lower boundary, "+
606 cvm::to_str(lower_boundary)+".\n",
607 INPUT_ERROR);
608 }
609 }
610
611 get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries);
612 if (expand_boundaries && periodic_boundaries()) {
613 error_code |= cvm::error("Error: trying to expand boundaries that already "
614 "cover a whole period of a periodic colvar.\n",
615 INPUT_ERROR);
616 }
617
618 if (expand_boundaries && is_enabled(f_cv_hard_lower_boundary) &&
619 is_enabled(f_cv_hard_upper_boundary)) {
620 error_code |= cvm::error("Error: inconsistent configuration "
621 "(trying to expand boundaries, but both "
622 "hardLowerBoundary and hardUpperBoundary "
623 "are enabled).\n", INPUT_ERROR);
624 }
625
626 return error_code;
627 }
628
629
init_extended_Lagrangian(std::string const & conf)630 int colvar::init_extended_Lagrangian(std::string const &conf)
631 {
632 get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false);
633
634 if (is_enabled(f_cv_extended_Lagrangian)) {
635 cvm::real temp, tolerance, extended_period;
636
637 cvm::log("Enabling the extended Lagrangian term for colvar \""+
638 this->name+"\".\n");
639
640 // Mark x_ext as uninitialized so we can initialize it to the colvar value when updating
641 x_ext.type(colvarvalue::type_notset);
642 v_ext.type(value());
643 fr.type(value());
644 const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
645 if (temp <= 0.0) {
646 if (found)
647 cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
648 else
649 cvm::error("Error: a positive temperature must be provided, either "
650 "by enabling a thermostat, or through \"extendedTemp\".\n",
651 INPUT_ERROR);
652 return INPUT_ERROR;
653 }
654
655 get_keyval(conf, "extendedFluctuation", tolerance);
656 if (tolerance <= 0.0) {
657 cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
658 return INPUT_ERROR;
659 }
660 ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
661 cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2\n");
662
663 get_keyval(conf, "extendedTimeConstant", extended_period, 200.0);
664 if (extended_period <= 0.0) {
665 cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
666 }
667 ext_mass = (cvm::boltzmann() * temp * extended_period * extended_period)
668 / (4.0 * PI * PI * tolerance * tolerance);
669 cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)\n");
670
671 {
672 bool b_output_energy;
673 get_keyval(conf, "outputEnergy", b_output_energy, false);
674 if (b_output_energy) {
675 enable(f_cv_output_energy);
676 }
677 }
678
679 get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
680 if (ext_gamma < 0.0) {
681 cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
682 return INPUT_ERROR;
683 }
684 if (ext_gamma != 0.0) {
685 enable(f_cv_Langevin);
686 ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs
687 // Adjust Langevin sigma for slow time step if time_step_factor != 1
688 ext_sigma = cvm::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
689 }
690
691 get_keyval_feature(this, conf, "reflectingLowerBoundary", f_cv_reflecting_lower_boundary, false);
692 get_keyval_feature(this, conf, "reflectingUpperBoundary", f_cv_reflecting_upper_boundary, false);
693 }
694
695 return COLVARS_OK;
696 }
697
698
init_output_flags(std::string const & conf)699 int colvar::init_output_flags(std::string const &conf)
700 {
701 {
702 bool b_output_value;
703 get_keyval(conf, "outputValue", b_output_value, true);
704 if (b_output_value) {
705 enable(f_cv_output_value);
706 }
707 }
708
709 {
710 bool b_output_velocity;
711 get_keyval(conf, "outputVelocity", b_output_velocity, false);
712 if (b_output_velocity) {
713 enable(f_cv_output_velocity);
714 }
715 }
716
717 {
718 bool temp;
719 if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
720 cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
721 "The two are NOT identical: see https://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
722 return INPUT_ERROR;
723 }
724 }
725
726 get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
727 get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
728 get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
729
730 return COLVARS_OK;
731 }
732
733 // read the configuration and set up corresponding instances, for
734 // each type of component implemented
init_components_type(std::string const & conf,char const *,char const * def_config_key)735 template<typename def_class_name> int colvar::init_components_type(std::string const &conf,
736 char const * /* def_desc */,
737 char const *def_config_key)
738 {
739 #if (__cplusplus >= 201103L)
740 global_cvc_map[def_config_key] = [](const std::string& cvc_conf){return new def_class_name(cvc_conf);};
741 #endif
742 size_t def_count = 0;
743 std::string def_conf = "";
744 size_t pos = 0;
745 while ( this->key_lookup(conf,
746 def_config_key,
747 &def_conf,
748 &pos) ) {
749 if (!def_conf.size()) continue;
750 cvm::log("Initializing "
751 "a new \""+std::string(def_config_key)+"\" component"+
752 (cvm::debug() ? ", with configuration:\n"+def_conf
753 : ".\n"));
754 cvm::increase_depth();
755 cvc *cvcp = new def_class_name(def_conf);
756 if (cvcp != NULL) {
757 cvcs.push_back(cvcp);
758 cvcp->check_keywords(def_conf, def_config_key);
759 cvcp->config_key = def_config_key;
760 if (cvm::get_error()) {
761 cvm::error("Error: in setting up component \""+
762 std::string(def_config_key)+"\".\n", INPUT_ERROR);
763 return INPUT_ERROR;
764 }
765 cvm::decrease_depth();
766 } else {
767 cvm::error("Error: in allocating component \""+
768 std::string(def_config_key)+"\".\n",
769 MEMORY_ERROR);
770 return MEMORY_ERROR;
771 }
772
773 if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) {
774 if ( (cvcp->function_type != std::string("distance_z")) &&
775 (cvcp->function_type != std::string("dihedral")) &&
776 (cvcp->function_type != std::string("polar_phi")) &&
777 (cvcp->function_type != std::string("spin_angle")) &&
778 (cvcp->function_type != std::string("euler_phi")) &&
779 (cvcp->function_type != std::string("euler_psi"))) {
780 cvm::error("Error: invalid use of period and/or "
781 "wrapAround in a \""+
782 std::string(def_config_key)+
783 "\" component.\n"+
784 "Period: "+cvm::to_str(cvcp->period) +
785 " wrapAround: "+cvm::to_str(cvcp->wrap_center),
786 INPUT_ERROR);
787 return INPUT_ERROR;
788 }
789 }
790
791 if ( ! cvcs.back()->name.size()) {
792 std::ostringstream s;
793 s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;
794 cvcs.back()->name = s.str();
795 /* pad cvc number for correct ordering when sorting by name */
796 }
797
798 cvcs.back()->setup();
799 if (cvm::debug()) {
800 cvm::log("Done initializing a \""+
801 std::string(def_config_key)+
802 "\" component"+
803 (cvm::debug() ?
804 ", named \""+cvcs.back()->name+"\""
805 : "")+".\n");
806 }
807 def_conf = "";
808 if (cvm::debug()) {
809 cvm::log("Parsed "+cvm::to_str(cvcs.size())+
810 " components at this time.\n");
811 }
812 }
813
814 return COLVARS_OK;
815 }
816
817
init_components(std::string const & conf)818 int colvar::init_components(std::string const &conf)
819 {
820 int error_code = COLVARS_OK;
821 size_t i = 0, j = 0;
822
823 error_code |= init_components_type<distance>(conf, "distance", "distance");
824 error_code |= init_components_type<distance_vec>(conf, "distance vector", "distanceVec");
825 error_code |= init_components_type<cartesian>(conf, "Cartesian coordinates", "cartesian");
826 error_code |= init_components_type<distance_dir>(conf, "distance vector "
827 "direction", "distanceDir");
828 error_code |= init_components_type<distance_z>(conf, "distance projection "
829 "on an axis", "distanceZ");
830 error_code |= init_components_type<distance_xy>(conf, "distance projection "
831 "on a plane", "distanceXY");
832 error_code |= init_components_type<polar_theta>(conf, "spherical polar angle theta",
833 "polarTheta");
834 error_code |= init_components_type<polar_phi>(conf, "spherical azimuthal angle phi",
835 "polarPhi");
836 error_code |= init_components_type<distance_inv>(conf, "average distance "
837 "weighted by inverse power", "distanceInv");
838 error_code |= init_components_type<distance_pairs>(conf, "N1xN2-long vector "
839 "of pairwise distances", "distancePairs");
840 error_code |= init_components_type<dipole_magnitude>(conf, "dipole magnitude",
841 "dipoleMagnitude");
842 error_code |= init_components_type<coordnum>(conf, "coordination "
843 "number", "coordNum");
844 error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
845 "number", "selfCoordNum");
846 error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
847 "number", "groupCoord");
848 error_code |= init_components_type<angle>(conf, "angle", "angle");
849 error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
850 error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
851 error_code |= init_components_type<h_bond>(conf, "hydrogen bond", "hBond");
852 error_code |= init_components_type<alpha_angles>(conf, "alpha helix", "alpha");
853 error_code |= init_components_type<dihedPC>(conf, "dihedral "
854 "principal component", "dihedralPC");
855 error_code |= init_components_type<orientation>(conf, "orientation", "orientation");
856 error_code |= init_components_type<orientation_angle>(conf, "orientation "
857 "angle", "orientationAngle");
858 error_code |= init_components_type<orientation_proj>(conf, "orientation "
859 "projection", "orientationProj");
860 error_code |= init_components_type<tilt>(conf, "tilt", "tilt");
861 error_code |= init_components_type<spin_angle>(conf, "spin angle", "spinAngle");
862 error_code |= init_components_type<rmsd>(conf, "RMSD", "rmsd");
863 error_code |= init_components_type<gyration>(conf, "radius of "
864 "gyration", "gyration");
865 error_code |= init_components_type<inertia>(conf, "moment of "
866 "inertia", "inertia");
867 error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");
868 error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
869 error_code |= init_components_type<alch_lambda>(conf, "alchemical coupling parameter", "alchLambda");
870 error_code |= init_components_type<gspath>(conf, "geometrical path collective variables (s)", "gspath");
871 error_code |= init_components_type<gzpath>(conf, "geometrical path collective variables (z)", "gzpath");
872 error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "linearCombination");
873 error_code |= init_components_type<gspathCV>(conf, "geometrical path collective variables (s) for other CVs", "gspathCV");
874 error_code |= init_components_type<gzpathCV>(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV");
875 error_code |= init_components_type<aspathCV>(conf, "arithmetic path collective variables (s) for other CVs", "aspathCV");
876 error_code |= init_components_type<azpathCV>(conf, "arithmetic path collective variables (s) for other CVs", "azpathCV");
877 error_code |= init_components_type<euler_phi>(conf, "euler phi angle of the optimal orientation", "eulerPhi");
878 error_code |= init_components_type<euler_psi>(conf, "euler psi angle of the optimal orientation", "eulerPsi");
879 error_code |= init_components_type<euler_theta>(conf, "euler theta angle of the optimal orientation", "eulerTheta");
880
881 error_code |= init_components_type<map_total>(conf, "total value of atomic map", "mapTotal");
882
883 if (!cvcs.size() || (error_code != COLVARS_OK)) {
884 cvm::error("Error: no valid components were provided "
885 "for this collective variable.\n",
886 INPUT_ERROR);
887 return INPUT_ERROR;
888 }
889
890 // Check for uniqueness of CVC names (esp. if user-provided)
891 for (i = 0; i < cvcs.size(); i++) {
892 for (j = i+1; j < cvcs.size(); j++) {
893 if (cvcs[i]->name == cvcs[j]->name) {
894 cvm::error("Components " + cvm::to_str(i) + " and " + cvm::to_str(j) +\
895 " cannot have the same name \"" + cvcs[i]->name+ "\".\n", INPUT_ERROR);
896 return INPUT_ERROR;
897 }
898 }
899 }
900
901 n_active_cvcs = cvcs.size();
902
903 // Store list of children cvcs for dependency checking purposes
904 for (i = 0; i < cvcs.size(); i++) {
905 add_child(cvcs[i]);
906 }
907
908 cvm::log("All components initialized.\n");
909
910 return COLVARS_OK;
911 }
912
913
do_feature_side_effects(int id)914 void colvar::do_feature_side_effects(int id)
915 {
916 switch (id) {
917 case f_cv_total_force_calc:
918 cvm::request_total_force();
919 break;
920 case f_cv_collect_gradient:
921 if (atom_ids.size() == 0) {
922 build_atom_list();
923 }
924 break;
925 }
926 }
927
928
build_atom_list(void)929 void colvar::build_atom_list(void)
930 {
931 // If atomic gradients are requested, build full list of atom ids from all cvcs
932 std::list<int> temp_id_list;
933
934 for (size_t i = 0; i < cvcs.size(); i++) {
935 for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
936 cvm::atom_group const &ag = *(cvcs[i]->atom_groups[j]);
937 for (size_t k = 0; k < ag.size(); k++) {
938 temp_id_list.push_back(ag[k].id);
939 }
940 if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) {
941 cvm::atom_group const &fg = *(ag.fitting_group);
942 for (size_t k = 0; k < fg.size(); k++) {
943 temp_id_list.push_back(fg[k].id);
944 }
945 }
946 }
947 }
948
949 temp_id_list.sort();
950 temp_id_list.unique();
951
952 std::list<int>::iterator li;
953 for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
954 atom_ids.push_back(*li);
955 }
956
957 temp_id_list.clear();
958
959 atomic_gradients.resize(atom_ids.size());
960 if (atom_ids.size()) {
961 if (cvm::debug())
962 cvm::log("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
963 } else {
964 cvm::log("Warning: colvar components communicated no atom IDs.\n");
965 }
966 }
967
968
parse_analysis(std::string const & conf)969 int colvar::parse_analysis(std::string const &conf)
970 {
971
972 // if (cvm::debug())
973 // cvm::log ("Parsing analysis flags for collective variable \""+
974 // this->name+"\".\n");
975
976 runave_length = 0;
977 bool b_runave = false;
978 if (get_keyval(conf, "runAve", b_runave) && b_runave) {
979
980 enable(f_cv_runave);
981
982 get_keyval(conf, "runAveLength", runave_length, 1000);
983 get_keyval(conf, "runAveStride", runave_stride, 1);
984
985 if ((cvm::restart_out_freq % runave_stride) != 0) {
986 cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
987 }
988
989 get_keyval(conf, "runAveOutputFile", runave_outfile, runave_outfile);
990 }
991
992 acf_length = 0;
993 bool b_acf = false;
994 if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {
995
996 enable(f_cv_corrfunc);
997
998 get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
999 if (acf_colvar_name == this->name) {
1000 cvm::log("Calculating auto-correlation function.\n");
1001 } else {
1002 cvm::log("Calculating correlation function with \""+
1003 this->name+"\".\n");
1004 }
1005
1006 std::string acf_type_str;
1007 get_keyval(conf, "corrFuncType", acf_type_str, to_lower_cppstr(std::string("velocity")));
1008 if (acf_type_str == to_lower_cppstr(std::string("coordinate"))) {
1009 acf_type = acf_coor;
1010 } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
1011 acf_type = acf_vel;
1012 enable(f_cv_fdiff_velocity);
1013 colvar *cv2 = cvm::colvar_by_name(acf_colvar_name);
1014 if (cv2 == NULL) {
1015 return cvm::error("Error: collective variable \""+acf_colvar_name+
1016 "\" is not defined at this time.\n", INPUT_ERROR);
1017 }
1018 cv2->enable(f_cv_fdiff_velocity); // Manual dependency to object of same type
1019 } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
1020 acf_type = acf_p2coor;
1021 } else {
1022 cvm::log("Unknown type of correlation function, \""+
1023 acf_type_str+"\".\n");
1024 cvm::set_error_bits(INPUT_ERROR);
1025 }
1026
1027 get_keyval(conf, "corrFuncOffset", acf_offset, 0);
1028 get_keyval(conf, "corrFuncLength", acf_length, 1000);
1029 get_keyval(conf, "corrFuncStride", acf_stride, 1);
1030
1031 if ((cvm::restart_out_freq % acf_stride) != 0) {
1032 cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
1033 }
1034
1035 get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
1036 get_keyval(conf, "corrFuncOutputFile", acf_outfile, acf_outfile);
1037 }
1038 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1039 }
1040
1041
init_dependencies()1042 int colvar::init_dependencies() {
1043 size_t i;
1044 if (features().size() == 0) {
1045 for (i = 0; i < f_cv_ntot; i++) {
1046 modify_features().push_back(new feature);
1047 }
1048
1049 init_feature(f_cv_active, "active", f_type_dynamic);
1050 // Do not require f_cvc_active in children, as some components may be disabled
1051 // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
1052 require_feature_alt(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
1053
1054 init_feature(f_cv_awake, "awake", f_type_static);
1055 require_feature_self(f_cv_awake, f_cv_active);
1056
1057 init_feature(f_cv_gradient, "gradient", f_type_dynamic);
1058 require_feature_children(f_cv_gradient, f_cvc_gradient);
1059
1060 init_feature(f_cv_collect_gradient, "collect_gradient", f_type_dynamic);
1061 require_feature_self(f_cv_collect_gradient, f_cv_gradient);
1062 require_feature_self(f_cv_collect_gradient, f_cv_scalar);
1063 // The following exlusion could be lifted by implementing the feature
1064 exclude_feature_self(f_cv_collect_gradient, f_cv_scripted);
1065 require_feature_children(f_cv_collect_gradient, f_cvc_explicit_gradient);
1066
1067 init_feature(f_cv_fdiff_velocity, "velocity_from_finite_differences", f_type_dynamic);
1068
1069 // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
1070 init_feature(f_cv_total_force, "total_force", f_type_dynamic);
1071 require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
1072
1073 // Deps for explicit total force calculation
1074 init_feature(f_cv_total_force_calc, "total_force_calculation", f_type_dynamic);
1075 require_feature_self(f_cv_total_force_calc, f_cv_scalar);
1076 require_feature_self(f_cv_total_force_calc, f_cv_linear);
1077 require_feature_children(f_cv_total_force_calc, f_cvc_inv_gradient);
1078 require_feature_self(f_cv_total_force_calc, f_cv_Jacobian);
1079
1080 init_feature(f_cv_Jacobian, "Jacobian_derivative", f_type_dynamic);
1081 require_feature_self(f_cv_Jacobian, f_cv_scalar);
1082 require_feature_self(f_cv_Jacobian, f_cv_linear);
1083 require_feature_children(f_cv_Jacobian, f_cvc_Jacobian);
1084
1085 init_feature(f_cv_hide_Jacobian, "hide_Jacobian_force", f_type_user);
1086 require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
1087 exclude_feature_self(f_cv_hide_Jacobian, f_cv_extended_Lagrangian);
1088
1089 init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user);
1090 require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
1091 require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient);
1092
1093 init_feature(f_cv_Langevin, "Langevin_dynamics", f_type_user);
1094 require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
1095
1096 init_feature(f_cv_external, "external", f_type_user);
1097 require_feature_self(f_cv_external, f_cv_single_cvc);
1098
1099 init_feature(f_cv_single_cvc, "single_component", f_type_static);
1100
1101 init_feature(f_cv_linear, "linear", f_type_static);
1102
1103 init_feature(f_cv_scalar, "scalar", f_type_static);
1104
1105 init_feature(f_cv_output_energy, "output_energy", f_type_user);
1106
1107 init_feature(f_cv_output_value, "output_value", f_type_user);
1108
1109 init_feature(f_cv_output_velocity, "output_velocity", f_type_user);
1110 require_feature_self(f_cv_output_velocity, f_cv_fdiff_velocity);
1111
1112 init_feature(f_cv_output_applied_force, "output_applied_force", f_type_user);
1113
1114 init_feature(f_cv_output_total_force, "output_total_force", f_type_user);
1115 require_feature_self(f_cv_output_total_force, f_cv_total_force);
1116
1117 init_feature(f_cv_subtract_applied_force, "subtract_applied_force_from_total_force", f_type_user);
1118 require_feature_self(f_cv_subtract_applied_force, f_cv_total_force);
1119
1120 init_feature(f_cv_lower_boundary, "lower_boundary", f_type_user);
1121 require_feature_self(f_cv_lower_boundary, f_cv_scalar);
1122
1123 init_feature(f_cv_upper_boundary, "upper_boundary", f_type_user);
1124 require_feature_self(f_cv_upper_boundary, f_cv_scalar);
1125
1126 init_feature(f_cv_hard_lower_boundary, "hard_lower_boundary", f_type_user);
1127 require_feature_self(f_cv_hard_lower_boundary, f_cv_lower_boundary);
1128
1129 init_feature(f_cv_hard_upper_boundary, "hard_upper_boundary", f_type_user);
1130 require_feature_self(f_cv_hard_upper_boundary, f_cv_upper_boundary);
1131
1132 init_feature(f_cv_reflecting_lower_boundary, "reflecting_lower_boundary", f_type_user);
1133 require_feature_self(f_cv_reflecting_lower_boundary, f_cv_lower_boundary);
1134 require_feature_self(f_cv_reflecting_lower_boundary, f_cv_extended_Lagrangian);
1135
1136 init_feature(f_cv_reflecting_upper_boundary, "reflecting_upper_boundary", f_type_user);
1137 require_feature_self(f_cv_reflecting_upper_boundary, f_cv_upper_boundary);
1138 require_feature_self(f_cv_reflecting_upper_boundary, f_cv_extended_Lagrangian);
1139
1140 init_feature(f_cv_grid, "grid", f_type_dynamic);
1141 require_feature_self(f_cv_grid, f_cv_lower_boundary);
1142 require_feature_self(f_cv_grid, f_cv_upper_boundary);
1143
1144 init_feature(f_cv_runave, "running_average", f_type_user);
1145
1146 init_feature(f_cv_corrfunc, "correlation_function", f_type_user);
1147
1148 init_feature(f_cv_scripted, "scripted", f_type_user);
1149
1150 init_feature(f_cv_custom_function, "custom_function", f_type_user);
1151 exclude_feature_self(f_cv_custom_function, f_cv_scripted);
1152
1153 init_feature(f_cv_periodic, "periodic", f_type_static);
1154 require_feature_self(f_cv_periodic, f_cv_scalar);
1155 init_feature(f_cv_scalar, "scalar", f_type_static);
1156 init_feature(f_cv_linear, "linear", f_type_static);
1157 init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
1158
1159 // because total forces are obtained from the previous time step,
1160 // we cannot (currently) have colvar values and total forces for the same timestep
1161 init_feature(f_cv_multiple_ts, "multiple_timestep", f_type_static);
1162 exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
1163
1164 // check that everything is initialized
1165 for (i = 0; i < colvardeps::f_cv_ntot; i++) {
1166 if (is_not_set(i)) {
1167 cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
1168 }
1169 }
1170 }
1171
1172 // Initialize feature_states for each instance
1173 feature_states.reserve(f_cv_ntot);
1174 for (i = 0; i < f_cv_ntot; i++) {
1175 feature_states.push_back(feature_state(true, false));
1176 // Most features are available, so we set them so
1177 // and list exceptions below
1178 }
1179
1180 feature_states[f_cv_fdiff_velocity].available =
1181 cvm::main()->proxy->simulation_running();
1182
1183 return COLVARS_OK;
1184 }
1185
1186
setup()1187 void colvar::setup()
1188 {
1189 // loop over all components to update masses and charges of all groups
1190 for (size_t i = 0; i < cvcs.size(); i++) {
1191 for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
1192 cvm::atom_group *atoms = cvcs[i]->atom_groups[ig];
1193 atoms->setup();
1194 atoms->print_properties(name, i, ig);
1195 atoms->read_positions();
1196 }
1197 }
1198 }
1199
1200
get_atom_lists()1201 std::vector<std::vector<int> > colvar::get_atom_lists()
1202 {
1203 std::vector<std::vector<int> > lists;
1204 for (size_t i = 0; i < cvcs.size(); i++) {
1205 std::vector<std::vector<int> > li = cvcs[i]->get_atom_lists();
1206 lists.insert(lists.end(), li.begin(), li.end());
1207 }
1208 return lists;
1209 }
1210
1211
get_volmap_ids()1212 std::vector<int> const &colvar::get_volmap_ids()
1213 {
1214 volmap_ids_.resize(cvcs.size());
1215 for (size_t i = 0; i < cvcs.size(); i++) {
1216 if (cvcs[i]->param_exists("mapID") == COLVARS_OK) {
1217 volmap_ids_[i] =
1218 *(reinterpret_cast<int const *>(cvcs[i]->get_param_ptr("mapID")));
1219 } else {
1220 volmap_ids_[i] = -1;
1221 }
1222 }
1223 return volmap_ids_;
1224 }
1225
1226
~colvar()1227 colvar::~colvar()
1228 {
1229 // There is no need to call free_children_deps() here
1230 // because the children are cvcs and will be deleted
1231 // just below
1232
1233 // Clear references to this colvar's cvcs as children
1234 // for dependency purposes
1235 remove_all_children();
1236
1237 for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();
1238 ci != cvcs.rend();
1239 ++ci) {
1240 // clear all children of this cvc (i.e. its atom groups)
1241 // because the cvc base class destructor can't do it early enough
1242 // and we don't want to have each cvc derived class do it separately
1243 (*ci)->remove_all_children();
1244 delete *ci;
1245 }
1246 cvcs.clear();
1247
1248 while (biases.size() > 0) {
1249 size_t const i = biases.size()-1;
1250 cvm::log("Warning: before deleting colvar " + name
1251 + ", deleting related bias " + biases[i]->name);
1252 delete biases[i];
1253 }
1254 biases.clear();
1255
1256 // remove reference to this colvar from the module
1257 colvarmodule *cv = cvm::main();
1258 for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
1259 cvi != cv->variables()->end();
1260 ++cvi) {
1261 if ( *cvi == this) {
1262 cv->variables()->erase(cvi);
1263 break;
1264 }
1265 }
1266
1267 cv->config_changed();
1268
1269 #ifdef LEPTON
1270 for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
1271 cei != value_evaluators.end();
1272 ++cei) {
1273 if (*cei != NULL) delete (*cei);
1274 }
1275 value_evaluators.clear();
1276
1277 for (std::vector<Lepton::CompiledExpression *>::iterator gei = gradient_evaluators.begin();
1278 gei != gradient_evaluators.end();
1279 ++gei) {
1280 if (*gei != NULL) delete (*gei);
1281 }
1282 gradient_evaluators.clear();
1283 #endif
1284 }
1285
1286
1287
1288 // ******************** CALC FUNCTIONS ********************
1289
1290
1291 // Default schedule (everything is serialized)
calc()1292 int colvar::calc()
1293 {
1294 // Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
1295 int error_code = COLVARS_OK;
1296 if (is_enabled(f_cv_active)) {
1297 error_code |= update_cvc_flags();
1298 if (error_code != COLVARS_OK) return error_code;
1299 error_code |= calc_cvcs();
1300 if (error_code != COLVARS_OK) return error_code;
1301 error_code |= collect_cvc_data();
1302 }
1303 return error_code;
1304 }
1305
1306
calc_cvcs(int first_cvc,size_t num_cvcs)1307 int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
1308 {
1309 if (cvm::debug())
1310 cvm::log("Calculating colvar \""+this->name+"\", components "+
1311 cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
1312
1313 colvarproxy *proxy = cvm::main()->proxy;
1314 int error_code = COLVARS_OK;
1315
1316 error_code |= check_cvc_range(first_cvc, num_cvcs);
1317 if (error_code != COLVARS_OK) {
1318 return error_code;
1319 }
1320
1321 if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
1322 // Use Jacobian derivative from previous timestep
1323 error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
1324 }
1325 // atom coordinates are updated by the next line
1326 error_code |= calc_cvc_values(first_cvc, num_cvcs);
1327 error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
1328 error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
1329 if (proxy->total_forces_same_step()){
1330 // Use Jacobian derivative from this timestep
1331 error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
1332 }
1333
1334 if (cvm::debug())
1335 cvm::log("Done calculating colvar \""+this->name+"\".\n");
1336
1337 return error_code;
1338 }
1339
1340
collect_cvc_data()1341 int colvar::collect_cvc_data()
1342 {
1343 if (cvm::debug())
1344 cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
1345
1346 colvarproxy *proxy = cvm::main()->proxy;
1347 int error_code = COLVARS_OK;
1348
1349 if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
1350 // Total force depends on Jacobian derivative from previous timestep
1351 // collect_cvc_total_forces() uses the previous value of jd
1352 error_code |= collect_cvc_total_forces();
1353 }
1354 error_code |= collect_cvc_values();
1355 error_code |= collect_cvc_gradients();
1356 error_code |= collect_cvc_Jacobians();
1357 if (proxy->total_forces_same_step()){
1358 // Use Jacobian derivative from this timestep
1359 error_code |= collect_cvc_total_forces();
1360 }
1361 error_code |= calc_colvar_properties();
1362
1363 if (cvm::debug())
1364 cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
1365
1366 return error_code;
1367 }
1368
1369
check_cvc_range(int first_cvc,size_t)1370 int colvar::check_cvc_range(int first_cvc, size_t /* num_cvcs */)
1371 {
1372 if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) {
1373 cvm::error("Error: trying to address a component outside the "
1374 "range defined for colvar \""+name+"\".\n", BUG_ERROR);
1375 return BUG_ERROR;
1376 }
1377 return COLVARS_OK;
1378 }
1379
1380
calc_cvc_values(int first_cvc,size_t num_cvcs)1381 int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)
1382 {
1383 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1384 size_t i, cvc_count;
1385
1386 // calculate the value of the colvar
1387
1388 if (cvm::debug())
1389 cvm::log("Calculating colvar components.\n");
1390
1391 // First, calculate component values
1392 cvm::increase_depth();
1393 for (i = first_cvc, cvc_count = 0;
1394 (i < cvcs.size()) && (cvc_count < cvc_max_count);
1395 i++) {
1396 if (!cvcs[i]->is_enabled()) continue;
1397 cvc_count++;
1398 (cvcs[i])->read_data();
1399 (cvcs[i])->calc_value();
1400 if (cvm::debug())
1401 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1402 " within colvar \""+this->name+"\" has value "+
1403 cvm::to_str((cvcs[i])->value(),
1404 cvm::cv_width, cvm::cv_prec)+".\n");
1405 }
1406 cvm::decrease_depth();
1407
1408 return COLVARS_OK;
1409 }
1410
1411
collect_cvc_values()1412 int colvar::collect_cvc_values()
1413 {
1414 x.reset();
1415
1416 // combine them appropriately, using either a scripted function or a polynomial
1417 if (is_enabled(f_cv_scripted)) {
1418 // cvcs combined by user script
1419 int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);
1420 if (res == COLVARS_NOT_IMPLEMENTED) {
1421 cvm::error("Scripted colvars are not implemented.");
1422 return COLVARS_NOT_IMPLEMENTED;
1423 }
1424 if (res != COLVARS_OK) {
1425 cvm::error("Error running scripted colvar");
1426 return COLVARS_OK;
1427 }
1428
1429 #ifdef LEPTON
1430 } else if (is_enabled(f_cv_custom_function)) {
1431
1432 size_t l = 0; // index in the vector of variable references
1433
1434 for (size_t i = 0; i < x.size(); i++) {
1435 // Fill Lepton evaluator variables with CVC values, serialized into scalars
1436 for (size_t j = 0; j < cvcs.size(); j++) {
1437 for (size_t k = 0; k < cvcs[j]->value().size(); k++) {
1438 *(value_eval_var_refs[l++]) = cvcs[j]->value()[k];
1439 }
1440 }
1441 x[i] = value_evaluators[i]->evaluate();
1442 }
1443 #endif
1444
1445 } else if (x.type() == colvarvalue::type_scalar) {
1446 // polynomial combination allowed
1447 for (size_t i = 0; i < cvcs.size(); i++) {
1448 if (!cvcs[i]->is_enabled()) continue;
1449 x += (cvcs[i])->sup_coeff *
1450 ( ((cvcs[i])->sup_np != 1) ?
1451 cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
1452 (cvcs[i])->value().real_value );
1453 }
1454 } else {
1455 for (size_t i = 0; i < cvcs.size(); i++) {
1456 if (!cvcs[i]->is_enabled()) continue;
1457 x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
1458 }
1459 }
1460
1461 if (cvm::debug())
1462 cvm::log("Colvar \""+this->name+"\" has value "+
1463 cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
1464
1465 if (after_restart) {
1466 if (cvm::proxy->simulation_running()) {
1467 cvm::real const jump2 = dist2(x, x_restart) / (width*width);
1468 if (jump2 > 0.25) {
1469 cvm::error("Error: the calculated value of colvar \""+name+
1470 "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
1471 "last read from the state file:\n"+cvm::to_str(x_restart)+
1472 "\nPossible causes are changes in configuration, "
1473 "wrong state file, or how PBC wrapping is handled.\n",
1474 INPUT_ERROR);
1475 return INPUT_ERROR;
1476 }
1477 }
1478 }
1479
1480 return COLVARS_OK;
1481 }
1482
1483
calc_cvc_gradients(int first_cvc,size_t num_cvcs)1484 int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
1485 {
1486 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1487 size_t i, cvc_count;
1488
1489 if (cvm::debug())
1490 cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
1491
1492 // calculate the gradients of each component
1493 cvm::increase_depth();
1494 for (i = first_cvc, cvc_count = 0;
1495 (i < cvcs.size()) && (cvc_count < cvc_max_count);
1496 i++) {
1497 if (!cvcs[i]->is_enabled()) continue;
1498 cvc_count++;
1499
1500 if ((cvcs[i])->is_enabled(f_cvc_gradient)) {
1501 (cvcs[i])->calc_gradients();
1502 // if requested, propagate (via chain rule) the gradients above
1503 // to the atoms used to define the roto-translation
1504 (cvcs[i])->calc_fit_gradients();
1505 if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
1506 (cvcs[i])->debug_gradients();
1507 }
1508
1509 cvm::decrease_depth();
1510
1511 if (cvm::debug())
1512 cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
1513 }
1514
1515 return COLVARS_OK;
1516 }
1517
1518
collect_cvc_gradients()1519 int colvar::collect_cvc_gradients()
1520 {
1521 size_t i;
1522 if (is_enabled(f_cv_collect_gradient)) {
1523 // Collect the atomic gradients inside colvar object
1524 for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
1525 atomic_gradients[a].reset();
1526 }
1527 for (i = 0; i < cvcs.size(); i++) {
1528 if (!cvcs[i]->is_enabled()) continue;
1529 cvcs[i]->collect_gradients(atom_ids, atomic_gradients);
1530 }
1531 }
1532 return COLVARS_OK;
1533 }
1534
1535
calc_cvc_total_force(int first_cvc,size_t num_cvcs)1536 int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
1537 {
1538 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1539 size_t i, cvc_count;
1540
1541 if (is_enabled(f_cv_total_force_calc)) {
1542 if (cvm::debug())
1543 cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
1544
1545 cvm::increase_depth();
1546
1547 for (i = first_cvc, cvc_count = 0;
1548 (i < cvcs.size()) && (cvc_count < cvc_max_count);
1549 i++) {
1550 if (!cvcs[i]->is_enabled()) continue;
1551 cvc_count++;
1552 (cvcs[i])->calc_force_invgrads();
1553 }
1554 cvm::decrease_depth();
1555
1556
1557 if (cvm::debug())
1558 cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
1559 }
1560
1561 return COLVARS_OK;
1562 }
1563
1564
collect_cvc_total_forces()1565 int colvar::collect_cvc_total_forces()
1566 {
1567 if (is_enabled(f_cv_total_force_calc)) {
1568 ft.reset();
1569
1570 if (cvm::step_relative() > 0) {
1571 // get from the cvcs the total forces from the PREVIOUS step
1572 for (size_t i = 0; i < cvcs.size(); i++) {
1573 if (!cvcs[i]->is_enabled()) continue;
1574 if (cvm::debug())
1575 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1576 " within colvar \""+this->name+"\" has total force "+
1577 cvm::to_str((cvcs[i])->total_force(),
1578 cvm::cv_width, cvm::cv_prec)+".\n");
1579 // linear combination is assumed
1580 ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
1581 }
1582 }
1583
1584 if (!(is_enabled(f_cv_hide_Jacobian) && is_enabled(f_cv_subtract_applied_force))) {
1585 // add the Jacobian force to the total force, and don't apply any silent
1586 // correction internally: biases such as colvarbias_abf will handle it
1587 // If f_cv_hide_Jacobian is enabled, a force of -fj is present in ft due to the
1588 // Jacobian-compensating force
1589 ft += fj;
1590 }
1591 }
1592
1593 return COLVARS_OK;
1594 }
1595
1596
calc_cvc_Jacobians(int first_cvc,size_t num_cvcs)1597 int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs)
1598 {
1599 size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1600
1601 if (is_enabled(f_cv_Jacobian)) {
1602 cvm::increase_depth();
1603 size_t i, cvc_count;
1604 for (i = first_cvc, cvc_count = 0;
1605 (i < cvcs.size()) && (cvc_count < cvc_max_count);
1606 i++) {
1607 if (!cvcs[i]->is_enabled()) continue;
1608 cvc_count++;
1609 (cvcs[i])->calc_Jacobian_derivative();
1610 }
1611 cvm::decrease_depth();
1612 }
1613
1614 return COLVARS_OK;
1615 }
1616
1617
collect_cvc_Jacobians()1618 int colvar::collect_cvc_Jacobians()
1619 {
1620 if (is_enabled(f_cv_Jacobian)) {
1621 fj.reset();
1622 for (size_t i = 0; i < cvcs.size(); i++) {
1623 if (!cvcs[i]->is_enabled()) continue;
1624 if (cvm::debug())
1625 cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1626 " within colvar \""+this->name+"\" has Jacobian derivative"+
1627 cvm::to_str((cvcs[i])->Jacobian_derivative(),
1628 cvm::cv_width, cvm::cv_prec)+".\n");
1629 // linear combination is assumed
1630 fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
1631 }
1632 fj *= cvm::boltzmann() * cvm::temperature();
1633 }
1634
1635 return COLVARS_OK;
1636 }
1637
1638
calc_colvar_properties()1639 int colvar::calc_colvar_properties()
1640 {
1641 if (is_enabled(f_cv_fdiff_velocity)) {
1642 // calculate the velocity by finite differences
1643 if (cvm::step_relative() == 0) {
1644 x_old = x;
1645 v_fdiff.reset(); // Do not pretend we know anything about the actual velocity
1646 // eg. upon restarting. That would require saving v_fdiff or x_old to the state file
1647 } else {
1648 v_fdiff = fdiff_velocity(x_old, x);
1649 v_reported = v_fdiff;
1650 }
1651 }
1652
1653 if (is_enabled(f_cv_extended_Lagrangian)) {
1654
1655 // initialize the restraint center in the first step to the value
1656 // just calculated from the cvcs
1657 if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset) {
1658 x_ext = x;
1659 if (is_enabled(f_cv_reflecting_lower_boundary) && x_ext < lower_boundary) {
1660 cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.");
1661 x_ext = lower_boundary;
1662 }
1663 if (is_enabled(f_cv_reflecting_upper_boundary) && x_ext > upper_boundary) {
1664 cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.");
1665 x_ext = upper_boundary;
1666 }
1667
1668 v_ext.reset(); // (already 0; added for clarity)
1669 }
1670
1671 // Special case of a repeated timestep (eg. multiple NAMD "run" statements)
1672 // revert values of the extended coordinate and velocity prior to latest integration
1673 if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) {
1674 x_ext = prev_x_ext;
1675 v_ext = prev_v_ext;
1676 }
1677 // report the restraint center as "value"
1678 // These position and velocities come from integration at the _previous timestep_ in update_forces_energy()
1679 // But we report values at the beginning of the timestep (value at t=0 on the first timestep)
1680 x_reported = x_ext;
1681 v_reported = v_ext;
1682 // the "total force" with the extended Lagrangian is
1683 // calculated in update_forces_energy() below
1684
1685 } else {
1686
1687 if (is_enabled(f_cv_subtract_applied_force)) {
1688 // correct the total force only if it has been measured
1689 // TODO add a specific test instead of relying on sq norm
1690 if (ft.norm2() > 0.0) {
1691 ft -= f_old;
1692 }
1693 }
1694
1695 x_reported = x;
1696 ft_reported = ft;
1697 }
1698
1699 // At the end of the first update after a restart, we can reset the flag
1700 after_restart = false;
1701 return COLVARS_OK;
1702 }
1703
1704
update_forces_energy()1705 cvm::real colvar::update_forces_energy()
1706 {
1707 if (cvm::debug())
1708 cvm::log("Updating colvar \""+this->name+"\".\n");
1709
1710 // set to zero the applied force
1711 f.type(value());
1712 f.reset();
1713 fr.reset();
1714
1715 // If we are not active at this timestep, that's all we have to do
1716 // return with energy == zero
1717 if (!is_enabled(f_cv_active)) return 0.;
1718
1719 // add the biases' force, which at this point should already have
1720 // been summed over each bias using this colvar
1721 // fb is already multiplied by the relevant time step factor for each bias
1722 f += fb;
1723
1724 if (is_enabled(f_cv_Jacobian)) {
1725 // the instantaneous Jacobian force was not included in the reported total force;
1726 // instead, it is subtracted from the applied force (silent Jacobian correction)
1727 // This requires the Jacobian term for the *current* timestep
1728 // Need to scale it for impulse MTS
1729 if (is_enabled(f_cv_hide_Jacobian))
1730 f -= fj * cvm::real(time_step_factor);
1731 }
1732
1733 // At this point f is the force f from external biases that will be applied to the
1734 // extended variable if there is one
1735
1736 if (is_enabled(f_cv_extended_Lagrangian)) {
1737 if (cvm::proxy->simulation_running()) {
1738 // Only integrate the extended equations of motion in running MD simulations
1739 if (cvm::debug()) {
1740 cvm::log("Updating extended-Lagrangian degree of freedom.\n");
1741 }
1742
1743 if (prev_timestep > -1L) {
1744 // Keep track of slow timestep to integrate MTS colvars
1745 // the colvar checks the interval after waking up twice
1746 cvm::step_number n_timesteps = cvm::step_relative() - prev_timestep;
1747 if (n_timesteps != 0 && n_timesteps != time_step_factor) {
1748 cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
1749 cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
1750 " steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " +
1751 cvm::to_str(cvm::step_relative()) + ").\n" +
1752 "Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n");
1753 return 0.;
1754 }
1755 }
1756
1757 // Integrate with slow timestep (if time_step_factor != 1)
1758 cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
1759
1760 colvarvalue f_ext(fr.type()); // force acting on the extended variable
1761 f_ext.reset();
1762
1763 if (is_enabled(f_cv_external)) {
1764 // There are no forces on the "actual colvar" bc there is no gradient wrt atomic coordinates
1765 // So we apply this to the extended DOF
1766 f += fb_actual;
1767 }
1768
1769 fr = f;
1770 // External force has been scaled for a 1-timestep impulse, scale it back because we will
1771 // integrate it with the colvar's own timestep factor
1772 f_ext = f / cvm::real(time_step_factor);
1773
1774 colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
1775
1776 if (is_enabled(f_cv_external)) {
1777 // Add "alchemical" force from external variable
1778 f_system = cvcs[0]->total_force();
1779 // f is now irrelevant because we are not applying atomic forces in the simulation
1780 // just driving the external variable lambda
1781 } else {
1782 // the total force is applied to the fictitious mass, while the
1783 // atoms only feel the harmonic force + wall force
1784 // fr: bias force on extended variable (without harmonic spring), for output in trajectory
1785 // f_ext: total force on extended variable (including harmonic spring)
1786 // f: - initially, external biasing force
1787 // - after this code block, colvar force to be applied to atomic coordinates
1788 // ie. spring force (fb_actual will be added just below)
1789 f_system = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
1790 f = -1.0 * f_system;
1791 // Coupling force is a slow force, to be applied to atomic coords impulse-style
1792 // over a single MD timestep
1793 f *= cvm::real(time_step_factor);
1794 }
1795 f_ext += f_system;
1796
1797 if (is_enabled(f_cv_subtract_applied_force)) {
1798 // Report a "system" force without the biases on this colvar
1799 // that is, just the spring force (or alchemical force)
1800 ft_reported = f_system;
1801 } else {
1802 // The total force acting on the extended variable is f_ext
1803 // This will be used in the next timestep
1804 ft_reported = f_ext;
1805 }
1806
1807 // backup in case we need to revert this integration timestep
1808 // if the same MD timestep is re-run
1809 prev_x_ext = x_ext;
1810 prev_v_ext = v_ext;
1811
1812 // leapfrog: starting from x_i, f_i, v_(i-1/2)
1813 v_ext += (0.5 * dt) * f_ext / ext_mass;
1814 // Because of leapfrog, kinetic energy at time i is approximate
1815 kinetic_energy = 0.5 * ext_mass * v_ext * v_ext;
1816 potential_energy = 0.5 * ext_force_k * this->dist2(x_ext, x);
1817 // leap to v_(i+1/2)
1818 if (is_enabled(f_cv_Langevin)) {
1819 v_ext -= dt * ext_gamma * v_ext;
1820 colvarvalue rnd(x);
1821 rnd.set_random();
1822 v_ext += dt * ext_sigma * rnd / ext_mass;
1823 }
1824 v_ext += (0.5 * dt) * f_ext / ext_mass;
1825 x_ext += dt * v_ext;
1826
1827 cvm::real delta = 0; // Length of overshoot past either reflecting boundary
1828 if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
1829 (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
1830 x_ext -= 2.0 * delta;
1831 v_ext *= -1.0;
1832 if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
1833 (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
1834 cvm::error("Error: extended coordinate value " + cvm::to_str(x_ext) + " is still outside boundaries after reflection.\n");
1835 }
1836 }
1837
1838 x_ext.apply_constraints();
1839 this->wrap(x_ext);
1840 if (is_enabled(f_cv_external)) {
1841 // Colvar value is constrained to the extended value
1842 x = x_ext;
1843 cvcs[0]->set_value(x_ext);
1844 }
1845 } else {
1846 // If this is a postprocessing run (eg. in VMD), the extended DOF
1847 // is equal to the actual coordinate
1848 x_ext = x;
1849 }
1850 }
1851
1852 if (!is_enabled(f_cv_external)) {
1853 // Now adding the force on the actual colvar (for those biases that
1854 // bypass the extended Lagrangian mass)
1855 f += fb_actual;
1856 }
1857
1858 if (cvm::debug())
1859 cvm::log("Done updating colvar \""+this->name+"\".\n");
1860 return (potential_energy + kinetic_energy);
1861 }
1862
1863
end_of_step()1864 int colvar::end_of_step()
1865 {
1866 if (cvm::debug())
1867 cvm::log("End of step for colvar \""+this->name+"\".\n");
1868
1869 if (is_enabled(f_cv_fdiff_velocity)) {
1870 x_old = x;
1871 }
1872
1873 if (is_enabled(f_cv_subtract_applied_force)) {
1874 f_old = f;
1875 }
1876
1877 prev_timestep = cvm::step_relative();
1878
1879 return COLVARS_OK;
1880 }
1881
1882
communicate_forces()1883 void colvar::communicate_forces()
1884 {
1885 size_t i;
1886 if (cvm::debug()) {
1887 cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
1888 cvm::log("Force to be applied: " + cvm::to_str(f) + "\n");
1889 }
1890
1891 if (is_enabled(f_cv_scripted)) {
1892 std::vector<cvm::matrix2d<cvm::real> > func_grads;
1893 func_grads.reserve(cvcs.size());
1894 for (i = 0; i < cvcs.size(); i++) {
1895 if (!cvcs[i]->is_enabled()) continue;
1896 func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
1897 cvcs[i]->value().size()));
1898 }
1899 int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
1900
1901 if (res != COLVARS_OK) {
1902 if (res == COLVARS_NOT_IMPLEMENTED) {
1903 cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);
1904 } else {
1905 cvm::error("Error running colvar gradient script");
1906 }
1907 return;
1908 }
1909
1910 int grad_index = 0; // index in the scripted gradients, to account for some components being disabled
1911 for (i = 0; i < cvcs.size(); i++) {
1912 if (!cvcs[i]->is_enabled()) continue;
1913 // cvc force is colvar force times colvar/cvc Jacobian
1914 // (vector-matrix product)
1915 (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
1916 cvcs[i]->value().type()));
1917 }
1918
1919 #ifdef LEPTON
1920 } else if (is_enabled(f_cv_custom_function)) {
1921
1922 size_t r = 0; // index in the vector of variable references
1923 size_t e = 0; // index of the gradient evaluator
1924
1925 for (size_t i = 0; i < cvcs.size(); i++) { // gradient with respect to cvc i
1926 cvm::matrix2d<cvm::real> jacobian (x.size(), cvcs[i]->value().size());
1927 for (size_t j = 0; j < cvcs[i]->value().size(); j++) { // j-th element
1928 for (size_t c = 0; c < x.size(); c++) { // derivative of scalar element c of the colvarvalue
1929
1930 // Feed cvc values to the evaluator
1931 for (size_t k = 0; k < cvcs.size(); k++) { //
1932 for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
1933 *(grad_eval_var_refs[r++]) = cvcs[k]->value()[l];
1934 }
1935 }
1936 jacobian[c][j] = gradient_evaluators[e++]->evaluate();
1937 }
1938 }
1939 // cvc force is colvar force times colvar/cvc Jacobian
1940 // (vector-matrix product)
1941 (cvcs[i])->apply_force(colvarvalue(f.as_vector() * jacobian,
1942 cvcs[i]->value().type()));
1943 }
1944 #endif
1945
1946 } else if (x.type() == colvarvalue::type_scalar) {
1947
1948 for (i = 0; i < cvcs.size(); i++) {
1949 if (!cvcs[i]->is_enabled()) continue;
1950 (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
1951 cvm::real((cvcs[i])->sup_np) *
1952 (cvm::integer_power((cvcs[i])->value().real_value,
1953 (cvcs[i])->sup_np-1)) );
1954 }
1955
1956 } else {
1957
1958 for (i = 0; i < cvcs.size(); i++) {
1959 if (!cvcs[i]->is_enabled()) continue;
1960 (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
1961 }
1962 }
1963
1964 if (cvm::debug())
1965 cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
1966 }
1967
1968
set_cvc_flags(std::vector<bool> const & flags)1969 int colvar::set_cvc_flags(std::vector<bool> const &flags)
1970 {
1971 if (flags.size() != cvcs.size()) {
1972 cvm::error("ERROR: Wrong number of CVC flags provided.");
1973 return COLVARS_ERROR;
1974 }
1975 // We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
1976 // so we store the flags that will be enforced at the next call to calc()
1977 cvc_flags = flags;
1978 return COLVARS_OK;
1979 }
1980
1981
update_active_cvc_square_norm()1982 void colvar::update_active_cvc_square_norm()
1983 {
1984 active_cvc_square_norm = 0.0;
1985 for (size_t i = 0; i < cvcs.size(); i++) {
1986 if (cvcs[i]->is_enabled()) {
1987 active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
1988 }
1989 }
1990 }
1991
1992
update_cvc_flags()1993 int colvar::update_cvc_flags()
1994 {
1995 // Update the enabled/disabled status of cvcs if necessary
1996 if (cvc_flags.size()) {
1997 n_active_cvcs = 0;
1998 for (size_t i = 0; i < cvcs.size(); i++) {
1999 cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
2000 if (cvcs[i]->is_enabled()) {
2001 n_active_cvcs++;
2002 }
2003 }
2004 if (!n_active_cvcs) {
2005 cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
2006 return COLVARS_ERROR;
2007 }
2008 cvc_flags.clear();
2009
2010 update_active_cvc_square_norm();
2011 }
2012
2013 return COLVARS_OK;
2014 }
2015
2016
update_cvc_config(std::vector<std::string> const & confs)2017 int colvar::update_cvc_config(std::vector<std::string> const &confs)
2018 {
2019 cvm::log("Updating configuration for colvar \""+name+"\"\n");
2020
2021 if (confs.size() != cvcs.size()) {
2022 return cvm::error("Error: Wrong number of CVC config strings. "
2023 "For those CVCs that are not being changed, try passing "
2024 "an empty string.", INPUT_ERROR);
2025 }
2026
2027 int error_code = COLVARS_OK;
2028 int num_changes = 0;
2029 for (size_t i = 0; i < cvcs.size(); i++) {
2030 if (confs[i].size()) {
2031 std::string conf(confs[i]);
2032 cvm::increase_depth();
2033 error_code |= cvcs[i]->colvar::cvc::init(conf);
2034 error_code |= cvcs[i]->check_keywords(conf,
2035 cvcs[i]->config_key.c_str());
2036 cvm::decrease_depth();
2037 num_changes++;
2038 }
2039 }
2040
2041 if (num_changes == 0) {
2042 cvm::log("Warning: no changes were applied through modifycvcs; "
2043 "please check that its argument is a list of strings.\n");
2044 }
2045
2046 update_active_cvc_square_norm();
2047
2048 return error_code;
2049 }
2050
2051
cvc_param_exists(std::string const & param_name)2052 int colvar::cvc_param_exists(std::string const ¶m_name)
2053 {
2054 if (is_enabled(f_cv_single_cvc)) {
2055 return cvcs[0]->param_exists(param_name);
2056 }
2057 return cvm::error("Error: calling colvar::cvc_param_exists() for a variable "
2058 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
2059 }
2060
2061
get_cvc_param(std::string const & param_name)2062 cvm::real colvar::get_cvc_param(std::string const ¶m_name)
2063 {
2064 if (is_enabled(f_cv_single_cvc)) {
2065 return cvcs[0]->get_param(param_name);
2066 }
2067 cvm::error("Error: calling colvar::get_cvc_param() for a variable "
2068 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
2069 return 0.0;
2070 }
2071
2072
get_cvc_param_ptr(std::string const & param_name)2073 void const *colvar::get_cvc_param_ptr(std::string const ¶m_name)
2074 {
2075 if (is_enabled(f_cv_single_cvc)) {
2076 return cvcs[0]->get_param_ptr(param_name);
2077 }
2078 cvm::error("Error: calling colvar::get_cvc_param() for a variable "
2079 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
2080 return NULL;
2081 }
2082
2083
get_cvc_param_grad(std::string const & param_name)2084 colvarvalue const *colvar::get_cvc_param_grad(std::string const ¶m_name)
2085 {
2086 if (is_enabled(f_cv_single_cvc)) {
2087 return cvcs[0]->get_param_grad(param_name);
2088 }
2089 cvm::error("Error: calling colvar::get_cvc_param_grad() for a variable "
2090 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
2091 return NULL;
2092 }
2093
2094
set_cvc_param(std::string const & param_name,void const * new_value)2095 int colvar::set_cvc_param(std::string const ¶m_name, void const *new_value)
2096 {
2097 if (is_enabled(f_cv_single_cvc)) {
2098 return cvcs[0]->set_param(param_name, new_value);
2099 }
2100 return cvm::error("Error: calling colvar::set_cvc_param() for a variable "
2101 "with more than one component.\n", COLVARS_NOT_IMPLEMENTED);
2102 }
2103
2104
2105 // ******************** METRIC FUNCTIONS ********************
2106 // Use the metrics defined by \link colvar::cvc \endlink objects
2107
2108
periodic_boundaries(colvarvalue const & lb,colvarvalue const & ub) const2109 bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
2110 {
2111 if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
2112 cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
2113 "requires lower and upper boundaries to be defined.\n");
2114 cvm::set_error_bits(INPUT_ERROR);
2115 }
2116
2117 if (period > 0.0) {
2118 if ( ((cvm::sqrt(this->dist2(lb, ub))) / this->width)
2119 < 1.0E-10 ) {
2120 return true;
2121 }
2122 }
2123
2124 return false;
2125 }
2126
periodic_boundaries() const2127 bool colvar::periodic_boundaries() const
2128 {
2129 if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
2130 cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
2131 "requires lower and upper boundaries to be defined.\n");
2132 }
2133
2134 return periodic_boundaries(lower_boundary, upper_boundary);
2135 }
2136
2137
dist2(colvarvalue const & x1,colvarvalue const & x2) const2138 cvm::real colvar::dist2(colvarvalue const &x1,
2139 colvarvalue const &x2) const
2140 {
2141 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
2142 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
2143 cvm::real diff = x1.real_value - x2.real_value;
2144 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
2145 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
2146 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
2147 return diff * diff;
2148 }
2149 }
2150 if (is_enabled(f_cv_homogeneous)) {
2151 return (cvcs[0])->dist2(x1, x2);
2152 } else {
2153 return x1.dist2(x2);
2154 }
2155 }
2156
dist2_lgrad(colvarvalue const & x1,colvarvalue const & x2) const2157 colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
2158 colvarvalue const &x2) const
2159 {
2160 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
2161 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
2162 cvm::real diff = x1.real_value - x2.real_value;
2163 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
2164 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
2165 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
2166 return 2.0 * diff;
2167 }
2168 }
2169 if (is_enabled(f_cv_homogeneous)) {
2170 return (cvcs[0])->dist2_lgrad(x1, x2);
2171 } else {
2172 return x1.dist2_grad(x2);
2173 }
2174 }
2175
dist2_rgrad(colvarvalue const & x1,colvarvalue const & x2) const2176 colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
2177 colvarvalue const &x2) const
2178 {
2179 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
2180 if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
2181 cvm::real diff = x1.real_value - x2.real_value;
2182 const cvm::real period_lower_boundary = wrap_center - period / 2.0;
2183 const cvm::real period_upper_boundary = wrap_center + period / 2.0;
2184 diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
2185 return (-2.0) * diff;
2186 }
2187 }
2188 if (is_enabled(f_cv_homogeneous)) {
2189 return (cvcs[0])->dist2_rgrad(x1, x2);
2190 } else {
2191 return x2.dist2_grad(x1);
2192 }
2193 }
2194
2195
wrap(colvarvalue & x_unwrapped) const2196 void colvar::wrap(colvarvalue &x_unwrapped) const
2197 {
2198 if (!is_enabled(f_cv_periodic)) {
2199 return;
2200 }
2201
2202 if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
2203 // Scripted functions do their own wrapping, as cvcs might not be periodic
2204 cvm::real shift = cvm::floor((x_unwrapped.real_value - wrap_center) /
2205 period + 0.5);
2206 x_unwrapped.real_value -= shift * period;
2207 } else {
2208 cvcs[0]->wrap(x_unwrapped);
2209 }
2210 }
2211
2212
2213 // ******************** INPUT FUNCTIONS ********************
2214
read_state(std::istream & is)2215 std::istream & colvar::read_state(std::istream &is)
2216 {
2217 std::streampos const start_pos = is.tellg();
2218
2219 std::string conf;
2220 if ( !(is >> colvarparse::read_block("colvar", &conf)) ) {
2221 // this is not a colvar block
2222 is.clear();
2223 is.seekg(start_pos, std::ios::beg);
2224 is.setstate(std::ios::failbit);
2225 return is;
2226 }
2227
2228 {
2229 std::string check_name = "";
2230 get_keyval(conf, "name", check_name,
2231 std::string(""), colvarparse::parse_silent);
2232 if (check_name.size() == 0) {
2233 cvm::error("Error: Collective variable in the "
2234 "restart file without any identifier.\n", INPUT_ERROR);
2235 is.clear();
2236 is.seekg(start_pos, std::ios::beg);
2237 is.setstate(std::ios::failbit);
2238 return is;
2239 }
2240
2241 if (check_name != name) {
2242 if (cvm::debug()) {
2243 cvm::log("Ignoring state of colvar \""+check_name+
2244 "\": this colvar is named \""+name+"\".\n");
2245 }
2246 is.seekg(start_pos, std::ios::beg);
2247 return is;
2248 }
2249 }
2250
2251 if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
2252 cvm::log("Error: restart file does not contain "
2253 "the value of the colvar \""+
2254 name+"\" .\n");
2255 } else {
2256 cvm::log("Restarting collective variable \""+name+"\" from value: "+
2257 cvm::to_str(x)+"\n");
2258 x_restart = x;
2259 after_restart = true;
2260 }
2261
2262 if (is_enabled(f_cv_extended_Lagrangian)) {
2263 if ( !(get_keyval(conf, "extended_x", x_ext,
2264 colvarvalue(x.type()), colvarparse::parse_silent)) ||
2265 !(get_keyval(conf, "extended_v", v_ext,
2266 colvarvalue(x.type()), colvarparse::parse_silent)) ) {
2267 cvm::log("Error: restart file does not contain "
2268 "\"extended_x\" or \"extended_v\" for the colvar \""+
2269 name+"\", but you requested \"extendedLagrangian\".\n");
2270 }
2271 x_reported = x_ext;
2272 } else {
2273 x_reported = x;
2274 }
2275
2276 if (is_enabled(f_cv_output_velocity)) {
2277
2278 if ( !(get_keyval(conf, "v", v_fdiff,
2279 colvarvalue(x.type()), colvarparse::parse_silent)) ) {
2280 cvm::log("Error: restart file does not contain "
2281 "the velocity for the colvar \""+
2282 name+"\", but you requested \"outputVelocity\".\n");
2283 }
2284
2285 if (is_enabled(f_cv_extended_Lagrangian)) {
2286 v_reported = v_ext;
2287 } else {
2288 v_reported = v_fdiff;
2289 }
2290 }
2291
2292 return is;
2293 }
2294
2295
read_traj(std::istream & is)2296 std::istream & colvar::read_traj(std::istream &is)
2297 {
2298 std::streampos const start_pos = is.tellg();
2299
2300 if (is_enabled(f_cv_output_value)) {
2301
2302 if (!(is >> x)) {
2303 cvm::log("Error: in reading the value of colvar \""+
2304 this->name+"\" from trajectory.\n");
2305 is.clear();
2306 is.seekg(start_pos, std::ios::beg);
2307 is.setstate(std::ios::failbit);
2308 return is;
2309 }
2310
2311 if (is_enabled(f_cv_extended_Lagrangian)) {
2312 is >> x_ext;
2313 x_reported = x_ext;
2314 } else {
2315 x_reported = x;
2316 }
2317 }
2318
2319 if (is_enabled(f_cv_output_velocity)) {
2320
2321 is >> v_fdiff;
2322
2323 if (is_enabled(f_cv_extended_Lagrangian)) {
2324 is >> v_ext;
2325 v_reported = v_ext;
2326 } else {
2327 v_reported = v_fdiff;
2328 }
2329 }
2330
2331 if (is_enabled(f_cv_output_total_force)) {
2332 is >> ft;
2333 ft_reported = ft;
2334 }
2335
2336 if (is_enabled(f_cv_output_applied_force)) {
2337 is >> f;
2338 }
2339
2340 return is;
2341 }
2342
2343
2344 // ******************** OUTPUT FUNCTIONS ********************
2345
write_state(std::ostream & os)2346 std::ostream & colvar::write_state(std::ostream &os) {
2347
2348 os << "colvar {\n"
2349 << " name " << name << "\n"
2350 << " x "
2351 << std::setprecision(cvm::cv_prec)
2352 << std::setw(cvm::cv_width)
2353 << x << "\n";
2354
2355 if (is_enabled(f_cv_output_velocity)) {
2356 os << " v "
2357 << std::setprecision(cvm::cv_prec)
2358 << std::setw(cvm::cv_width)
2359 << v_reported << "\n";
2360 }
2361
2362 if (is_enabled(f_cv_extended_Lagrangian)) {
2363 os << " extended_x "
2364 << std::setprecision(cvm::cv_prec)
2365 << std::setw(cvm::cv_width)
2366 << x_reported << "\n"
2367 << " extended_v "
2368 << std::setprecision(cvm::cv_prec)
2369 << std::setw(cvm::cv_width)
2370 << v_reported << "\n";
2371 }
2372
2373 os << "}\n\n";
2374
2375 if (runave_os) {
2376 cvm::main()->proxy->flush_output_stream(runave_os);
2377 }
2378
2379 return os;
2380 }
2381
2382
write_traj_label(std::ostream & os)2383 std::ostream & colvar::write_traj_label(std::ostream & os)
2384 {
2385 size_t const this_cv_width = x.output_width(cvm::cv_width);
2386
2387 os << " ";
2388
2389 if (is_enabled(f_cv_output_value)) {
2390
2391 os << " "
2392 << cvm::wrap_string(this->name, this_cv_width);
2393
2394 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
2395 // extended DOF
2396 os << " r_"
2397 << cvm::wrap_string(this->name, this_cv_width-2);
2398 }
2399 }
2400
2401 if (is_enabled(f_cv_output_velocity)) {
2402
2403 os << " v_"
2404 << cvm::wrap_string(this->name, this_cv_width-2);
2405
2406 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
2407 // extended DOF
2408 os << " vr_"
2409 << cvm::wrap_string(this->name, this_cv_width-3);
2410 }
2411 }
2412
2413 if (is_enabled(f_cv_output_energy)) {
2414 os << " Ep_"
2415 << cvm::wrap_string(this->name, this_cv_width-3)
2416 << " Ek_"
2417 << cvm::wrap_string(this->name, this_cv_width-3);
2418 }
2419
2420 if (is_enabled(f_cv_output_total_force)) {
2421 os << " ft_"
2422 << cvm::wrap_string(this->name, this_cv_width-3);
2423 }
2424
2425 if (is_enabled(f_cv_output_applied_force)) {
2426 os << " fa_"
2427 << cvm::wrap_string(this->name, this_cv_width-3);
2428 }
2429
2430 return os;
2431 }
2432
2433
write_traj(std::ostream & os)2434 std::ostream & colvar::write_traj(std::ostream &os)
2435 {
2436 os << " ";
2437 if (is_enabled(f_cv_output_value)) {
2438
2439 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
2440 os << " "
2441 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2442 << x;
2443 }
2444
2445 os << " "
2446 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2447 << x_reported;
2448 }
2449
2450 if (is_enabled(f_cv_output_velocity)) {
2451
2452 if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
2453 os << " "
2454 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2455 << v_fdiff;
2456 }
2457
2458 os << " "
2459 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2460 << v_reported;
2461 }
2462
2463 if (is_enabled(f_cv_output_energy)) {
2464 os << " "
2465 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2466 << potential_energy
2467 << " "
2468 << kinetic_energy;
2469 }
2470
2471
2472 if (is_enabled(f_cv_output_total_force)) {
2473 os << " "
2474 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2475 << ft_reported;
2476 }
2477
2478 if (is_enabled(f_cv_output_applied_force)) {
2479 os << " "
2480 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2481 << applied_force();
2482 }
2483
2484 return os;
2485 }
2486
2487
write_output_files()2488 int colvar::write_output_files()
2489 {
2490 int error_code = COLVARS_OK;
2491
2492 if (is_enabled(f_cv_corrfunc)) {
2493 if (acf.size()) {
2494 if (acf_outfile.size() == 0) {
2495 acf_outfile = std::string(cvm::output_prefix()+"."+this->name+
2496 ".corrfunc.dat");
2497 }
2498 cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n");
2499 cvm::backup_file(acf_outfile.c_str());
2500 std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
2501 if (!acf_os) return cvm::get_error();
2502 error_code |= write_acf(*acf_os);
2503 cvm::proxy->close_output_stream(acf_outfile);
2504 }
2505 }
2506
2507 return error_code;
2508 }
2509
2510
2511
2512 // ******************** ANALYSIS FUNCTIONS ********************
2513
analyze()2514 int colvar::analyze()
2515 {
2516 int error_code = COLVARS_OK;
2517
2518 if (is_enabled(f_cv_runave)) {
2519 error_code |= calc_runave();
2520 }
2521
2522 if (is_enabled(f_cv_corrfunc)) {
2523 error_code |= calc_acf();
2524 }
2525
2526 return error_code;
2527 }
2528
2529
history_add_value(size_t const & history_length,std::list<colvarvalue> & history,colvarvalue const & new_value)2530 inline void history_add_value(size_t const &history_length,
2531 std::list<colvarvalue> &history,
2532 colvarvalue const &new_value)
2533 {
2534 history.push_front(new_value);
2535 if (history.size() > history_length)
2536 history.pop_back();
2537 }
2538
2539
history_incr(std::list<std::list<colvarvalue>> & history,std::list<std::list<colvarvalue>>::iterator & history_p)2540 inline void history_incr(std::list< std::list<colvarvalue> > &history,
2541 std::list< std::list<colvarvalue> >::iterator &history_p)
2542 {
2543 if ((++history_p) == history.end())
2544 history_p = history.begin();
2545 }
2546
2547
calc_acf()2548 int colvar::calc_acf()
2549 {
2550 // using here an acf_stride-long list of vectors for either
2551 // coordinates (acf_x_history) or velocities (acf_v_history); each vector can
2552 // contain up to acf_length values, which are contiguous in memory
2553 // representation but separated by acf_stride in the time series;
2554 // the pointer to each vector is changed at every step
2555
2556 colvar const *cfcv = cvm::colvar_by_name(acf_colvar_name);
2557 if (cfcv == NULL) {
2558 return cvm::error("Error: collective variable \""+acf_colvar_name+
2559 "\" is not defined at this time.\n", INPUT_ERROR);
2560 }
2561
2562 if (acf_x_history.empty() && acf_v_history.empty()) {
2563
2564 // first-step operations
2565
2566 if (colvarvalue::check_types(cfcv->value(), value())) {
2567 cvm::error("Error: correlation function between \""+cfcv->name+
2568 "\" and \""+this->name+"\" cannot be calculated, "
2569 "because their value types are different.\n",
2570 INPUT_ERROR);
2571 }
2572 acf_nframes = 0;
2573
2574 cvm::log("Colvar \""+this->name+"\": initializing correlation function "
2575 "calculation.\n");
2576
2577 if (acf.size() < acf_length+1)
2578 acf.resize(acf_length+1, 0.0);
2579
2580 size_t i;
2581 switch (acf_type) {
2582
2583 case acf_vel:
2584 // allocate space for the velocities history
2585 for (i = 0; i < acf_stride; i++) {
2586 acf_v_history.push_back(std::list<colvarvalue>());
2587 }
2588 acf_v_history_p = acf_v_history.begin();
2589 break;
2590
2591 case acf_coor:
2592 case acf_p2coor:
2593 // allocate space for the coordinates history
2594 for (i = 0; i < acf_stride; i++) {
2595 acf_x_history.push_back(std::list<colvarvalue>());
2596 }
2597 acf_x_history_p = acf_x_history.begin();
2598 break;
2599
2600 case acf_notset:
2601 default:
2602 break;
2603 }
2604
2605 } else if (cvm::step_relative() > prev_timestep) {
2606
2607 switch (acf_type) {
2608
2609 case acf_vel:
2610
2611 calc_vel_acf((*acf_v_history_p), cfcv->velocity());
2612 history_add_value(acf_length+acf_offset, *acf_v_history_p,
2613 cfcv->velocity());
2614 history_incr(acf_v_history, acf_v_history_p);
2615 break;
2616
2617 case acf_coor:
2618
2619 calc_coor_acf((*acf_x_history_p), cfcv->value());
2620 history_add_value(acf_length+acf_offset, *acf_x_history_p,
2621 cfcv->value());
2622 history_incr(acf_x_history, acf_x_history_p);
2623 break;
2624
2625 case acf_p2coor:
2626
2627 calc_p2coor_acf((*acf_x_history_p), cfcv->value());
2628 history_add_value(acf_length+acf_offset, *acf_x_history_p,
2629 cfcv->value());
2630 history_incr(acf_x_history, acf_x_history_p);
2631 break;
2632
2633 case acf_notset:
2634 default:
2635 break;
2636 }
2637 }
2638
2639 return COLVARS_OK;
2640 }
2641
2642
calc_vel_acf(std::list<colvarvalue> & v_list,colvarvalue const & v)2643 void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
2644 colvarvalue const &v)
2645 {
2646 // loop over stored velocities and add to the ACF, but only if the
2647 // length is sufficient to hold an entire row of ACF values
2648 if (v_list.size() >= acf_length+acf_offset) {
2649 std::list<colvarvalue>::iterator vs_i = v_list.begin();
2650 std::vector<cvm::real>::iterator acf_i = acf.begin();
2651
2652 for (size_t i = 0; i < acf_offset; i++)
2653 ++vs_i;
2654
2655 // current vel with itself
2656 *(acf_i) += v.norm2();
2657 ++acf_i;
2658
2659 // inner products of previous velocities with current (acf_i and
2660 // vs_i are updated)
2661 colvarvalue::inner_opt(v, vs_i, v_list.end(), acf_i);
2662
2663 acf_nframes++;
2664 }
2665 }
2666
2667
calc_coor_acf(std::list<colvarvalue> & x_list,colvarvalue const & x_now)2668 void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
2669 colvarvalue const &x_now)
2670 {
2671 // same as above but for coordinates
2672 if (x_list.size() >= acf_length+acf_offset) {
2673 std::list<colvarvalue>::iterator xs_i = x_list.begin();
2674 std::vector<cvm::real>::iterator acf_i = acf.begin();
2675
2676 for (size_t i = 0; i < acf_offset; i++)
2677 ++xs_i;
2678
2679 *(acf_i++) += x.norm2();
2680
2681 colvarvalue::inner_opt(x_now, xs_i, x_list.end(), acf_i);
2682
2683 acf_nframes++;
2684 }
2685 }
2686
2687
calc_p2coor_acf(std::list<colvarvalue> & x_list,colvarvalue const & x_now)2688 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
2689 colvarvalue const &x_now)
2690 {
2691 // same as above but with second order Legendre polynomial instead
2692 // of just the scalar product
2693 if (x_list.size() >= acf_length+acf_offset) {
2694 std::list<colvarvalue>::iterator xs_i = x_list.begin();
2695 std::vector<cvm::real>::iterator acf_i = acf.begin();
2696
2697 for (size_t i = 0; i < acf_offset; i++)
2698 ++xs_i;
2699
2700 // value of P2(0) = 1
2701 *(acf_i++) += 1.0;
2702
2703 colvarvalue::p2leg_opt(x_now, xs_i, x_list.end(), acf_i);
2704
2705 acf_nframes++;
2706 }
2707 }
2708
2709
write_acf(std::ostream & os)2710 int colvar::write_acf(std::ostream &os)
2711 {
2712 if (!acf_nframes) {
2713 return COLVARS_OK;
2714 }
2715
2716 os.setf(std::ios::scientific, std::ios::floatfield);
2717 os << "# ";
2718 switch (acf_type) {
2719 case acf_vel:
2720 os << "Velocity";
2721 break;
2722 case acf_coor:
2723 os << "Coordinate";
2724 break;
2725 case acf_p2coor:
2726 os << "Coordinate (2nd Legendre poly)";
2727 break;
2728 case acf_notset:
2729 default:
2730 break;
2731 }
2732
2733 if (acf_colvar_name == name) {
2734 os << " autocorrelation function for variable \""
2735 << this->name << "\"\n";
2736 } else {
2737 os << " correlation function between variables \"" //
2738 << this->name << "\" and \"" << acf_colvar_name << "\"\n";
2739 }
2740
2741 os << "# Number of samples = ";
2742 if (acf_normalize) {
2743 os << (acf_nframes-1) << " (one DoF is used for normalization)\n";
2744 } else {
2745 os << acf_nframes << "\n";
2746 }
2747
2748 os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "
2749 << cvm::wrap_string("corrfunc(step)", cvm::cv_width) << "\n";
2750
2751 cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
2752
2753 std::vector<cvm::real>::iterator acf_i;
2754 size_t it = acf_offset;
2755 for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
2756 os << std::setw(cvm::it_width) << acf_stride * (it++) << " "
2757 << std::setprecision(cvm::cv_prec)
2758 << std::setw(cvm::cv_width)
2759 << ( acf_normalize ?
2760 (*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
2761 (*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
2762 }
2763
2764 return os.good() ? COLVARS_OK : FILE_ERROR;
2765 }
2766
2767
calc_runave()2768 int colvar::calc_runave()
2769 {
2770 int error_code = COLVARS_OK;
2771
2772 if (x_history.empty()) {
2773
2774 runave.type(value().type());
2775 runave.reset();
2776
2777 // first-step operationsf
2778
2779 if (cvm::debug())
2780 cvm::log("Colvar \""+this->name+
2781 "\": initializing running average calculation.\n");
2782
2783 acf_nframes = 0;
2784
2785 x_history.push_back(std::list<colvarvalue>());
2786 x_history_p = x_history.begin();
2787
2788 } else {
2789
2790 if ( (cvm::step_relative() % runave_stride) == 0 &&
2791 (cvm::step_relative() > prev_timestep) ) {
2792
2793 if ((*x_history_p).size() >= runave_length-1) {
2794
2795 if (runave_os == NULL) {
2796 if (runave_outfile.size() == 0) {
2797 runave_outfile = std::string(cvm::output_prefix()+"."+
2798 this->name+".runave.traj");
2799 }
2800
2801 size_t const this_cv_width = x.output_width(cvm::cv_width);
2802 cvm::proxy->backup_file(runave_outfile);
2803 runave_os = cvm::proxy->output_stream(runave_outfile);
2804 runave_os->setf(std::ios::scientific, std::ios::floatfield);
2805 *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
2806 << " "
2807 << cvm::wrap_string("running average", this_cv_width)
2808 << " "
2809 << cvm::wrap_string("running stddev", this_cv_width)
2810 << "\n";
2811 }
2812
2813 runave = x;
2814 std::list<colvarvalue>::iterator xs_i;
2815 for (xs_i = (*x_history_p).begin();
2816 xs_i != (*x_history_p).end(); ++xs_i) {
2817 runave += (*xs_i);
2818 }
2819 runave *= 1.0 / cvm::real(runave_length);
2820 runave.apply_constraints();
2821
2822 runave_variance = 0.0;
2823 runave_variance += this->dist2(x, runave);
2824 for (xs_i = (*x_history_p).begin();
2825 xs_i != (*x_history_p).end(); ++xs_i) {
2826 runave_variance += this->dist2(x, (*xs_i));
2827 }
2828 runave_variance *= 1.0 / cvm::real(runave_length-1);
2829
2830 *runave_os << std::setw(cvm::it_width) << cvm::step_relative()
2831 << " "
2832 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2833 << runave << " "
2834 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2835 << cvm::sqrt(runave_variance) << "\n";
2836 }
2837
2838 history_add_value(runave_length, *x_history_p, x);
2839 }
2840 }
2841
2842 return error_code;
2843 }
2844
2845 // Static members
2846
2847 std::vector<colvardeps::feature *> colvar::cv_features;
2848