Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvardeps.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
11 #include "colvarmodule.h"
12 #include "colvarproxy.h"
13 #include "colvardeps.h"
14
15
16 colvardeps::colvardeps()
17 {
18   time_step_factor = 1;
19 }
20
21
22 colvardeps::~colvardeps() {
23   size_t i;
24
25   // Protest if we are deleting an object while a parent object may still depend on it
26   if (parents.size()) {
27     cvm::log("Warning: destroying \"" + description + "\" before its parents objects:");
28     for (i=0; i<parents.size(); i++) {
29       cvm::log(parents[i]->description);
30     }
31   }
32
33   // Do not delete features if it's a static object
34   // may change in the future though
35 //     for (i=0; i<features.size(); i++) {
36 //       if (features[i] != NULL) delete features[i];
37 //     }
38
39   remove_all_children();
40 }
41
42
43 void colvardeps::free_children_deps() {
44   // Dereference children requirements of all enabled features
45   // Useful when object is destroyed or set inactive
46   // CAUTION: when setting the parent object inactive, disable "active" first
47   // then call this, to avoid double-dereferencing the deps of "active"
48
49   // Cannot be in the base class destructor because it needs the derived class features()
50   size_t i,j,fid;
51
52   if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description);
53
54   cvm::increase_depth();
55   for (fid = 0; fid < feature_states.size(); fid++) {
56     if (is_enabled(fid)) {
57       for (i=0; i<features()[fid]->requires_children.size(); i++) {
58         int g = features()[fid]->requires_children[i];
59         for (j=0; j<children.size(); j++) {
60           if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
61             + children[j]->features()[g]->description);
62           children[j]->decr_ref_count(g);
63         }
64       }
65     }
66   }
67   cvm::decrease_depth();
68 }
69
70
71 // re-enable children features (and increase ref count accordingly)
72 // So free_children_deps() can be called whenever an object becomes inactive
73 void colvardeps::restore_children_deps() {
74   size_t i,j,fid;
75
76   cvm::increase_depth();
77   for (fid = 0; fid < feature_states.size(); fid++) {
78     if (is_enabled(fid)) {
79       for (i=0; i<features()[fid]->requires_children.size(); i++) {
80         int g = features()[fid]->requires_children[i];
81         for (j=0; j<children.size(); j++) {
82           if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
83             + children[j]->features()[g]->description);
84           children[j]->enable(g, false, false);
85         }
86       }
87     }
88   }
89   cvm::decrease_depth();
90 }
91
92
93 void colvardeps::provide(int feature_id, bool truefalse) {
94   feature_states[feature_id].available = truefalse;
95 }
96
97
98 void colvardeps::set_enabled(int feature_id, bool truefalse) {
99   if (truefalse) {
100     enable(feature_id);
101   } else {
102     disable(feature_id);
103   }
104 }
105
106
107 bool colvardeps::get_keyval_feature(colvarparse *cvp,
108                                     std::string const &conf, char const *key,
109                                     int feature_id, bool const &def_value,
110                                     colvarparse::Parse_Mode const parse_mode)
111 {
112   if (!is_user(feature_id)) {
113     cvm::error("Cannot set feature \"" + features()[feature_id]->description + "\" from user input in \"" + description + "\".\n");
114     return false;
115   }
116   bool value;
117   bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode);
118   if (value) enable(feature_id);
119   return found;
120 }
121
122
123 int colvardeps::enable(int feature_id,
124                        bool dry_run /* default: false */,
125                        // dry_run: fail silently, do not enable if available
126                        // flag is passed recursively to deps of this feature
127                        bool toplevel /* default: true */)
128 // toplevel: false if this is called as part of a chain of dependency resolution
129 // this is used to diagnose failed dependencies by displaying the full stack
130 // only the toplevel dependency will throw a fatal error
131 {
132   int res;
133   size_t i, j;
134   bool ok;
135   feature *f = features()[feature_id];
136   feature_state *fs = &feature_states[feature_id];
137
138   if (cvm::debug()) {
139     cvm::log("DEPS: " + description +
140       (dry_run ? " testing " : " enabling ") +
141       "\"" + f->description +"\"");
142   }
143
144   if (fs->enabled) {
145     if (!(dry_run || toplevel)) {
146       // This is a dependency
147       // Prevent disabling this feature as long
148       // as requirement is enabled
149       fs->ref_count++;
150       if (cvm::debug())
151         cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count));
152     }
153     // Do not try to further resolve deps
154     return COLVARS_OK;
155   }
156
157   std::string feature_type_descr = is_static(feature_id) ? "Static" :
158     (is_dynamic(feature_id) ? "Dynamic" : "User-controlled");
159
160   if (!fs->available) {
161     if (!dry_run) {
162       if (toplevel) {
163         cvm::error("Error: " + feature_type_descr + " feature unavailable: \""
164           + f->description + "\" in " + description + ".");
165       } else {
166         cvm::log(feature_type_descr + " feature unavailable: \""
167           + f->description + "\" in " + description + ".");
168       }
169     }
170     return COLVARS_ERROR;
171   }
172
173   if (!toplevel && !is_dynamic(feature_id)) {
174     if (!dry_run) {
175       cvm::log(feature_type_descr + " feature \"" + f->description
176         + "\" may not be enabled as a dependency in " + description + ".\n");
177     }
178     return COLVARS_ERROR;
179   }
180
181   // 1) enforce exclusions
182   // reminder: exclusions must be mutual for this to work
183   for (i=0; i<f->requires_exclude.size(); i++) {
184     feature *g = features()[f->requires_exclude[i]];
185     if (cvm::debug())
186       cvm::log(f->description + " requires exclude " + g->description);
187     if (is_enabled(f->requires_exclude[i])) {
188       if (!dry_run) {
189         cvm::log("Feature \"" + f->description + "\" is incompatible with \""
190         + g->description + "\" in " + description + ".");
191         if (toplevel) {
192           cvm::error("Error: Failed dependency in " + description + ".");
193         }
194       }
195       return COLVARS_ERROR;
196     }
197   }
198
199   // 2) solve internal deps (self)
200   for (i=0; i<f->requires_self.size(); i++) {
201     if (cvm::debug())
202       cvm::log(f->description + " requires self " + features()[f->requires_self[i]]->description);
203     res = enable(f->requires_self[i], dry_run, false);
204     if (res != COLVARS_OK) {
205       if (!dry_run) {
206         cvm::log("...required by \"" + f->description + "\" in " + description);
207         if (toplevel) {
208           cvm::error("Error: Failed dependency in " + description + ".");
209         }
210       }
211       return res;
212     }
213   }
214
215   // 3) solve internal alternate deps
216   for (i=0; i<f->requires_alt.size(); i++) {
217
218     // test if one is available; if yes, enable and exit w/ success
219     ok = false;
220     for (j=0; j<f->requires_alt[i].size(); j++) {
221       int g = f->requires_alt[i][j];
222       if (cvm::debug())
223         cvm::log(f->description + " requires alt " + features()[g]->description);
224       res = enable(g, true, false);  // see if available
225       if (res == COLVARS_OK) {
226         ok = true;
227         if (!dry_run) {
228           enable(g, false, false); // Require again, for real
229           fs->alternate_refs.push_back(g); // We remember we enabled this
230           // so we can free it if this feature gets disabled
231         }
232         break;
233       }
234     }
235     if (!ok) {
236       if (!dry_run) {
237         cvm::log("\"" + f->description + "\" in " + description
238           + " requires one of the following features, none of which can be enabled:\n");
239         cvm::log("-----------------------------------------\n");
240         cvm::increase_depth();
241         for (j=0; j<f->requires_alt[i].size(); j++) {
242           int g = f->requires_alt[i][j];
243           cvm::log(cvm::to_str(j+1) + ". " + features()[g]->description);
244           enable(g, false, false); // Just for printing error output
245         }
246         cvm::decrease_depth();
247         cvm::log("-----------------------------------------");
248         if (toplevel) {
249           cvm::error("Error: Failed dependency in " + description + ".");
250         }
251       }
252       return COLVARS_ERROR;
253     }
254   }
255
256   // 4) solve deps in children
257   // if the object is inactive, we solve but do not enable: will be enabled
258   // when the object becomes active
259   cvm::increase_depth();
260   for (i=0; i<f->requires_children.size(); i++) {
261     int g = f->requires_children[i];
262     for (j=0; j<children.size(); j++) {
263       res = children[j]->enable(g, dry_run || !is_enabled(), false);
264       if (res != COLVARS_OK) {
265         if (!dry_run) {
266           cvm::log("...required by \"" + f->description + "\" in " + description);
267           if (toplevel) {
268             cvm::error("Error: Failed dependency in " + description + ".");
269           }
270         }
271         return res;
272       }
273     }
274   }
275   cvm::decrease_depth();
276
277   // Actually enable feature only once everything checks out
278   if (!dry_run) {
279     fs->enabled = true;
280     // This should be the only reference
281     if (!toplevel) fs->ref_count = 1;
282     if (feature_id == 0) {
283       // Waking up this object, enable all deps in children
284       restore_children_deps();
285     }
286     do_feature_side_effects(feature_id);
287     if (cvm::debug())
288       cvm::log("DEPS: feature \"" + f->description + "\" in "
289         + description + " enabled, ref_count = 1.");
290   }
291   return COLVARS_OK;
292 }
293
294
295 int colvardeps::disable(int feature_id) {
296   size_t i, j;
297   feature *f = features()[feature_id];
298   feature_state *fs = &feature_states[feature_id];
299
300   if (cvm::debug()) cvm::log("DEPS: disabling feature \""
301       + f->description + "\" in " + description);
302
303   if (fs->enabled == false) {
304     return COLVARS_OK;
305   }
306
307   if (fs->ref_count > 1) {
308     cvm::error("Error: cannot disable feature \"" + f->description
309      + "\" in " + description + " because of " + cvm::to_str(fs->ref_count-1)
310      + " remaining references.\n" );
311     return COLVARS_ERROR;
312   }
313
314   // internal deps (self)
315   for (i=0; i<f->requires_self.size(); i++) {
316     if (cvm::debug()) cvm::log("DEPS: dereferencing self "
317       + features()[f->requires_self[i]]->description);
318     decr_ref_count(f->requires_self[i]);
319   }
320
321   // alternates
322   for (i=0; i<fs->alternate_refs.size(); i++) {
323     if (cvm::debug()) cvm::log("DEPS: dereferencing alt "
324       + features()[fs->alternate_refs[i]]->description);
325     decr_ref_count(fs->alternate_refs[i]);
326   }
327   // Forget these, now that they are dereferenced
328   fs->alternate_refs.clear();
329
330   // deps in children
331   // except if the object is inactive, then children dependencies
332   // have already been dereferenced by this function
333   // (or never referenced if feature was enabled while the object
334   // was inactive)
335   if (is_enabled()) {
336     cvm::increase_depth();
337     for (i=0; i<f->requires_children.size(); i++) {
338       int g = f->requires_children[i];
339       for (j=0; j<children.size(); j++) {
340         if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
341           + children[j]->features()[g]->description);
342         children[j]->decr_ref_count(g);
343       }
344     }
345     cvm::decrease_depth();
346   }
347
348   fs->enabled = false;
349   fs->ref_count = 0;
350   if (feature_id == 0) {
351     // Putting this object to sleep
352     free_children_deps();
353   }
354   return COLVARS_OK;
355 }
356
357 int colvardeps::decr_ref_count(int feature_id) {
358   int &rc = feature_states[feature_id].ref_count;
359   feature *f = features()[feature_id];
360
361   if (cvm::debug())
362       cvm::log("DEPS: decreasing reference count of \"" + f->description
363         + "\" in " + description + ".\n");
364
365   if (rc <= 0) {
366     cvm::error("Error: cannot decrease reference count of feature \"" + f->description
367       +  "\" in " + description + ", which is " + cvm::to_str(rc) + ".\n");
368     return COLVARS_ERROR;
369   }
370
371   rc--;
372   if (rc == 0 && f->is_dynamic()) {
373     // we can auto-disable this feature
374     if (cvm::debug())
375       cvm::log("DEPS will now auto-disable dynamic feature \"" + f->description
376      + "\" in " + description + ".\n");
377     disable(feature_id);
378   }
379   return COLVARS_OK;
380 }
381
382 void colvardeps::init_feature(int feature_id, const char *description, feature_type type) {
383   modify_features()[feature_id]->description = description;
384   modify_features()[feature_id]->type = type;
385 }
386
387 // Shorthand macros for describing dependencies
388 #define f_req_self(f, g) features()[f]->requires_self.push_back(g)
389 // This macro ensures that exclusions are symmetric
390 #define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \
391                             features()[g]->requires_exclude.push_back(f)
392 #define f_req_children(f, g) features()[f]->requires_children.push_back(g)
393 #define f_req_alt2(f, g, h) features()[f]->requires_alt.push_back(std::vector<int>(2));\
394   features()[f]->requires_alt.back()[0] = g;                                           \
395   features()[f]->requires_alt.back()[1] = h
396 #define f_req_alt3(f, g, h, i) features()[f]->requires_alt.push_back(std::vector<int>(3));\
397   features()[f]->requires_alt.back()[0] = g;                                           \
398   features()[f]->requires_alt.back()[1] = h;                                           \
399   features()[f]->requires_alt.back()[2] = i
400 #define f_req_alt4(f, g, h, i, j) features()[f]->requires_alt.push_back(std::vector<int>(4));\
401   features()[f]->requires_alt.back()[0] = g;                                           \
402   features()[f]->requires_alt.back()[1] = h;                                           \
403   features()[f]->requires_alt.back()[2] = i;                                           \
404   features()[f]->requires_alt.back()[3] = j
405
406 void colvardeps::init_cvb_requires() {
407   int i;
408   if (features().size() == 0) {
409     for (i = 0; i < f_cvb_ntot; i++) {
410       modify_features().push_back(new feature);
411     }
412
413     init_feature(f_cvb_active, "active", f_type_dynamic);
414     f_req_children(f_cvb_active, f_cv_active);
415
416     init_feature(f_cvb_awake, "awake", f_type_static);
417     f_req_self(f_cvb_awake, f_cvb_active);
418
419     init_feature(f_cvb_apply_force, "apply force", f_type_user);
420     f_req_children(f_cvb_apply_force, f_cv_gradient);
421
422     init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
423     f_req_children(f_cvb_get_total_force, f_cv_total_force);
424
425     init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user);
426     f_req_self(f_cvb_output_acc_work, f_cvb_apply_force);
427
428     init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
429
430     init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
431
432     init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
433     f_req_children(f_cvb_scalar_variables, f_cv_scalar);
434
435     init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
436
437     init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic);
438     f_req_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
439     f_req_children(f_cvb_calc_ti_samples, f_cv_grid);
440
441     init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user);
442     f_req_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
443
444     init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user);
445     f_req_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
446   }
447
448   // Initialize feature_states for each instance
449   feature_states.reserve(f_cvb_ntot);
450   for (i = 0; i < f_cvb_ntot; i++) {
451     feature_states.push_back(feature_state(true, false));
452     // Most features are available, so we set them so
453     // and list exceptions below
454   }
455
456   // only compute TI samples when deriving from colvarbias_ti
457   feature_states[f_cvb_calc_ti_samples].available = false;
458 }
459
460
461 void colvardeps::init_cv_requires() {
462   size_t i;
463   if (features().size() == 0) {
464     for (i = 0; i < f_cv_ntot; i++) {
465       modify_features().push_back(new feature);
466     }
467
468     init_feature(f_cv_active, "active", f_type_dynamic);
469     // Do not require f_cvc_active in children, as some components may be disabled
470     // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
471     f_req_alt4(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
472
473     init_feature(f_cv_awake, "awake", f_type_static);
474     f_req_self(f_cv_awake, f_cv_active);
475
476     init_feature(f_cv_gradient, "gradient", f_type_dynamic);
477     f_req_children(f_cv_gradient, f_cvc_gradient);
478
479     init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
480     f_req_self(f_cv_collect_gradient, f_cv_gradient);
481     f_req_self(f_cv_collect_gradient, f_cv_scalar);
482     // The following exlusion could be lifted by implementing the feature
483     f_req_exclude(f_cv_collect_gradient, f_cv_scripted);
484
485     init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic);
486
487     // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
488     init_feature(f_cv_total_force, "total force", f_type_dynamic);
489     f_req_alt2(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
490
491     // Deps for explicit total force calculation
492     init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic);
493     f_req_self(f_cv_total_force_calc, f_cv_scalar);
494     f_req_self(f_cv_total_force_calc, f_cv_linear);
495     f_req_children(f_cv_total_force_calc, f_cvc_inv_gradient);
496     f_req_self(f_cv_total_force_calc, f_cv_Jacobian);
497
498     init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic);
499     f_req_self(f_cv_Jacobian, f_cv_scalar);
500     f_req_self(f_cv_Jacobian, f_cv_linear);
501     f_req_children(f_cv_Jacobian, f_cvc_Jacobian);
502
503     init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user);
504     f_req_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
505
506     init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user);
507     f_req_self(f_cv_extended_Lagrangian, f_cv_scalar);
508     f_req_self(f_cv_extended_Lagrangian, f_cv_gradient);
509
510     init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user);
511     f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian);
512
513     init_feature(f_cv_linear, "linear", f_type_static);
514
515     init_feature(f_cv_scalar, "scalar", f_type_static);
516
517     init_feature(f_cv_output_energy, "output energy", f_type_user);
518
519     init_feature(f_cv_output_value, "output value", f_type_user);
520
521     init_feature(f_cv_output_velocity, "output velocity", f_type_user);
522     f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity);
523
524     init_feature(f_cv_output_applied_force, "output applied force", f_type_user);
525
526     init_feature(f_cv_output_total_force, "output total force", f_type_user);
527     f_req_self(f_cv_output_total_force, f_cv_total_force);
528
529     init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
530     f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
531
532     init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
533     f_req_self(f_cv_lower_boundary, f_cv_scalar);
534
535     init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
536     f_req_self(f_cv_upper_boundary, f_cv_scalar);
537
538     init_feature(f_cv_grid, "grid", f_type_dynamic);
539     f_req_self(f_cv_grid, f_cv_lower_boundary);
540     f_req_self(f_cv_grid, f_cv_upper_boundary);
541
542     init_feature(f_cv_runave, "running average", f_type_user);
543
544     init_feature(f_cv_corrfunc, "correlation function", f_type_user);
545
546     init_feature(f_cv_scripted, "scripted", f_type_user);
547
548     init_feature(f_cv_custom_function, "custom function", f_type_user);
549     f_req_exclude(f_cv_custom_function, f_cv_scripted);
550
551     init_feature(f_cv_periodic, "periodic", f_type_static);
552     f_req_self(f_cv_periodic, f_cv_scalar);
553     init_feature(f_cv_scalar, "scalar", f_type_static);
554     init_feature(f_cv_linear, "linear", f_type_static);
555     init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
556
557     // because total forces are obtained from the previous time step,
558     // we cannot (currently) have colvar values and total forces for the same timestep
559     init_feature(f_cv_multiple_ts, "multiple timestep colvar");
560     f_req_exclude(f_cv_multiple_ts, f_cv_total_force_calc);
561   }
562
563   // Initialize feature_states for each instance
564   feature_states.reserve(f_cv_ntot);
565   for (i = 0; i < f_cv_ntot; i++) {
566     feature_states.push_back(feature_state(true, false));
567     // Most features are available, so we set them so
568     // and list exceptions below
569    }
570
571   feature_states[f_cv_fdiff_velocity].available =
572     cvm::main()->proxy->simulation_running();
573 }
574
575
576 void colvardeps::init_cvc_requires() {
577   size_t i;
578   // Initialize static array once and for all
579   if (features().size() == 0) {
580     for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
581       modify_features().push_back(new feature);
582     }
583
584     init_feature(f_cvc_active, "active", f_type_dynamic);
585 //     The dependency below may become useful if we use dynamic atom groups
586 //     f_req_children(f_cvc_active, f_ag_active);
587
588     init_feature(f_cvc_scalar, "scalar", f_type_static);
589
590     init_feature(f_cvc_gradient, "gradient", f_type_dynamic);
591
592     init_feature(f_cvc_implicit_gradient, "implicit gradient", f_type_static);
593     f_req_children(f_cvc_implicit_gradient, f_ag_implicit_gradient);
594
595     init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
596     f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
597
598     init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
599     f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
600     f_req_exclude(f_cvc_debug_gradient, f_cvc_implicit_gradient);
601
602     init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
603     f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
604
605     init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static);
606
607     // init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user);
608
609     // Compute total force on first site only to avoid unwanted
610     // coupling to other colvars (see e.g. Ciccotti et al., 2005)
611     init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user);
612     f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
613
614     init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
615     f_req_self(f_cvc_scalable, f_cvc_scalable_com);
616
617     init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static);
618     f_req_self(f_cvc_scalable_com, f_cvc_com_based);
619
620
621     // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group
622     // f_req_children(f_cvc_scalable, f_ag_scalable);
623     // f_req_children(f_cvc_scalable_com, f_ag_scalable_com);
624   }
625
626   // Initialize feature_states for each instance
627   // default as available, not enabled
628   // except dynamic features which default as unavailable
629   feature_states.reserve(f_cvc_ntot);
630   for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
631     bool avail = is_dynamic(i) ? false : true;
632     feature_states.push_back(feature_state(avail, false));
633   }
634
635   // CVCs are enabled from the start - get disabled based on flags
636   feature_states[f_cvc_active].enabled = true;
637
638   // Features that are implemented by all cvcs by default
639   // Each cvc specifies what other features are available
640   feature_states[f_cvc_active].available = true;
641   feature_states[f_cvc_gradient].available = true;
642
643   // Use minimum-image distances by default
644   feature_states[f_cvc_pbc_minimum_image].enabled = true;
645
646   // Features that are implemented by default if their requirements are
647   feature_states[f_cvc_one_site_total_force].available = true;
648
649   // Features That are implemented only for certain simulation engine configurations
650   feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
651   feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available;
652 }
653
654
655 void colvardeps::init_ag_requires() {
656   size_t i;
657   // Initialize static array once and for all
658   if (features().size() == 0) {
659     for (i = 0; i < f_ag_ntot; i++) {
660       modify_features().push_back(new feature);
661     }
662
663     init_feature(f_ag_active, "active", f_type_dynamic);
664     init_feature(f_ag_center, "translational fit", f_type_static);
665     init_feature(f_ag_rotate, "rotational fit", f_type_static);
666     init_feature(f_ag_fitting_group, "reference positions group", f_type_static);
667     init_feature(f_ag_implicit_gradient, "implicit atom gradient", f_type_dynamic);
668     init_feature(f_ag_fit_gradients, "fit gradients", f_type_user);
669     f_req_exclude(f_ag_fit_gradients, f_ag_implicit_gradient);
670
671     init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
672
673     // parallel calculation implies that we have at least a scalable center of mass,
674     // but f_ag_scalable is kept as a separate feature to deal with future dependencies
675     init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
676     init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static);
677     f_req_self(f_ag_scalable, f_ag_scalable_com);
678
679 //     init_feature(f_ag_min_msd_fit, "minimum MSD fit")
680 //     f_req_self(f_ag_min_msd_fit, f_ag_center)
681 //     f_req_self(f_ag_min_msd_fit, f_ag_rotate)
682 //     f_req_exclude(f_ag_min_msd_fit, f_ag_fitting_group)
683   }
684
685   // Initialize feature_states for each instance
686   // default as unavailable, not enabled
687   feature_states.reserve(f_ag_ntot);
688   for (i = 0; i < colvardeps::f_ag_ntot; i++) {
689     feature_states.push_back(feature_state(false, false));
690   }
691
692   // Features that are implemented (or not) by all atom groups
693   feature_states[f_ag_active].available = true;
694   // f_ag_scalable_com is provided by the CVC iff it is COM-based
695   feature_states[f_ag_scalable_com].available = false;
696   // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
697   feature_states[f_ag_scalable].available = true;
698   feature_states[f_ag_fit_gradients].available = true;
699   feature_states[f_ag_implicit_gradient].available = true;
700 }
701
702
703 void colvardeps::print_state() {
704   size_t i;
705   cvm::log("Enabled features of \"" + description + "\" (with reference count)");
706   for (i = 0; i < feature_states.size(); i++) {
707     if (is_enabled(i))
708       cvm::log("- " + features()[i]->description + " ("
709         + cvm::to_str(feature_states[i].ref_count) + ")");
710   }
711   cvm::increase_depth();
712   for (i=0; i<children.size(); i++) {
713     cvm::log("* child " + cvm::to_str(i+1));
714     children[i]->print_state();
715   }
716   cvm::decrease_depth();
717 }
718
719
720 void colvardeps::add_child(colvardeps *child) {
721
722   children.push_back(child);
723   child->parents.push_back((colvardeps *)this);
724
725   // Solve dependencies of already enabled parent features
726   // in the new child
727
728   size_t i, fid;
729   cvm::increase_depth();
730   for (fid = 0; fid < feature_states.size(); fid++) {
731     if (is_enabled(fid)) {
732       for (i=0; i<features()[fid]->requires_children.size(); i++) {
733         int g = features()[fid]->requires_children[i];
734         if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
735           + child->features()[g]->description);
736         child->enable(g, false, false);
737       }
738     }
739   }
740   cvm::decrease_depth();
741 }
742
743
744 void colvardeps::remove_child(colvardeps *child) {
745   int i;
746   bool found = false;
747
748   for (i = children.size()-1; i>=0; --i) {
749     if (children[i] == child) {
750       children.erase(children.begin() + i);
751       found = true;
752       break;
753     }
754   }
755   if (!found) {
756     cvm::error("Trying to remove missing child reference from " + description + "\n");
757   }
758   found = false;
759   for (i = child->parents.size()-1; i>=0; --i) {
760     if (child->parents[i] == this) {
761       child->parents.erase(child->parents.begin() + i);
762       found = true;
763       break;
764     }
765   }
766   if (!found) {
767     cvm::error("Trying to remove missing parent reference from " + child->description + "\n");
768   }
769 }
770
771
772 void colvardeps::remove_all_children() {
773   size_t i;
774   int j;
775   bool found;
776
777   for (i = 0; i < children.size(); ++i) {
778     found = false;
779     for (j = children[i]->parents.size()-1; j>=0; --j) {
780       if (children[i]->parents[j] == this) {
781         children[i]->parents.erase(children[i]->parents.begin() + j);
782         found = true;
783         break;
784       }
785     }
786     if (!found) {
787       cvm::error("Trying to remove missing parent reference from " + children[i]->description + "\n");
788     }
789   }
790   children.clear();
791 }