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