Update Colvars to version 2018-12-14
[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<coordnum>(conf, "coordination "
765     "number", "coordNum");
766   error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
767     "number", "selfCoordNum");
768   error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
769     "number", "groupCoord");
770   error_code |= init_components_type<angle>(conf, "angle", "angle");
771   error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
772   error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
773   error_code |= init_components_type<h_bond>(conf, "hydrogen bond", "hBond");
774   error_code |= init_components_type<alpha_angles>(conf, "alpha helix", "alpha");
775   error_code |= init_components_type<dihedPC>(conf, "dihedral "
776     "principal component", "dihedralPC");
777   error_code |= init_components_type<orientation>(conf, "orientation", "orientation");
778   error_code |= init_components_type<orientation_angle>(conf, "orientation "
779     "angle", "orientationAngle");
780   error_code |= init_components_type<orientation_proj>(conf, "orientation "
781     "projection", "orientationProj");
782   error_code |= init_components_type<tilt>(conf, "tilt", "tilt");
783   error_code |= init_components_type<spin_angle>(conf, "spin angle", "spinAngle");
784   error_code |= init_components_type<rmsd>(conf, "RMSD", "rmsd");
785   error_code |= init_components_type<gyration>(conf, "radius of "
786     "gyration", "gyration");
787   error_code |= init_components_type<inertia>(conf, "moment of "
788     "inertia", "inertia");
789   error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");
790   error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
791
792   if (!cvcs.size() || (error_code != COLVARS_OK)) {
793     cvm::error("Error: no valid components were provided "
794                "for this collective variable.\n",
795                INPUT_ERROR);
796     return INPUT_ERROR;
797   }
798
799   n_active_cvcs = cvcs.size();
800
801   cvm::log("All components initialized.\n");
802
803   // Store list of children cvcs for dependency checking purposes
804   for (size_t i = 0; i < cvcs.size(); i++) {
805     add_child(cvcs[i]);
806   }
807
808   return COLVARS_OK;
809 }
810
811
812 void colvar::do_feature_side_effects(int id)
813 {
814   switch (id) {
815     case f_cv_total_force_calc:
816       cvm::request_total_force();
817       break;
818     case f_cv_collect_gradient:
819       if (atom_ids.size() == 0) {
820         build_atom_list();
821       }
822       break;
823   }
824 }
825
826
827 void colvar::build_atom_list(void)
828 {
829   // If atomic gradients are requested, build full list of atom ids from all cvcs
830   std::list<int> temp_id_list;
831
832   for (size_t i = 0; i < cvcs.size(); i++) {
833     for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
834       cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
835       for (size_t k = 0; k < ag.size(); k++) {
836         temp_id_list.push_back(ag[k].id);
837       }
838     }
839   }
840
841   temp_id_list.sort();
842   temp_id_list.unique();
843
844   // atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
845   unsigned int id_i = 0;
846   std::list<int>::iterator li;
847   for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
848     atom_ids[id_i] = *li;
849     id_i++;
850   }
851
852   temp_id_list.clear();
853
854   atomic_gradients.resize(atom_ids.size());
855   if (atom_ids.size()) {
856     if (cvm::debug())
857       cvm::log("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
858   } else {
859     cvm::log("Warning: colvar components communicated no atom IDs.\n");
860   }
861 }
862
863
864 int colvar::parse_analysis(std::string const &conf)
865 {
866
867   //   if (cvm::debug())
868   //     cvm::log ("Parsing analysis flags for collective variable \""+
869   //               this->name+"\".\n");
870
871   runave_length = 0;
872   bool b_runave = false;
873   if (get_keyval(conf, "runAve", b_runave) && b_runave) {
874
875     enable(f_cv_runave);
876
877     get_keyval(conf, "runAveLength", runave_length, 1000);
878     get_keyval(conf, "runAveStride", runave_stride, 1);
879
880     if ((cvm::restart_out_freq % runave_stride) != 0) {
881       cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
882     }
883
884     get_keyval(conf, "runAveOutputFile", runave_outfile, runave_outfile);
885   }
886
887   acf_length = 0;
888   bool b_acf = false;
889   if (get_keyval(conf, "corrFunc", b_acf) && b_acf) {
890
891     enable(f_cv_corrfunc);
892
893     get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
894     if (acf_colvar_name == this->name) {
895       cvm::log("Calculating auto-correlation function.\n");
896     } else {
897       cvm::log("Calculating correlation function with \""+
898                 this->name+"\".\n");
899     }
900
901     std::string acf_type_str;
902     get_keyval(conf, "corrFuncType", acf_type_str, to_lower_cppstr(std::string("velocity")));
903     if (acf_type_str == to_lower_cppstr(std::string("coordinate"))) {
904       acf_type = acf_coor;
905     } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
906       acf_type = acf_vel;
907       enable(f_cv_fdiff_velocity);
908       colvar *cv2 = cvm::colvar_by_name(acf_colvar_name);
909       if (cv2 == NULL) {
910         return cvm::error("Error: collective variable \""+acf_colvar_name+
911                           "\" is not defined at this time.\n", INPUT_ERROR);
912       }
913       cv2->enable(f_cv_fdiff_velocity);
914     } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
915       acf_type = acf_p2coor;
916     } else {
917       cvm::log("Unknown type of correlation function, \""+
918                         acf_type_str+"\".\n");
919       cvm::set_error_bits(INPUT_ERROR);
920     }
921
922     get_keyval(conf, "corrFuncOffset", acf_offset, 0);
923     get_keyval(conf, "corrFuncLength", acf_length, 1000);
924     get_keyval(conf, "corrFuncStride", acf_stride, 1);
925
926     if ((cvm::restart_out_freq % acf_stride) != 0) {
927       cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
928     }
929
930     get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
931     get_keyval(conf, "corrFuncOutputFile", acf_outfile, acf_outfile);
932   }
933   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
934 }
935
936
937 void colvar::setup() {
938   // loop over all components to reset masses of all groups
939   for (size_t i = 0; i < cvcs.size(); i++) {
940     for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
941       cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
942       atoms.setup();
943       atoms.reset_mass(name,i,ig);
944       atoms.read_positions();
945     }
946   }
947 }
948
949
950 colvar::~colvar()
951 {
952   // There is no need to call free_children_deps() here
953   // because the children are cvcs and will be deleted
954   // just below
955
956 //   Clear references to this colvar's cvcs as children
957 //   for dependency purposes
958   remove_all_children();
959
960   for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();
961       ci != cvcs.rend();
962       ++ci) {
963     // clear all children of this cvc (i.e. its atom groups)
964     // because the cvc base class destructor can't do it early enough
965     // and we don't want to have each cvc derived class do it separately
966     (*ci)->remove_all_children();
967     delete *ci;
968   }
969
970   // remove reference to this colvar from the CVM
971   colvarmodule *cv = cvm::main();
972   for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
973        cvi != cv->variables()->end();
974        ++cvi) {
975     if ( *cvi == this) {
976       cv->variables()->erase(cvi);
977       break;
978     }
979   }
980
981 #ifdef LEPTON
982   for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
983        cei != value_evaluators.end();
984        ++cei) {
985     if (*cei != NULL) delete (*cei);
986   }
987   value_evaluators.clear();
988
989   for (std::vector<Lepton::CompiledExpression *>::iterator gei = gradient_evaluators.begin();
990        gei != gradient_evaluators.end();
991        ++gei) {
992     if (*gei != NULL) delete (*gei);
993   }
994   gradient_evaluators.clear();
995 #endif
996 }
997
998
999
1000 // ******************** CALC FUNCTIONS ********************
1001
1002
1003 // Default schedule (everything is serialized)
1004 int colvar::calc()
1005 {
1006   // Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
1007   int error_code = COLVARS_OK;
1008   if (is_enabled(f_cv_active)) {
1009     error_code |= update_cvc_flags();
1010     if (error_code != COLVARS_OK) return error_code;
1011     error_code |= calc_cvcs();
1012     if (error_code != COLVARS_OK) return error_code;
1013     error_code |= collect_cvc_data();
1014   }
1015   return error_code;
1016 }
1017
1018
1019 int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
1020 {
1021   if (cvm::debug())
1022     cvm::log("Calculating colvar \""+this->name+"\", components "+
1023              cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
1024
1025   colvarproxy *proxy = cvm::main()->proxy;
1026   int error_code = COLVARS_OK;
1027
1028   error_code |= check_cvc_range(first_cvc, num_cvcs);
1029   if (error_code != COLVARS_OK) {
1030     return error_code;
1031   }
1032
1033   if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
1034     // Use Jacobian derivative from previous timestep
1035     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
1036   }
1037   // atom coordinates are updated by the next line
1038   error_code |= calc_cvc_values(first_cvc, num_cvcs);
1039   error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
1040   error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
1041   if (proxy->total_forces_same_step()){
1042     // Use Jacobian derivative from this timestep
1043     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
1044   }
1045
1046   if (cvm::debug())
1047     cvm::log("Done calculating colvar \""+this->name+"\".\n");
1048
1049   return error_code;
1050 }
1051
1052
1053 int colvar::collect_cvc_data()
1054 {
1055   if (cvm::debug())
1056     cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
1057
1058   colvarproxy *proxy = cvm::main()->proxy;
1059   int error_code = COLVARS_OK;
1060
1061   if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
1062     // Total force depends on Jacobian derivative from previous timestep
1063     // collect_cvc_total_forces() uses the previous value of jd
1064     error_code |= collect_cvc_total_forces();
1065   }
1066   error_code |= collect_cvc_values();
1067   error_code |= collect_cvc_gradients();
1068   error_code |= collect_cvc_Jacobians();
1069   if (proxy->total_forces_same_step()){
1070     // Use Jacobian derivative from this timestep
1071     error_code |= collect_cvc_total_forces();
1072   }
1073   error_code |= calc_colvar_properties();
1074
1075   if (cvm::debug())
1076     cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
1077
1078   return error_code;
1079 }
1080
1081
1082 int colvar::check_cvc_range(int first_cvc, size_t num_cvcs)
1083 {
1084   if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) {
1085     cvm::error("Error: trying to address a component outside the "
1086                "range defined for colvar \""+name+"\".\n", BUG_ERROR);
1087     return BUG_ERROR;
1088   }
1089   return COLVARS_OK;
1090 }
1091
1092
1093 int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)
1094 {
1095   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1096   size_t i, cvc_count;
1097
1098   // calculate the value of the colvar
1099
1100   if (cvm::debug())
1101     cvm::log("Calculating colvar components.\n");
1102
1103   // First, calculate component values
1104   cvm::increase_depth();
1105   for (i = first_cvc, cvc_count = 0;
1106        (i < cvcs.size()) && (cvc_count < cvc_max_count);
1107        i++) {
1108     if (!cvcs[i]->is_enabled()) continue;
1109     cvc_count++;
1110     (cvcs[i])->read_data();
1111     (cvcs[i])->calc_value();
1112     if (cvm::debug())
1113       cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1114                 " within colvar \""+this->name+"\" has value "+
1115                 cvm::to_str((cvcs[i])->value(),
1116                 cvm::cv_width, cvm::cv_prec)+".\n");
1117   }
1118   cvm::decrease_depth();
1119
1120   return COLVARS_OK;
1121 }
1122
1123
1124 int colvar::collect_cvc_values()
1125 {
1126   x.reset();
1127
1128   // combine them appropriately, using either a scripted function or a polynomial
1129   if (is_enabled(f_cv_scripted)) {
1130     // cvcs combined by user script
1131     int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);
1132     if (res == COLVARS_NOT_IMPLEMENTED) {
1133       cvm::error("Scripted colvars are not implemented.");
1134       return COLVARS_NOT_IMPLEMENTED;
1135     }
1136     if (res != COLVARS_OK) {
1137       cvm::error("Error running scripted colvar");
1138       return COLVARS_OK;
1139     }
1140
1141 #ifdef LEPTON
1142   } else if (is_enabled(f_cv_custom_function)) {
1143
1144     size_t l = 0; // index in the vector of variable references
1145
1146     for (size_t i = 0; i < x.size(); i++) {
1147       // Fill Lepton evaluator variables with CVC values, serialized into scalars
1148       for (size_t j = 0; j < cvcs.size(); j++) {
1149         for (size_t k = 0; k < cvcs[j]->value().size(); k++) {
1150           *(value_eval_var_refs[l++]) = cvcs[j]->value()[k];
1151         }
1152       }
1153       x[i] = value_evaluators[i]->evaluate();
1154     }
1155 #endif
1156
1157   } else if (x.type() == colvarvalue::type_scalar) {
1158     // polynomial combination allowed
1159     for (size_t i = 0; i < cvcs.size(); i++) {
1160       if (!cvcs[i]->is_enabled()) continue;
1161       x += (cvcs[i])->sup_coeff *
1162       ( ((cvcs[i])->sup_np != 1) ?
1163         cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
1164         (cvcs[i])->value().real_value );
1165     }
1166   } else {
1167     for (size_t i = 0; i < cvcs.size(); i++) {
1168       if (!cvcs[i]->is_enabled()) continue;
1169       x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
1170     }
1171   }
1172
1173   if (cvm::debug())
1174     cvm::log("Colvar \""+this->name+"\" has value "+
1175               cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
1176
1177   if (after_restart) {
1178     if (cvm::proxy->simulation_running()) {
1179       cvm::real const jump2 = dist2(x, x_restart) / (width*width);
1180       if (jump2 > 0.25) {
1181         cvm::error("Error: the calculated value of colvar \""+name+
1182                    "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
1183                    "last read from the state file:\n"+cvm::to_str(x_restart)+
1184                    "\nPossible causes are changes in configuration, "
1185                    "wrong state file, or how PBC wrapping is handled.\n",
1186                    INPUT_ERROR);
1187         return INPUT_ERROR;
1188       }
1189     }
1190   }
1191
1192   return COLVARS_OK;
1193 }
1194
1195
1196 int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
1197 {
1198   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1199   size_t i, cvc_count;
1200
1201   if (cvm::debug())
1202     cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
1203
1204   // calculate the gradients of each component
1205   cvm::increase_depth();
1206   for (i = first_cvc, cvc_count = 0;
1207       (i < cvcs.size()) && (cvc_count < cvc_max_count);
1208       i++) {
1209     if (!cvcs[i]->is_enabled()) continue;
1210     cvc_count++;
1211
1212     if ((cvcs[i])->is_enabled(f_cvc_gradient)) {
1213       (cvcs[i])->calc_gradients();
1214       // if requested, propagate (via chain rule) the gradients above
1215       // to the atoms used to define the roto-translation
1216      (cvcs[i])->calc_fit_gradients();
1217       if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
1218         (cvcs[i])->debug_gradients();
1219     }
1220
1221     cvm::decrease_depth();
1222
1223     if (cvm::debug())
1224       cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
1225   }
1226
1227   return COLVARS_OK;
1228 }
1229
1230
1231 int colvar::collect_cvc_gradients()
1232 {
1233   size_t i;
1234
1235   if (is_enabled(f_cv_collect_gradient)) {
1236     // Collect the atomic gradients inside colvar object
1237     for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
1238       atomic_gradients[a].reset();
1239     }
1240     for (i = 0; i < cvcs.size(); i++) {
1241       if (!cvcs[i]->is_enabled()) continue;
1242       // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
1243       cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
1244         cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
1245
1246       for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
1247
1248         cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
1249
1250         // If necessary, apply inverse rotation to get atomic
1251         // gradient in the laboratory frame
1252         if (ag.b_rotate) {
1253           cvm::rotation const rot_inv = ag.rot.inverse();
1254
1255           for (size_t k = 0; k < ag.size(); k++) {
1256             size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
1257                                         ag[k].id) - atom_ids.begin();
1258             atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
1259           }
1260
1261         } else {
1262
1263           for (size_t k = 0; k < ag.size(); k++) {
1264             size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
1265                                         ag[k].id) - atom_ids.begin();
1266             atomic_gradients[a] += coeff * ag[k].grad;
1267           }
1268         }
1269       }
1270     }
1271   }
1272   return COLVARS_OK;
1273 }
1274
1275
1276 int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
1277 {
1278   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1279   size_t i, cvc_count;
1280
1281   if (is_enabled(f_cv_total_force_calc)) {
1282     if (cvm::debug())
1283       cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
1284
1285     cvm::increase_depth();
1286
1287     for (i = first_cvc, cvc_count = 0;
1288         (i < cvcs.size()) && (cvc_count < cvc_max_count);
1289         i++) {
1290       if (!cvcs[i]->is_enabled()) continue;
1291       cvc_count++;
1292       (cvcs[i])->calc_force_invgrads();
1293     }
1294     cvm::decrease_depth();
1295
1296
1297     if (cvm::debug())
1298       cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
1299   }
1300
1301   return COLVARS_OK;
1302 }
1303
1304
1305 int colvar::collect_cvc_total_forces()
1306 {
1307   if (is_enabled(f_cv_total_force_calc)) {
1308     ft.reset();
1309
1310     if (cvm::step_relative() > 0) {
1311       // get from the cvcs the total forces from the PREVIOUS step
1312       for (size_t i = 0; i < cvcs.size();  i++) {
1313         if (!cvcs[i]->is_enabled()) continue;
1314             if (cvm::debug())
1315             cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1316                 " within colvar \""+this->name+"\" has total force "+
1317                 cvm::to_str((cvcs[i])->total_force(),
1318                 cvm::cv_width, cvm::cv_prec)+".\n");
1319         // linear combination is assumed
1320         ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
1321       }
1322     }
1323
1324     if (!is_enabled(f_cv_hide_Jacobian)) {
1325       // add the Jacobian force to the total force, and don't apply any silent
1326       // correction internally: biases such as colvarbias_abf will handle it
1327       ft += fj;
1328     }
1329   }
1330
1331   return COLVARS_OK;
1332 }
1333
1334
1335 int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs)
1336 {
1337   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
1338
1339   if (is_enabled(f_cv_Jacobian)) {
1340     cvm::increase_depth();
1341     size_t i, cvc_count;
1342     for (i = first_cvc, cvc_count = 0;
1343          (i < cvcs.size()) && (cvc_count < cvc_max_count);
1344          i++) {
1345       if (!cvcs[i]->is_enabled()) continue;
1346       cvc_count++;
1347       (cvcs[i])->calc_Jacobian_derivative();
1348     }
1349     cvm::decrease_depth();
1350   }
1351
1352   return COLVARS_OK;
1353 }
1354
1355
1356 int colvar::collect_cvc_Jacobians()
1357 {
1358   if (is_enabled(f_cv_Jacobian)) {
1359     fj.reset();
1360     for (size_t i = 0; i < cvcs.size(); i++) {
1361       if (!cvcs[i]->is_enabled()) continue;
1362         if (cvm::debug())
1363           cvm::log("Colvar component no. "+cvm::to_str(i+1)+
1364             " within colvar \""+this->name+"\" has Jacobian derivative"+
1365             cvm::to_str((cvcs[i])->Jacobian_derivative(),
1366             cvm::cv_width, cvm::cv_prec)+".\n");
1367       // linear combination is assumed
1368       fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
1369     }
1370     fj *= cvm::boltzmann() * cvm::temperature();
1371   }
1372
1373   return COLVARS_OK;
1374 }
1375
1376
1377 int colvar::calc_colvar_properties()
1378 {
1379   if (is_enabled(f_cv_fdiff_velocity)) {
1380     // calculate the velocity by finite differences
1381     if (cvm::step_relative() == 0) {
1382       x_old = x;
1383     } else {
1384       v_fdiff = fdiff_velocity(x_old, x);
1385       v_reported = v_fdiff;
1386     }
1387   }
1388
1389   if (is_enabled(f_cv_extended_Lagrangian)) {
1390
1391     // initialize the restraint center in the first step to the value
1392     // just calculated from the cvcs
1393     if (cvm::step_relative() == 0 && !after_restart) {
1394       xr = x;
1395       vr.reset(); // (already 0; added for clarity)
1396     }
1397
1398     // Special case of a repeated timestep (eg. multiple NAMD "run" statements)
1399     // revert values of the extended coordinate and velocity prior to latest integration
1400     if (cvm::step_relative() == prev_timestep) {
1401       xr = prev_xr;
1402       vr = prev_vr;
1403     }
1404
1405     // report the restraint center as "value"
1406     x_reported = xr;
1407     v_reported = vr;
1408     // the "total force" with the extended Lagrangian is
1409     // calculated in update_forces_energy() below
1410
1411   } else {
1412
1413     if (is_enabled(f_cv_subtract_applied_force)) {
1414       // correct the total force only if it has been measured
1415       // TODO add a specific test instead of relying on sq norm
1416       if (ft.norm2() > 0.0) {
1417         ft -= f_old;
1418       }
1419     }
1420
1421     x_reported = x;
1422     ft_reported = ft;
1423   }
1424
1425   // At the end of the first update after a restart, we can reset the flag
1426   after_restart = false;
1427   return COLVARS_OK;
1428 }
1429
1430
1431 cvm::real colvar::update_forces_energy()
1432 {
1433   if (cvm::debug())
1434     cvm::log("Updating colvar \""+this->name+"\".\n");
1435
1436   // set to zero the applied force
1437   f.type(value());
1438   f.reset();
1439   fr.reset();
1440
1441   // If we are not active at this timestep, that's all we have to do
1442   // return with energy == zero
1443   if (!is_enabled(f_cv_active)) return 0.;
1444
1445   // add the biases' force, which at this point should already have
1446   // been summed over each bias using this colvar
1447   f += fb;
1448
1449   if (is_enabled(f_cv_Jacobian)) {
1450     // the instantaneous Jacobian force was not included in the reported total force;
1451     // instead, it is subtracted from the applied force (silent Jacobian correction)
1452     // This requires the Jacobian term for the *current* timestep
1453     if (is_enabled(f_cv_hide_Jacobian))
1454       f -= fj;
1455   }
1456
1457   // At this point f is the force f from external biases that will be applied to the
1458   // extended variable if there is one
1459
1460   if (is_enabled(f_cv_extended_Lagrangian)) {
1461     if (cvm::debug()) {
1462       cvm::log("Updating extended-Lagrangian degree of freedom.\n");
1463     }
1464
1465     if (prev_timestep > -1) {
1466       // Keep track of slow timestep to integrate MTS colvars
1467       // the colvar checks the interval after waking up twice
1468       int n_timesteps = cvm::step_relative() - prev_timestep;
1469       if (n_timesteps != 0 && n_timesteps != time_step_factor) {
1470         cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
1471           cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
1472           " steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " +
1473           cvm::to_str(cvm::step_relative()) + ").\n" +
1474           "Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n");
1475         return 0.;
1476       }
1477     }
1478
1479     // Integrate with slow timestep (if time_step_factor != 1)
1480     cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
1481
1482     colvarvalue f_ext(fr.type()); // force acting on the extended variable
1483     f_ext.reset();
1484
1485     // the total force is applied to the fictitious mass, while the
1486     // atoms only feel the harmonic force + wall force
1487     // fr: bias force on extended variable (without harmonic spring), for output in trajectory
1488     // f_ext: total force on extended variable (including harmonic spring)
1489     // f: - initially, external biasing force
1490     //    - after this code block, colvar force to be applied to atomic coordinates
1491     //      ie. spring force (fb_actual will be added just below)
1492     fr    = f;
1493     // External force has been scaled for a 1-timestep impulse, scale it back because we will
1494     // integrate it with the colvar's own timestep factor
1495     f_ext = f / cvm::real(time_step_factor);
1496     f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
1497     f      = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
1498     // Coupling force is a slow force, to be applied to atomic coords impulse-style
1499     f *= cvm::real(time_step_factor);
1500
1501     if (is_enabled(f_cv_subtract_applied_force)) {
1502       // Report a "system" force without the biases on this colvar
1503       // that is, just the spring force
1504       ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
1505     } else {
1506       // The total force acting on the extended variable is f_ext
1507       // This will be used in the next timestep
1508       ft_reported = f_ext;
1509     }
1510
1511     // backup in case we need to revert this integration timestep
1512     // if the same MD timestep is re-run
1513     prev_xr = xr;
1514     prev_vr = vr;
1515
1516     // leapfrog: starting from x_i, f_i, v_(i-1/2)
1517     vr  += (0.5 * dt) * f_ext / ext_mass;
1518     // Because of leapfrog, kinetic energy at time i is approximate
1519     kinetic_energy = 0.5 * ext_mass * vr * vr;
1520     potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
1521     // leap to v_(i+1/2)
1522     if (is_enabled(f_cv_Langevin)) {
1523       vr -= dt * ext_gamma * vr;
1524       colvarvalue rnd(x);
1525       rnd.set_random();
1526       vr += dt * ext_sigma * rnd / ext_mass;
1527     }
1528     vr  += (0.5 * dt) * f_ext / ext_mass;
1529     xr  += dt * vr;
1530     xr.apply_constraints();
1531     this->wrap(xr);
1532   }
1533
1534   // Now adding the force on the actual colvar (for those biases that
1535   // bypass the extended Lagrangian mass)
1536   f += fb_actual;
1537
1538   if (cvm::debug())
1539     cvm::log("Done updating colvar \""+this->name+"\".\n");
1540   return (potential_energy + kinetic_energy);
1541 }
1542
1543
1544 int colvar::end_of_step()
1545 {
1546   if (cvm::debug())
1547     cvm::log("End of step for colvar \""+this->name+"\".\n");
1548
1549   if (is_enabled(f_cv_fdiff_velocity)) {
1550     x_old = x;
1551   }
1552
1553   if (is_enabled(f_cv_subtract_applied_force)) {
1554     f_old = f;
1555   }
1556
1557   prev_timestep = cvm::step_relative();
1558
1559   return COLVARS_OK;
1560 }
1561
1562
1563 void colvar::communicate_forces()
1564 {
1565   size_t i;
1566   if (cvm::debug()) {
1567     cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
1568     cvm::log("Force to be applied: " + cvm::to_str(f) + "\n");
1569   }
1570
1571   if (is_enabled(f_cv_scripted)) {
1572     std::vector<cvm::matrix2d<cvm::real> > func_grads;
1573     func_grads.reserve(cvcs.size());
1574     for (i = 0; i < cvcs.size(); i++) {
1575       if (!cvcs[i]->is_enabled()) continue;
1576       func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
1577                                                      cvcs[i]->value().size()));
1578     }
1579     int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
1580
1581     if (res != COLVARS_OK) {
1582       if (res == COLVARS_NOT_IMPLEMENTED) {
1583         cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);
1584       } else {
1585         cvm::error("Error running colvar gradient script");
1586       }
1587       return;
1588     }
1589
1590     int grad_index = 0; // index in the scripted gradients, to account for some components being disabled
1591     for (i = 0; i < cvcs.size(); i++) {
1592       if (!cvcs[i]->is_enabled()) continue;
1593       // cvc force is colvar force times colvar/cvc Jacobian
1594       // (vector-matrix product)
1595       (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
1596                              cvcs[i]->value().type()));
1597     }
1598
1599 #ifdef LEPTON
1600   } else if (is_enabled(f_cv_custom_function)) {
1601
1602     size_t r = 0; // index in the vector of variable references
1603     size_t e = 0; // index of the gradient evaluator
1604
1605     for (size_t i = 0; i < cvcs.size(); i++) {  // gradient with respect to cvc i
1606       cvm::matrix2d<cvm::real> jacobian (x.size(), cvcs[i]->value().size());
1607       for (size_t j = 0; j < cvcs[i]->value().size(); j++) { // j-th element
1608         for (size_t c = 0; c < x.size(); c++) { // derivative of scalar element c of the colvarvalue
1609
1610           // Feed cvc values to the evaluator
1611           for (size_t k = 0; k < cvcs.size(); k++) { //
1612             for (size_t l = 0; l < cvcs[k]->value().size(); l++) {
1613               *(grad_eval_var_refs[r++]) = cvcs[k]->value()[l];
1614             }
1615           }
1616           jacobian[c][j] = gradient_evaluators[e++]->evaluate();
1617         }
1618       }
1619       // cvc force is colvar force times colvar/cvc Jacobian
1620       // (vector-matrix product)
1621       (cvcs[i])->apply_force(colvarvalue(f.as_vector() * jacobian,
1622                              cvcs[i]->value().type()));
1623     }
1624 #endif
1625
1626   } else if (x.type() == colvarvalue::type_scalar) {
1627
1628     for (i = 0; i < cvcs.size(); i++) {
1629       if (!cvcs[i]->is_enabled()) continue;
1630       (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
1631                              cvm::real((cvcs[i])->sup_np) *
1632                              (cvm::integer_power((cvcs[i])->value().real_value,
1633                                                  (cvcs[i])->sup_np-1)) );
1634     }
1635
1636   } else {
1637
1638     for (i = 0; i < cvcs.size(); i++) {
1639       if (!cvcs[i]->is_enabled()) continue;
1640       (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
1641     }
1642   }
1643
1644   if (cvm::debug())
1645     cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
1646 }
1647
1648
1649 int colvar::set_cvc_flags(std::vector<bool> const &flags)
1650 {
1651   if (flags.size() != cvcs.size()) {
1652     cvm::error("ERROR: Wrong number of CVC flags provided.");
1653     return COLVARS_ERROR;
1654   }
1655   // We cannot enable or disable cvcs in the middle of a timestep or colvar evaluation sequence
1656   // so we store the flags that will be enforced at the next call to calc()
1657   cvc_flags = flags;
1658   return COLVARS_OK;
1659 }
1660
1661
1662 void colvar::update_active_cvc_square_norm()
1663 {
1664   active_cvc_square_norm = 0.0;
1665   for (size_t i = 0; i < cvcs.size(); i++) {
1666     if (cvcs[i]->is_enabled()) {
1667       active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
1668     }
1669   }
1670 }
1671
1672
1673 int colvar::update_cvc_flags()
1674 {
1675   // Update the enabled/disabled status of cvcs if necessary
1676   if (cvc_flags.size()) {
1677     n_active_cvcs = 0;
1678     for (size_t i = 0; i < cvcs.size(); i++) {
1679       cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
1680       if (cvcs[i]->is_enabled()) {
1681         n_active_cvcs++;
1682       }
1683     }
1684     if (!n_active_cvcs) {
1685       cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
1686       return COLVARS_ERROR;
1687     }
1688     cvc_flags.clear();
1689
1690     update_active_cvc_square_norm();
1691   }
1692
1693   return COLVARS_OK;
1694 }
1695
1696
1697 int colvar::update_cvc_config(std::vector<std::string> const &confs)
1698 {
1699   cvm::log("Updating configuration for colvar \""+name+"\"");
1700
1701   if (confs.size() != cvcs.size()) {
1702     return cvm::error("Error: Wrong number of CVC config strings.  "
1703                       "For those CVCs that are not being changed, try passing "
1704                       "an empty string.", INPUT_ERROR);
1705   }
1706
1707   int error_code = COLVARS_OK;
1708   int num_changes = 0;
1709   for (size_t i = 0; i < cvcs.size(); i++) {
1710     if (confs[i].size()) {
1711       std::string conf(confs[i]);
1712       cvm::increase_depth();
1713       error_code |= cvcs[i]->colvar::cvc::init(conf);
1714       error_code |= cvcs[i]->check_keywords(conf,
1715                                             cvcs[i]->config_key.c_str());
1716       cvm::decrease_depth();
1717       num_changes++;
1718     }
1719   }
1720
1721   if (num_changes == 0) {
1722     cvm::log("Warning: no changes were applied through modifycvcs; "
1723              "please check that its argument is a list of strings.\n");
1724   }
1725
1726   update_active_cvc_square_norm();
1727
1728   return error_code;
1729 }
1730
1731
1732 // ******************** METRIC FUNCTIONS ********************
1733 // Use the metrics defined by \link cvc \endlink objects
1734
1735
1736 bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
1737 {
1738   if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
1739     cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
1740                     "requires lower and upper boundaries to be defined.\n");
1741     cvm::set_error_bits(INPUT_ERROR);
1742   }
1743
1744   if (period > 0.0) {
1745     if ( ((std::sqrt(this->dist2(lb, ub))) / this->width)
1746          < 1.0E-10 ) {
1747       return true;
1748     }
1749   }
1750
1751   return false;
1752 }
1753
1754 bool colvar::periodic_boundaries() const
1755 {
1756   if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
1757     cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
1758                     "requires lower and upper boundaries to be defined.\n");
1759   }
1760
1761   return periodic_boundaries(lower_boundary, upper_boundary);
1762 }
1763
1764
1765 cvm::real colvar::dist2(colvarvalue const &x1,
1766                          colvarvalue const &x2) const
1767 {
1768   if (is_enabled(f_cv_homogeneous)) {
1769     return (cvcs[0])->dist2(x1, x2);
1770   } else {
1771     return x1.dist2(x2);
1772   }
1773 }
1774
1775 colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
1776                                  colvarvalue const &x2) const
1777 {
1778   if (is_enabled(f_cv_homogeneous)) {
1779     return (cvcs[0])->dist2_lgrad(x1, x2);
1780   } else {
1781     return x1.dist2_grad(x2);
1782   }
1783 }
1784
1785 colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
1786                                  colvarvalue const &x2) const
1787 {
1788   if (is_enabled(f_cv_homogeneous)) {
1789     return (cvcs[0])->dist2_rgrad(x1, x2);
1790   } else {
1791     return x2.dist2_grad(x1);
1792   }
1793 }
1794
1795 void colvar::wrap(colvarvalue &x) const
1796 {
1797   if ( !is_enabled(f_cv_periodic) ) {
1798     return;
1799   }
1800
1801   if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
1802     // Scripted functions do their own wrapping, as cvcs might not be periodic
1803     cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
1804     x.real_value -= shift * period;
1805   } else {
1806     cvcs[0]->wrap(x);
1807   }
1808
1809   return;
1810 }
1811
1812
1813 // ******************** INPUT FUNCTIONS ********************
1814
1815 std::istream & colvar::read_restart(std::istream &is)
1816 {
1817   size_t const start_pos = is.tellg();
1818
1819   std::string conf;
1820   if ( !(is >> colvarparse::read_block("colvar", conf)) ) {
1821     // this is not a colvar block
1822     is.clear();
1823     is.seekg(start_pos, std::ios::beg);
1824     is.setstate(std::ios::failbit);
1825     return is;
1826   }
1827
1828   {
1829     std::string check_name = "";
1830     if ( (get_keyval(conf, "name", check_name,
1831                      std::string(""), colvarparse::parse_silent)) &&
1832          (check_name != name) )  {
1833       cvm::error("Error: the state file does not match the "
1834                  "configuration file, at colvar \""+name+"\".\n");
1835     }
1836     if (check_name.size() == 0) {
1837       cvm::error("Error: Collective variable in the "
1838                  "restart file without any identifier.\n");
1839     }
1840   }
1841
1842   if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
1843     cvm::log("Error: restart file does not contain "
1844              "the value of the colvar \""+
1845              name+"\" .\n");
1846   } else {
1847     cvm::log("Restarting collective variable \""+name+"\" from value: "+
1848              cvm::to_str(x)+"\n");
1849     x_restart = x;
1850     after_restart = true;
1851   }
1852
1853   if (is_enabled(f_cv_extended_Lagrangian)) {
1854
1855     if ( !(get_keyval(conf, "extended_x", xr,
1856                       colvarvalue(x.type()), colvarparse::parse_silent)) &&
1857          !(get_keyval(conf, "extended_v", vr,
1858                       colvarvalue(x.type()), colvarparse::parse_silent)) ) {
1859       cvm::log("Error: restart file does not contain "
1860                "\"extended_x\" or \"extended_v\" for the colvar \""+
1861                name+"\", but you requested \"extendedLagrangian\".\n");
1862     }
1863     x_reported = xr;
1864   } else {
1865     x_reported = x;
1866   }
1867
1868   if (is_enabled(f_cv_output_velocity)) {
1869
1870     if ( !(get_keyval(conf, "v", v_fdiff,
1871                       colvarvalue(x.type()), colvarparse::parse_silent)) ) {
1872       cvm::log("Error: restart file does not contain "
1873                "the velocity for the colvar \""+
1874                name+"\", but you requested \"outputVelocity\".\n");
1875     }
1876
1877     if (is_enabled(f_cv_extended_Lagrangian)) {
1878       v_reported = vr;
1879     } else {
1880       v_reported = v_fdiff;
1881     }
1882   }
1883
1884   return is;
1885 }
1886
1887
1888 std::istream & colvar::read_traj(std::istream &is)
1889 {
1890   size_t const start_pos = is.tellg();
1891
1892   if (is_enabled(f_cv_output_value)) {
1893
1894     if (!(is >> x)) {
1895       cvm::log("Error: in reading the value of colvar \""+
1896                 this->name+"\" from trajectory.\n");
1897       is.clear();
1898       is.seekg(start_pos, std::ios::beg);
1899       is.setstate(std::ios::failbit);
1900       return is;
1901     }
1902
1903     if (is_enabled(f_cv_extended_Lagrangian)) {
1904       is >> xr;
1905       x_reported = xr;
1906     } else {
1907       x_reported = x;
1908     }
1909   }
1910
1911   if (is_enabled(f_cv_output_velocity)) {
1912
1913     is >> v_fdiff;
1914
1915     if (is_enabled(f_cv_extended_Lagrangian)) {
1916       is >> vr;
1917       v_reported = vr;
1918     } else {
1919       v_reported = v_fdiff;
1920     }
1921   }
1922
1923   if (is_enabled(f_cv_output_total_force)) {
1924     is >> ft;
1925     ft_reported = ft;
1926   }
1927
1928   if (is_enabled(f_cv_output_applied_force)) {
1929     is >> f;
1930   }
1931
1932   return is;
1933 }
1934
1935
1936 // ******************** OUTPUT FUNCTIONS ********************
1937
1938 std::ostream & colvar::write_restart(std::ostream &os) {
1939
1940   os << "colvar {\n"
1941      << "  name " << name << "\n"
1942      << "  x "
1943      << std::setprecision(cvm::cv_prec)
1944      << std::setw(cvm::cv_width)
1945      << x << "\n";
1946
1947   if (is_enabled(f_cv_output_velocity)) {
1948     os << "  v "
1949        << std::setprecision(cvm::cv_prec)
1950        << std::setw(cvm::cv_width)
1951        << v_reported << "\n";
1952   }
1953
1954   if (is_enabled(f_cv_extended_Lagrangian)) {
1955     os << "  extended_x "
1956        << std::setprecision(cvm::cv_prec)
1957        << std::setw(cvm::cv_width)
1958        << xr << "\n"
1959        << "  extended_v "
1960        << std::setprecision(cvm::cv_prec)
1961        << std::setw(cvm::cv_width)
1962        << vr << "\n";
1963   }
1964
1965   os << "}\n\n";
1966
1967   if (runave_os) {
1968     cvm::main()->proxy->flush_output_stream(runave_os);
1969   }
1970
1971   return os;
1972 }
1973
1974
1975 std::ostream & colvar::write_traj_label(std::ostream & os)
1976 {
1977   size_t const this_cv_width = x.output_width(cvm::cv_width);
1978
1979   os << " ";
1980
1981   if (is_enabled(f_cv_output_value)) {
1982
1983     os << " "
1984        << cvm::wrap_string(this->name, this_cv_width);
1985
1986     if (is_enabled(f_cv_extended_Lagrangian)) {
1987       // extended DOF
1988       os << " r_"
1989          << cvm::wrap_string(this->name, this_cv_width-2);
1990     }
1991   }
1992
1993   if (is_enabled(f_cv_output_velocity)) {
1994
1995     os << " v_"
1996        << cvm::wrap_string(this->name, this_cv_width-2);
1997
1998     if (is_enabled(f_cv_extended_Lagrangian)) {
1999       // extended DOF
2000       os << " vr_"
2001          << cvm::wrap_string(this->name, this_cv_width-3);
2002     }
2003   }
2004
2005   if (is_enabled(f_cv_output_energy)) {
2006     os << " Ep_"
2007        << cvm::wrap_string(this->name, this_cv_width-3)
2008        << " Ek_"
2009        << cvm::wrap_string(this->name, this_cv_width-3);
2010   }
2011
2012   if (is_enabled(f_cv_output_total_force)) {
2013     os << " ft_"
2014        << cvm::wrap_string(this->name, this_cv_width-3);
2015   }
2016
2017   if (is_enabled(f_cv_output_applied_force)) {
2018     os << " fa_"
2019        << cvm::wrap_string(this->name, this_cv_width-3);
2020   }
2021
2022   return os;
2023 }
2024
2025
2026 std::ostream & colvar::write_traj(std::ostream &os)
2027 {
2028   os << " ";
2029
2030   if (is_enabled(f_cv_output_value)) {
2031
2032     if (is_enabled(f_cv_extended_Lagrangian)) {
2033       os << " "
2034          << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2035          << x;
2036     }
2037
2038     os << " "
2039        << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2040        << x_reported;
2041   }
2042
2043   if (is_enabled(f_cv_output_velocity)) {
2044
2045     if (is_enabled(f_cv_extended_Lagrangian)) {
2046       os << " "
2047          << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2048          << v_fdiff;
2049     }
2050
2051     os << " "
2052        << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2053        << v_reported;
2054   }
2055
2056   if (is_enabled(f_cv_output_energy)) {
2057     os << " "
2058        << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2059        << potential_energy
2060        << " "
2061        << kinetic_energy;
2062   }
2063
2064
2065   if (is_enabled(f_cv_output_total_force)) {
2066     os << " "
2067        << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2068        << ft_reported;
2069   }
2070
2071   if (is_enabled(f_cv_output_applied_force)) {
2072     os << " "
2073        << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2074        << applied_force();
2075   }
2076
2077   return os;
2078 }
2079
2080
2081 int colvar::write_output_files()
2082 {
2083   int error_code = COLVARS_OK;
2084
2085   if (is_enabled(f_cv_corrfunc)) {
2086     if (acf.size()) {
2087       if (acf_outfile.size() == 0) {
2088         acf_outfile = std::string(cvm::output_prefix()+"."+this->name+
2089                                   ".corrfunc.dat");
2090       }
2091       cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n");
2092       cvm::backup_file(acf_outfile.c_str());
2093       std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
2094       if (!acf_os) return cvm::get_error();
2095       error_code |= write_acf(*acf_os);
2096       cvm::proxy->close_output_stream(acf_outfile);
2097     }
2098   }
2099
2100   return error_code;
2101 }
2102
2103
2104
2105 // ******************** ANALYSIS FUNCTIONS ********************
2106
2107 int colvar::analyze()
2108 {
2109   int error_code = COLVARS_OK;
2110
2111   if (is_enabled(f_cv_runave)) {
2112     error_code |= calc_runave();
2113   }
2114
2115   if (is_enabled(f_cv_corrfunc)) {
2116     error_code |= calc_acf();
2117   }
2118
2119   return error_code;
2120 }
2121
2122
2123 inline void history_add_value(size_t const           &history_length,
2124                               std::list<colvarvalue> &history,
2125                               colvarvalue const      &new_value)
2126 {
2127   history.push_front(new_value);
2128   if (history.size() > history_length)
2129     history.pop_back();
2130 }
2131
2132
2133 inline void history_incr(std::list< std::list<colvarvalue> >           &history,
2134                          std::list< std::list<colvarvalue> >::iterator &history_p)
2135 {
2136   if ((++history_p) == history.end())
2137     history_p = history.begin();
2138 }
2139
2140
2141 int colvar::calc_acf()
2142 {
2143   // using here an acf_stride-long list of vectors for either
2144   // coordinates (acf_x_history) or velocities (acf_v_history); each vector can
2145   // contain up to acf_length values, which are contiguous in memory
2146   // representation but separated by acf_stride in the time series;
2147   // the pointer to each vector is changed at every step
2148
2149   colvar const *cfcv = cvm::colvar_by_name(acf_colvar_name);
2150   if (cfcv == NULL) {
2151     return cvm::error("Error: collective variable \""+acf_colvar_name+
2152                       "\" is not defined at this time.\n", INPUT_ERROR);
2153   }
2154
2155   if (acf_x_history.empty() && acf_v_history.empty()) {
2156
2157     // first-step operations
2158
2159     if (colvarvalue::check_types(cfcv->value(), value())) {
2160       cvm::error("Error: correlation function between \""+cfcv->name+
2161                  "\" and \""+this->name+"\" cannot be calculated, "
2162                  "because their value types are different.\n",
2163                  INPUT_ERROR);
2164     }
2165     acf_nframes = 0;
2166
2167     cvm::log("Colvar \""+this->name+"\": initializing correlation function "
2168              "calculation.\n");
2169
2170     if (acf.size() < acf_length+1)
2171       acf.resize(acf_length+1, 0.0);
2172
2173     size_t i;
2174     switch (acf_type) {
2175
2176     case acf_vel:
2177       // allocate space for the velocities history
2178       for (i = 0; i < acf_stride; i++) {
2179         acf_v_history.push_back(std::list<colvarvalue>());
2180       }
2181       acf_v_history_p = acf_v_history.begin();
2182       break;
2183
2184     case acf_coor:
2185     case acf_p2coor:
2186       // allocate space for the coordinates history
2187       for (i = 0; i < acf_stride; i++) {
2188         acf_x_history.push_back(std::list<colvarvalue>());
2189       }
2190       acf_x_history_p = acf_x_history.begin();
2191       break;
2192
2193     default:
2194       break;
2195     }
2196
2197   } else if (cvm::step_relative() > prev_timestep) {
2198
2199     switch (acf_type) {
2200
2201     case acf_vel:
2202
2203       calc_vel_acf((*acf_v_history_p), cfcv->velocity());
2204       history_add_value(acf_length+acf_offset, *acf_v_history_p,
2205                         cfcv->velocity());
2206       history_incr(acf_v_history, acf_v_history_p);
2207       break;
2208
2209     case acf_coor:
2210
2211       calc_coor_acf((*acf_x_history_p), cfcv->value());
2212       history_add_value(acf_length+acf_offset, *acf_x_history_p,
2213                         cfcv->value());
2214       history_incr(acf_x_history, acf_x_history_p);
2215       break;
2216
2217     case acf_p2coor:
2218
2219       calc_p2coor_acf((*acf_x_history_p), cfcv->value());
2220       history_add_value(acf_length+acf_offset, *acf_x_history_p,
2221                         cfcv->value());
2222       history_incr(acf_x_history, acf_x_history_p);
2223       break;
2224
2225     default:
2226       break;
2227     }
2228   }
2229
2230   return COLVARS_OK;
2231 }
2232
2233
2234 void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
2235                           colvarvalue const      &v)
2236 {
2237   // loop over stored velocities and add to the ACF, but only if the
2238   // length is sufficient to hold an entire row of ACF values
2239   if (v_list.size() >= acf_length+acf_offset) {
2240     std::list<colvarvalue>::iterator  vs_i = v_list.begin();
2241     std::vector<cvm::real>::iterator acf_i = acf.begin();
2242
2243     for (size_t i = 0; i < acf_offset; i++)
2244       ++vs_i;
2245
2246     // current vel with itself
2247     *(acf_i) += v.norm2();
2248     ++acf_i;
2249
2250     // inner products of previous velocities with current (acf_i and
2251     // vs_i are updated)
2252     colvarvalue::inner_opt(v, vs_i, v_list.end(), acf_i);
2253
2254     acf_nframes++;
2255   }
2256 }
2257
2258
2259 void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
2260                             colvarvalue const      &x)
2261 {
2262   // same as above but for coordinates
2263   if (x_list.size() >= acf_length+acf_offset) {
2264     std::list<colvarvalue>::iterator  xs_i = x_list.begin();
2265     std::vector<cvm::real>::iterator acf_i = acf.begin();
2266
2267     for (size_t i = 0; i < acf_offset; i++)
2268       ++xs_i;
2269
2270     *(acf_i++) += x.norm2();
2271
2272     colvarvalue::inner_opt(x, xs_i, x_list.end(), acf_i);
2273
2274     acf_nframes++;
2275   }
2276 }
2277
2278
2279 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
2280                              colvarvalue const      &x)
2281 {
2282   // same as above but with second order Legendre polynomial instead
2283   // of just the scalar product
2284   if (x_list.size() >= acf_length+acf_offset) {
2285     std::list<colvarvalue>::iterator  xs_i = x_list.begin();
2286     std::vector<cvm::real>::iterator acf_i = acf.begin();
2287
2288     for (size_t i = 0; i < acf_offset; i++)
2289       ++xs_i;
2290
2291     // value of P2(0) = 1
2292     *(acf_i++) += 1.0;
2293
2294     colvarvalue::p2leg_opt(x, xs_i, x_list.end(), acf_i);
2295
2296     acf_nframes++;
2297   }
2298 }
2299
2300
2301 int colvar::write_acf(std::ostream &os)
2302 {
2303   if (!acf_nframes) {
2304     return COLVARS_OK;
2305   }
2306
2307   os.setf(std::ios::scientific, std::ios::floatfield);
2308   os << "# ";
2309   switch (acf_type) {
2310   case acf_vel:
2311     os << "Velocity";
2312     break;
2313   case acf_coor:
2314     os << "Coordinate";
2315     break;
2316   case acf_p2coor:
2317     os << "Coordinate (2nd Legendre poly)";
2318     break;
2319   }
2320
2321   if (acf_colvar_name == name) {
2322     os << " autocorrelation function for variable \""
2323        << this->name << "\"\n";
2324   } else {
2325     os << " correlation function between variables \"" //
2326        << this->name << "\" and \"" << acf_colvar_name << "\"\n";
2327   }
2328
2329   os << "# Number of samples = ";
2330   if (acf_normalize) {
2331     os << (acf_nframes-1) << " (one DoF is used for normalization)\n";
2332   } else {
2333     os << acf_nframes << "\n";
2334   }
2335
2336   os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "
2337      << cvm::wrap_string("corrfunc(step)", cvm::cv_width) << "\n";
2338
2339   cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
2340
2341   std::vector<cvm::real>::iterator acf_i;
2342   size_t it = acf_offset;
2343   for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
2344     os << std::setw(cvm::it_width) << acf_stride * (it++) << " "
2345        << std::setprecision(cvm::cv_prec)
2346        << std::setw(cvm::cv_width)
2347        << ( acf_normalize ?
2348             (*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
2349             (*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
2350   }
2351
2352   return os.good() ? COLVARS_OK : FILE_ERROR;
2353 }
2354
2355
2356 int colvar::calc_runave()
2357 {
2358   int error_code = COLVARS_OK;
2359
2360   if (x_history.empty()) {
2361
2362     runave.type(value().type());
2363     runave.reset();
2364
2365     // first-step operationsf
2366
2367     if (cvm::debug())
2368       cvm::log("Colvar \""+this->name+
2369                 "\": initializing running average calculation.\n");
2370
2371     acf_nframes = 0;
2372
2373     x_history.push_back(std::list<colvarvalue>());
2374     x_history_p = x_history.begin();
2375
2376   } else {
2377
2378     if ( (cvm::step_relative() % runave_stride) == 0 &&
2379          (cvm::step_relative() > prev_timestep) ) {
2380
2381       if ((*x_history_p).size() >= runave_length-1) {
2382
2383         if (runave_os == NULL) {
2384           if (runave_outfile.size() == 0) {
2385             runave_outfile = std::string(cvm::output_prefix()+"."+
2386                                          this->name+".runave.traj");
2387           }
2388
2389           size_t const this_cv_width = x.output_width(cvm::cv_width);
2390           cvm::proxy->backup_file(runave_outfile);
2391           runave_os = cvm::proxy->output_stream(runave_outfile);
2392           runave_os->setf(std::ios::scientific, std::ios::floatfield);
2393           *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
2394                      << "   "
2395                      << cvm::wrap_string("running average", this_cv_width)
2396                      << " "
2397                      << cvm::wrap_string("running stddev", this_cv_width)
2398                      << "\n";
2399         }
2400
2401         runave = x;
2402         std::list<colvarvalue>::iterator xs_i;
2403         for (xs_i = (*x_history_p).begin();
2404              xs_i != (*x_history_p).end(); ++xs_i) {
2405           runave += (*xs_i);
2406         }
2407         runave *= 1.0 / cvm::real(runave_length);
2408         runave.apply_constraints();
2409
2410         runave_variance = 0.0;
2411         runave_variance += this->dist2(x, runave);
2412         for (xs_i = (*x_history_p).begin();
2413              xs_i != (*x_history_p).end(); ++xs_i) {
2414           runave_variance += this->dist2(x, (*xs_i));
2415         }
2416         runave_variance *= 1.0 / cvm::real(runave_length-1);
2417
2418         *runave_os << std::setw(cvm::it_width) << cvm::step_relative()
2419                    << "   "
2420                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2421                    << runave << " "
2422                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
2423                    << std::sqrt(runave_variance) << "\n";
2424       }
2425
2426       history_add_value(runave_length, *x_history_p, x);
2427     }
2428   }
2429
2430   return error_code;
2431 }
2432
2433 // Static members
2434
2435 std::vector<colvardeps::feature *> colvar::cv_features;