Merge "Actual example for constant pH using a protein with multiple sites"
[namd.git] / colvars / src / colvarmodule.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 <sstream>
11 #include <string.h>
12
13 #include "colvarmodule.h"
14 #include "colvarparse.h"
15 #include "colvarproxy.h"
16 #include "colvar.h"
17 #include "colvarbias.h"
18 #include "colvarbias_abf.h"
19 #include "colvarbias_alb.h"
20 #include "colvarbias_histogram.h"
21 #include "colvarbias_meta.h"
22 #include "colvarbias_restraint.h"
23 #include "colvarscript.h"
24 #include "colvaratoms.h"
25
26
27 colvarmodule::colvarmodule(colvarproxy *proxy_in)
28 {
29   depth_s = 0;
30   cv_traj_os = NULL;
31
32   // pointer to the proxy object
33   if (proxy == NULL) {
34     proxy = proxy_in;
35     parse = new colvarparse();
36   } else {
37     // TODO relax this error to handle multiple molecules in VMD
38     // once the module is not static anymore
39     cvm::error("Error: trying to allocate the collective "
40                "variable module twice.\n", BUG_ERROR);
41     return;
42   }
43
44   cvm::log(cvm::line_marker);
45   cvm::log("Initializing the collective variables module, version "+
46            cvm::to_str(COLVARS_VERSION)+".\n");
47   cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n "
48            "http://dx.doi.org/10.1080/00268976.2013.813594\n"
49            "in any publication based on this calculation.\n");
50
51   if (proxy->smp_enabled() == COLVARS_OK) {
52     cvm::log("SMP parallelism is available.\n");
53   }
54
55   // set initial default values
56
57   // "it_restart" will be set by the input state file, if any;
58   // "it" should be updated by the proxy
59   colvarmodule::it = colvarmodule::it_restart = 0;
60   colvarmodule::it_restart_from_state_file = true;
61
62   colvarmodule::use_scripted_forces = false;
63
64   colvarmodule::b_analysis = false;
65
66   colvarmodule::debug_gradients_step_size = 1.0e-07;
67
68   colvarmodule::rotation::monitor_crossings = false;
69   colvarmodule::rotation::crossing_threshold = 1.0e-02;
70
71   colvarmodule::cv_traj_freq = 100;
72   colvarmodule::restart_out_freq = proxy->restart_frequency();
73
74   // by default overwrite the existing trajectory file
75   colvarmodule::cv_traj_append = false;
76 }
77
78
79 colvarmodule * colvarmodule::main()
80 {
81   return proxy->colvars;
82 }
83
84
85 std::vector<colvar *> *colvarmodule::variables()
86 {
87   return &colvars;
88 }
89
90
91 std::vector<colvar *> *colvarmodule::variables_active()
92 {
93   return &colvars_active;
94 }
95
96
97 std::vector<colvar *> *colvarmodule::variables_active_smp()
98 {
99   return &colvars_smp;
100 }
101
102
103 std::vector<int> *colvarmodule::variables_active_smp_items()
104 {
105   return &colvars_smp_items;
106 }
107
108
109 std::vector<colvarbias *> *colvarmodule::biases_active()
110 {
111   return &(biases_active_);
112 }
113
114
115 size_t colvarmodule::size() const
116 {
117   return colvars.size() + biases.size();
118 }
119
120
121 int colvarmodule::read_config_file(char const  *config_filename)
122 {
123   cvm::log(cvm::line_marker);
124   cvm::log("Reading new configuration from file \""+
125            std::string(config_filename)+"\":\n");
126
127   // open the configfile
128   config_s.open(config_filename);
129   if (!config_s.is_open()) {
130     cvm::error("Error: in opening configuration file \""+
131                std::string(config_filename)+"\".\n",
132                FILE_ERROR);
133     return COLVARS_ERROR;
134   }
135
136   // read the config file into a string
137   std::string conf = "";
138   std::string line;
139   while (colvarparse::getline_nocomments(config_s, line)) {
140     // Delete lines that contain only white space after removing comments
141     if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
142       conf.append(line+"\n");
143   }
144   config_s.close();
145
146   return parse_config(conf);
147 }
148
149
150 int colvarmodule::read_config_string(std::string const &config_str)
151 {
152   cvm::log(cvm::line_marker);
153   cvm::log("Reading new configuration:\n");
154   std::istringstream config_s(config_str);
155
156   // strip the comments away
157   std::string conf = "";
158   std::string line;
159   while (colvarparse::getline_nocomments(config_s, line)) {
160     // Delete lines that contain only white space after removing comments
161     if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
162       conf.append(line+"\n");
163   }
164   return parse_config(conf);
165 }
166
167
168 std::istream & colvarmodule::getline(std::istream &is, std::string &line)
169 {
170   std::string l;
171   if (std::getline(is, l)) {
172     size_t const sz = l.size();
173     if (sz > 0) {
174       if (l[sz-1] == '\r' ) {
175         line = l.substr(0, sz-1);
176       } else {
177         line = l;
178       }
179     } else {
180       line.clear();
181     }
182   }
183   return is;
184 }
185
186
187 int colvarmodule::parse_config(std::string &conf)
188 {
189   extra_conf.clear();
190
191   // parse global options
192   if (catch_input_errors(parse_global_params(conf))) {
193     return get_error();
194   }
195
196   // parse the options for collective variables
197   if (catch_input_errors(parse_colvars(conf))) {
198     return get_error();
199   }
200
201   // parse the options for biases
202   if (catch_input_errors(parse_biases(conf))) {
203     return get_error();
204   }
205
206   // done parsing known keywords, check that all keywords found were valid ones
207   if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
208     return get_error();
209   }
210
211   if (extra_conf.size()) {
212     catch_input_errors(parse_global_params(extra_conf));
213     catch_input_errors(parse_colvars(extra_conf));
214     catch_input_errors(parse_biases(extra_conf));
215     parse->check_keywords(extra_conf, "colvarmodule");
216     extra_conf.clear();
217     if (get_error() != COLVARS_OK) return get_error();
218   }
219
220   cvm::log(cvm::line_marker);
221   cvm::log("Collective variables module (re)initialized.\n");
222   cvm::log(cvm::line_marker);
223
224   // update any necessary proxy data
225   proxy->setup();
226
227   if (cv_traj_os != NULL) {
228     // configuration might have changed, better redo the labels
229     write_traj_label(*cv_traj_os);
230   }
231
232   return get_error();
233 }
234
235
236 int colvarmodule::append_new_config(std::string const &new_conf)
237 {
238   extra_conf += new_conf;
239   return COLVARS_OK;
240 }
241
242
243 int colvarmodule::parse_global_params(std::string const &conf)
244 {
245   colvarmodule *cvm = cvm::main();
246
247   std::string index_file_name;
248   if (parse->get_keyval(conf, "indexFile", index_file_name)) {
249     cvm->read_index_file(index_file_name.c_str());
250   }
251
252   if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
253     if (proxy->b_smp_active == false) {
254       cvm::log("SMP parallelism has been disabled.\n");
255     }
256   }
257
258   parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
259
260   parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
261                     debug_gradients_step_size,
262                     colvarparse::parse_silent);
263
264   parse->get_keyval(conf, "monitorEigenvalueCrossing",
265                     colvarmodule::rotation::monitor_crossings,
266                     colvarmodule::rotation::monitor_crossings,
267                     colvarparse::parse_silent);
268   parse->get_keyval(conf, "eigenvalueCrossingThreshold",
269                     colvarmodule::rotation::crossing_threshold,
270                     colvarmodule::rotation::crossing_threshold,
271                     colvarparse::parse_silent);
272
273   parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq);
274   parse->get_keyval(conf, "colvarsRestartFrequency",
275                     restart_out_freq, restart_out_freq);
276
277   // if this is true when initializing, it means
278   // we are continuing after a reset(): default to true
279   parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append);
280
281   parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);
282
283   parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true);
284
285   if (use_scripted_forces && !proxy->force_script_defined) {
286     cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
287   }
288
289   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
290 }
291
292
293 int colvarmodule::parse_colvars(std::string const &conf)
294 {
295   if (cvm::debug())
296     cvm::log("Initializing the collective variables.\n");
297
298   std::string colvar_conf = "";
299   size_t pos = 0;
300   while (parse->key_lookup(conf, "colvar", &colvar_conf, &pos)) {
301
302     if (colvar_conf.size()) {
303       cvm::log(cvm::line_marker);
304       cvm::increase_depth();
305       colvars.push_back(new colvar());
306       if (((colvars.back())->init(colvar_conf) != COLVARS_OK) ||
307           ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
308         cvm::log("Error while constructing colvar number " +
309                  cvm::to_str(colvars.size()) + " : deleting.");
310         delete colvars.back();  // the colvar destructor updates the colvars array
311         return COLVARS_ERROR;
312       }
313       cvm::decrease_depth();
314     } else {
315       cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR);
316       return COLVARS_ERROR;
317     }
318     cvm::decrease_depth();
319     colvar_conf = "";
320   }
321
322   if (!colvars.size()) {
323     cvm::log("Warning: no collective variables defined.\n");
324   }
325
326   if (colvars.size())
327     cvm::log(cvm::line_marker);
328   cvm::log("Collective variables initialized, "+
329            cvm::to_str(colvars.size())+
330            " in total.\n");
331
332   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
333 }
334
335
336 bool colvarmodule::check_new_bias(std::string &conf, char const *key)
337 {
338   if (cvm::get_error() ||
339       (biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
340     cvm::log("Error while constructing bias number " +
341              cvm::to_str(biases.size()) + " : deleting.\n");
342     delete biases.back(); // the bias destructor updates the biases array
343     return true;
344   }
345   return false;
346 }
347
348
349 template <class bias_type>
350 int colvarmodule::parse_biases_type(std::string const &conf,
351                                     char const *keyword)
352 {
353   std::string bias_conf = "";
354   size_t conf_saved_pos = 0;
355   while (parse->key_lookup(conf, keyword, &bias_conf, &conf_saved_pos)) {
356     if (bias_conf.size()) {
357       cvm::log(cvm::line_marker);
358       cvm::increase_depth();
359       biases.push_back(new bias_type(keyword));
360       biases.back()->init(bias_conf);
361       if (cvm::check_new_bias(bias_conf, keyword) != COLVARS_OK) {
362         return COLVARS_ERROR;
363       }
364       cvm::decrease_depth();
365     } else {
366       cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n",
367                  INPUT_ERROR);
368       return COLVARS_ERROR;
369     }
370     bias_conf = "";
371   }
372   return COLVARS_OK;
373 }
374
375
376 int colvarmodule::parse_biases(std::string const &conf)
377 {
378   if (cvm::debug())
379     cvm::log("Initializing the collective variables biases.\n");
380
381   /// initialize ABF instances
382   parse_biases_type<colvarbias_abf>(conf, "abf");
383
384   /// initialize adaptive linear biases
385   parse_biases_type<colvarbias_alb>(conf, "ALB");
386
387   /// initialize harmonic restraints
388   parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic");
389
390   /// initialize harmonic walls restraints
391   parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls");
392
393   /// initialize histograms
394   parse_biases_type<colvarbias_histogram>(conf, "histogram");
395
396   /// initialize histogram restraints
397   parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint");
398
399   /// initialize linear restraints
400   parse_biases_type<colvarbias_restraint_linear>(conf, "linear");
401
402   /// initialize metadynamics instances
403   parse_biases_type<colvarbias_meta>(conf, "metadynamics");
404
405   if (use_scripted_forces) {
406     cvm::log(cvm::line_marker);
407     cvm::increase_depth();
408     cvm::log("User forces script will be run at each bias update.");
409     cvm::decrease_depth();
410   }
411
412   size_t i;
413
414   size_t n_hist_dep_biases = 0;
415   std::vector<std::string> hist_dep_biases_names;
416   for (i = 0; i < biases.size(); i++) {
417     if (biases[i]->is_enabled(colvardeps::f_cvb_apply_force) &&
418         biases[i]->is_enabled(colvardeps::f_cvb_history_dependent)) {
419       n_hist_dep_biases++;
420       hist_dep_biases_names.push_back(biases[i]->name);
421     }
422   }
423   if (n_hist_dep_biases > 1) {
424     cvm::log("WARNING: there are "+cvm::to_str(n_hist_dep_biases)+
425              " history-dependent biases with non-zero force parameters:\n"+
426              cvm::to_str(hist_dep_biases_names)+"\n"+
427              "Please make sure that their forces do not counteract each other.\n");
428   }
429
430   if (biases.size() || use_scripted_forces) {
431     cvm::log(cvm::line_marker);
432     cvm::log("Collective variables biases initialized, "+
433              cvm::to_str(biases.size())+" in total.\n");
434   } else {
435     if (!use_scripted_forces) {
436       cvm::log("No collective variables biases were defined.\n");
437     }
438   }
439
440   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
441 }
442
443
444 int colvarmodule::num_biases_feature(int feature_id)
445 {
446   colvarmodule *cv = cvm::main();
447   size_t n = 0;
448   for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
449        bi != cv->biases.end();
450        bi++) {
451     if ((*bi)->is_enabled(feature_id)) {
452       n++;
453     }
454   }
455   return n;
456 }
457
458
459 int colvarmodule::num_biases_type(std::string const &type)
460 {
461   colvarmodule *cv = cvm::main();
462   size_t n = 0;
463   for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
464        bi != cv->biases.end();
465        bi++) {
466     if ((*bi)->bias_type == type) {
467       n++;
468     }
469   }
470   return n;
471 }
472
473
474 int colvarmodule::catch_input_errors(int result)
475 {
476   if (result != COLVARS_OK || get_error()) {
477     set_error_bits(result);
478     set_error_bits(INPUT_ERROR);
479     parse->init();
480     return get_error();
481   }
482   return COLVARS_OK;
483 }
484
485
486 colvarbias * colvarmodule::bias_by_name(std::string const &name)
487 {
488   colvarmodule *cv = cvm::main();
489   for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
490        bi != cv->biases.end();
491        bi++) {
492     if ((*bi)->name == name) {
493       return (*bi);
494     }
495   }
496   return NULL;
497 }
498
499
500 colvar *colvarmodule::colvar_by_name(std::string const &name)
501 {
502   colvarmodule *cv = cvm::main();
503   for (std::vector<colvar *>::iterator cvi = cv->colvars.begin();
504        cvi != cv->colvars.end();
505        cvi++) {
506     if ((*cvi)->name == name) {
507       return (*cvi);
508     }
509   }
510   return NULL;
511 }
512
513
514 cvm::atom_group *colvarmodule::atom_group_by_name(std::string const &name)
515 {
516   colvarmodule *cv = cvm::main();
517   for (std::vector<cvm::atom_group *>::iterator agi = cv->named_atom_groups.begin();
518        agi != cv->named_atom_groups.end();
519        agi++) {
520     if ((*agi)->name == name) {
521       return (*agi);
522     }
523   }
524   return NULL;
525 }
526
527
528 int colvarmodule::change_configuration(std::string const &bias_name,
529                                        std::string const &conf)
530 {
531   // This is deprecated; supported strategy is to delete the bias
532   // and parse the new config
533   cvm::increase_depth();
534   colvarbias *b;
535   b = bias_by_name(bias_name);
536   if (b == NULL) {
537     cvm::error("Error: bias not found: " + bias_name);
538     return COLVARS_ERROR;
539   }
540   b->change_configuration(conf);
541   cvm::decrease_depth();
542   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
543 }
544
545
546 std::string colvarmodule::read_colvar(std::string const &name)
547 {
548   cvm::increase_depth();
549   colvar *c;
550   std::stringstream ss;
551   c = colvar_by_name(name);
552   if (c == NULL) {
553     cvm::error("Error: colvar not found: " + name);
554     return std::string();
555   }
556   ss << c->value();
557   cvm::decrease_depth();
558   return ss.str();
559 }
560
561 cvm::real colvarmodule::energy_difference(std::string const &bias_name,
562                                           std::string const &conf)
563 {
564   cvm::increase_depth();
565   colvarbias *b;
566   cvm::real energy_diff = 0.;
567   b = bias_by_name(bias_name);
568   if (b == NULL) {
569     cvm::error("Error: bias not found: " + bias_name);
570     return 0.;
571   }
572   energy_diff = b->energy_difference(conf);
573   cvm::decrease_depth();
574   return energy_diff;
575 }
576
577 int colvarmodule::bias_current_bin(std::string const &bias_name)
578 {
579   cvm::increase_depth();
580   int ret;
581   colvarbias *b = bias_by_name(bias_name);
582
583   if (b != NULL) {
584     ret = b->current_bin();
585   } else {
586     cvm::error("Error: bias not found.\n");
587     ret = COLVARS_ERROR;
588   }
589
590   cvm::decrease_depth();
591   return ret;
592 }
593
594 int colvarmodule::bias_bin_num(std::string const &bias_name)
595 {
596   cvm::increase_depth();
597   int ret;
598   colvarbias *b = bias_by_name(bias_name);
599
600   if (b != NULL) {
601     ret = b->bin_num();
602   } else {
603     cvm::error("Error: bias not found.\n");
604     ret = COLVARS_ERROR;
605   }
606
607   cvm::decrease_depth();
608   return ret;
609 }
610
611 int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index)
612 {
613   cvm::increase_depth();
614   int ret;
615   colvarbias *b = bias_by_name(bias_name);
616
617   if (b != NULL) {
618     ret = b->bin_count(bin_index);
619   } else {
620     cvm::error("Error: bias not found.\n");
621     ret = COLVARS_ERROR;
622   }
623
624   cvm::decrease_depth();
625   return ret;
626 }
627
628 int colvarmodule::bias_share(std::string const &bias_name)
629 {
630   cvm::increase_depth();
631   int ret;
632   colvarbias *b = bias_by_name(bias_name);
633
634   if (b != NULL) {
635     b->replica_share();
636     ret = COLVARS_OK;
637   } else {
638     cvm::error("Error: bias not found.\n");
639     ret = COLVARS_ERROR;
640   }
641
642   cvm::decrease_depth();
643   return ret;
644 }
645
646
647 int colvarmodule::calc()
648 {
649   int error_code = COLVARS_OK;
650
651   if (cvm::debug()) {
652     cvm::log(cvm::line_marker);
653     cvm::log("Collective variables module, step no. "+
654              cvm::to_str(cvm::step_absolute())+"\n");
655   }
656
657   error_code |= calc_colvars();
658   // set biasing forces to zero before biases are calculated and summed over
659   for (std::vector<colvar *>::iterator cvi = colvars.begin();
660        cvi != colvars.end(); cvi++) {
661     (*cvi)->reset_bias_force();
662   }
663   error_code |= calc_biases();
664   error_code |= update_colvar_forces();
665
666   if (cvm::b_analysis) {
667     error_code |= analyze();
668   }
669
670   // write trajectory files, if needed
671   if (cv_traj_freq && cv_traj_name.size()) {
672     error_code |= write_traj_files();
673   }
674
675   // write restart files, if needed
676   if (restart_out_freq && restart_out_name.size()) {
677     error_code |= write_restart_files();
678   }
679
680   return error_code;
681 }
682
683
684 int colvarmodule::calc_colvars()
685 {
686   if (cvm::debug())
687     cvm::log("Calculating collective variables.\n");
688   // calculate collective variables and their gradients
689
690   // First, we need to decide which biases are awake
691   // so they can activate colvars as needed
692   std::vector<colvarbias *>::iterator bi;
693   for (bi = biases.begin(); bi != biases.end(); bi++) {
694     int tsf = (*bi)->get_time_step_factor();
695     if (tsf > 0 && (step_absolute() % tsf == 0)) {
696       (*bi)->enable(colvardeps::f_cvb_awake);
697     } else {
698       (*bi)->disable(colvardeps::f_cvb_awake);
699     }
700   }
701
702   int error_code = COLVARS_OK;
703   std::vector<colvar *>::iterator cvi;
704
705   // Determine which colvars are active at this iteration
706   variables_active()->clear();
707   variables_active()->reserve(variables()->size());
708   for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
709     // Wake up or put to sleep variables
710     int tsf = (*cvi)->get_time_step_factor();
711     if (tsf > 0 && (step_absolute() % tsf == 0)) {
712       (*cvi)->enable(colvardeps::f_cv_awake);
713     } else {
714       (*cvi)->disable(colvardeps::f_cv_awake);
715     }
716
717     if ((*cvi)->is_enabled()) {
718       variables_active()->push_back(*cvi);
719     }
720   }
721
722   // if SMP support is available, split up the work
723   if (proxy->smp_enabled() == COLVARS_OK) {
724
725     // first, calculate how much work (currently, how many active CVCs) each colvar has
726
727     variables_active_smp()->clear();
728     variables_active_smp_items()->clear();
729
730     variables_active_smp()->reserve(variables_active()->size());
731     variables_active_smp_items()->reserve(variables_active()->size());
732
733     // set up a vector containing all components
734     cvm::increase_depth();
735     for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
736
737       error_code |= (*cvi)->update_cvc_flags();
738
739       size_t num_items = (*cvi)->num_active_cvcs();
740       variables_active_smp()->reserve(variables_active_smp()->size() + num_items);
741       variables_active_smp_items()->reserve(variables_active_smp_items()->size() + num_items);
742       for (size_t icvc = 0; icvc < num_items; icvc++) {
743         variables_active_smp()->push_back(*cvi);
744         variables_active_smp_items()->push_back(icvc);
745       }
746     }
747     cvm::decrease_depth();
748
749     // calculate colvar components in parallel
750     error_code |= proxy->smp_colvars_loop();
751
752     cvm::increase_depth();
753     for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
754       error_code |= (*cvi)->collect_cvc_data();
755     }
756     cvm::decrease_depth();
757
758   } else {
759
760     // calculate colvars one at a time
761     cvm::increase_depth();
762     for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
763       error_code |= (*cvi)->calc();
764       if (cvm::get_error()) {
765         return COLVARS_ERROR;
766       }
767     }
768     cvm::decrease_depth();
769   }
770
771   error_code |= cvm::get_error();
772   return error_code;
773 }
774
775
776 int colvarmodule::calc_biases()
777 {
778   // update the biases and communicate their forces to the collective
779   // variables
780   if (cvm::debug() && biases.size())
781     cvm::log("Updating collective variable biases.\n");
782
783   std::vector<colvarbias *>::iterator bi;
784   int error_code = COLVARS_OK;
785
786   // Total bias energy is reset before calling scripted biases
787   total_bias_energy = 0.0;
788
789   // update the list of active biases
790   // which may have changed based on f_cvb_awake in calc_colvars()
791   biases_active()->clear();
792   biases_active()->reserve(biases.size());
793   for (bi = biases.begin(); bi != biases.end(); bi++) {
794     if ((*bi)->is_enabled()) {
795       biases_active()->push_back(*bi);
796     }
797   }
798
799   // if SMP support is available, split up the work
800   if (proxy->smp_enabled() == COLVARS_OK) {
801
802     if (use_scripted_forces && !scripting_after_biases) {
803       // calculate biases and scripted forces in parallel
804       error_code |= proxy->smp_biases_script_loop();
805     } else {
806       // calculate biases in parallel
807       error_code |= proxy->smp_biases_loop();
808     }
809
810   } else {
811
812     if (use_scripted_forces && !scripting_after_biases) {
813       error_code |= calc_scripted_forces();
814     }
815
816     cvm::increase_depth();
817     for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
818       error_code |= (*bi)->update();
819       if (cvm::get_error()) {
820         return error_code;
821       }
822     }
823     cvm::decrease_depth();
824   }
825
826   for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
827     total_bias_energy += (*bi)->get_energy();
828   }
829
830   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
831 }
832
833
834 int colvarmodule::update_colvar_forces()
835 {
836   int error_code = COLVARS_OK;
837
838   std::vector<colvar *>::iterator cvi;
839   std::vector<colvarbias *>::iterator bi;
840
841   // sum the forces from all biases for each collective variable
842   if (cvm::debug() && biases.size())
843     cvm::log("Collecting forces from all biases.\n");
844   cvm::increase_depth();
845   for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
846     (*bi)->communicate_forces();
847     if (cvm::get_error()) {
848       return COLVARS_ERROR;
849     }
850   }
851   cvm::decrease_depth();
852
853   if (use_scripted_forces && scripting_after_biases) {
854     error_code |= calc_scripted_forces();
855   }
856
857   // Now we have collected energies from both built-in and scripted biases
858   if (cvm::debug())
859     cvm::log("Adding total bias energy: " + cvm::to_str(total_bias_energy) + "\n");
860   proxy->add_energy(total_bias_energy);
861
862   cvm::real total_colvar_energy = 0.0;
863   // sum up the forces for each colvar, including wall forces
864   // and integrate any internal
865   // equation of motion (extended system)
866   if (cvm::debug())
867     cvm::log("Updating the internal degrees of freedom "
868              "of colvars (if they have any).\n");
869   cvm::increase_depth();
870   for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
871     // Inactive colvars will only reset their forces and return 0 energy
872     total_colvar_energy += (*cvi)->update_forces_energy();
873     if (cvm::get_error()) {
874       return COLVARS_ERROR;
875     }
876   }
877   cvm::decrease_depth();
878   if (cvm::debug())
879     cvm::log("Adding total colvar energy: " + cvm::to_str(total_colvar_energy) + "\n");
880   proxy->add_energy(total_colvar_energy);
881
882   // make collective variables communicate their forces to their
883   // coupled degrees of freedom (i.e. atoms)
884   if (cvm::debug())
885     cvm::log("Communicating forces from the colvars to the atoms.\n");
886   cvm::increase_depth();
887   for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
888     if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) {
889       (*cvi)->communicate_forces();
890       if (cvm::get_error()) {
891         return COLVARS_ERROR;
892       }
893     }
894   }
895   cvm::decrease_depth();
896
897   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
898 }
899
900
901 int colvarmodule::calc_scripted_forces()
902 {
903   // Run user force script, if provided,
904   // potentially adding scripted forces to the colvars
905   int res;
906   res = proxy->run_force_callback();
907   if (res == COLVARS_NOT_IMPLEMENTED) {
908     cvm::error("Colvar forces scripts are not implemented.");
909     return COLVARS_NOT_IMPLEMENTED;
910   }
911   if (res != COLVARS_OK) {
912     cvm::error("Error running user colvar forces script");
913     return COLVARS_ERROR;
914   }
915   return COLVARS_OK;
916 }
917
918
919 int colvarmodule::write_restart_files()
920 {
921   if ( (cvm::step_relative() > 0) &&
922        ((cvm::step_absolute() % restart_out_freq) == 0) ) {
923     cvm::log("Writing the state file \""+
924              restart_out_name+"\".\n");
925     proxy->backup_file(restart_out_name);
926     std::ostream *restart_out_os = proxy->output_stream(restart_out_name);
927     if (!restart_out_os) return cvm::get_error();
928     if (!write_restart(*restart_out_os)) {
929       return cvm::error("Error: in writing restart file.\n", FILE_ERROR);
930     }
931     proxy->close_output_stream(restart_out_name);
932   }
933
934   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
935 }
936
937
938 int colvarmodule::write_traj_files()
939 {
940   if (cv_traj_os == NULL) {
941     open_traj_file(cv_traj_name);
942   }
943
944   // write labels in the traj file every 1000 lines and at first timestep
945   if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
946     write_traj_label(*cv_traj_os);
947   }
948
949   if ((cvm::step_absolute() % cv_traj_freq) == 0) {
950     write_traj(*cv_traj_os);
951   }
952
953   if (restart_out_freq && (cv_traj_os != NULL)) {
954     // flush the trajectory file if we are at the restart frequency
955     if ( (cvm::step_relative() > 0) &&
956          ((cvm::step_absolute() % restart_out_freq) == 0) ) {
957       cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+
958                cv_traj_name+"\".\n");
959       proxy->flush_output_stream(cv_traj_os);
960     }
961   }
962
963   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
964 }
965
966
967 int colvarmodule::analyze()
968 {
969   if (cvm::debug()) {
970     cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n");
971   }
972
973   if (cvm::step_relative() == 0)
974     cvm::log("Performing analysis.\n");
975
976   // perform colvar-specific analysis
977   for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
978        cvi != variables_active()->end();
979        cvi++) {
980     cvm::increase_depth();
981     (*cvi)->analyze();
982     cvm::decrease_depth();
983   }
984
985   // perform bias-specific analysis
986   for (std::vector<colvarbias *>::iterator bi = biases.begin();
987        bi != biases.end();
988        bi++) {
989     cvm::increase_depth();
990     (*bi)->analyze();
991     cvm::decrease_depth();
992   }
993
994   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
995 }
996
997
998 int colvarmodule::setup()
999 {
1000   if (this->size() == 0) return cvm::get_error();
1001   // loop over all components of all colvars to reset masses of all groups
1002   for (std::vector<colvar *>::iterator cvi = variables()->begin();
1003        cvi != variables()->end();  cvi++) {
1004     (*cvi)->setup();
1005   }
1006   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1007 }
1008
1009
1010 colvarmodule::~colvarmodule()
1011 {
1012   if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) ||
1013       (proxy->smp_thread_id() == 0)) {
1014     reset();
1015     delete parse;
1016     parse = NULL;
1017     proxy = NULL;
1018   }
1019 }
1020
1021
1022 int colvarmodule::reset()
1023 {
1024   parse->init();
1025
1026   cvm::log("Resetting the Collective Variables Module.\n");
1027
1028   // Iterate backwards because we are deleting the elements as we go
1029   for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
1030        bi != biases.rend();
1031        bi++) {
1032     delete *bi; // the bias destructor updates the biases array
1033   }
1034   biases.clear();
1035
1036   // Iterate backwards because we are deleting the elements as we go
1037   for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
1038        cvi != colvars.rend();
1039        cvi++) {
1040     delete *cvi; // the colvar destructor updates the colvars array
1041   }
1042   colvars.clear();
1043
1044   index_groups.clear();
1045   index_group_names.clear();
1046
1047   proxy->reset();
1048
1049   if (cv_traj_os != NULL) {
1050     // Do not close file here, as we might not be done with it yet.
1051     proxy->flush_output_stream(cv_traj_os);
1052   }
1053
1054   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1055 }
1056
1057
1058 int colvarmodule::setup_input()
1059 {
1060   if (this->size() == 0) return cvm::get_error();
1061
1062   std::string restart_in_name("");
1063
1064   // read the restart configuration, if available
1065   if (proxy->input_prefix().size()) {
1066     // read the restart file
1067     restart_in_name = proxy->input_prefix();
1068     std::ifstream input_is(restart_in_name.c_str());
1069     if (!input_is.good()) {
1070       // try by adding the suffix
1071       input_is.clear();
1072       restart_in_name = restart_in_name+std::string(".colvars.state");
1073       input_is.open(restart_in_name.c_str());
1074     }
1075
1076     if (!input_is.good()) {
1077       cvm::error("Error: in opening input file \""+
1078                  std::string(restart_in_name)+"\".\n",
1079                  FILE_ERROR);
1080       return COLVARS_ERROR;
1081     } else {
1082       cvm::log(cvm::line_marker);
1083       cvm::log("Restarting from file \""+restart_in_name+"\".\n");
1084       read_restart(input_is);
1085       if (cvm::get_error() != COLVARS_OK) {
1086         return COLVARS_ERROR;
1087       } else {
1088         proxy->input_prefix().clear();
1089       }
1090       cvm::log(cvm::line_marker);
1091     }
1092   }
1093
1094   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1095 }
1096
1097
1098 int colvarmodule::setup_output()
1099 {
1100   if (this->size() == 0) return cvm::get_error();
1101
1102   int error_code = COLVARS_OK;
1103
1104   // output state file (restart)
1105   restart_out_name = proxy->restart_output_prefix().size() ?
1106     std::string(proxy->restart_output_prefix()+".colvars.state") :
1107     std::string("");
1108
1109   if (restart_out_name.size()) {
1110     cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
1111   }
1112
1113   output_prefix() = proxy->output_prefix();
1114   if (output_prefix().size()) {
1115     cvm::log("The final output state file will be \""+
1116              (output_prefix().size() ?
1117               std::string(output_prefix()+".colvars.state") :
1118               std::string("colvars.state"))+"\".\n");
1119     // cvm::log (cvm::line_marker);
1120   }
1121
1122   cv_traj_name =
1123     (output_prefix().size() ?
1124      std::string(output_prefix()+".colvars.traj") :
1125      std::string(""));
1126
1127   if (cv_traj_freq && cv_traj_name.size()) {
1128     error_code |= open_traj_file(cv_traj_name);
1129   }
1130
1131   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1132        bi != biases.end();
1133        bi++) {
1134     error_code |= (*bi)->setup_output();
1135   }
1136
1137   if (error_code != COLVARS_OK || cvm::get_error()) {
1138     set_error_bits(FILE_ERROR);
1139   }
1140
1141   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1142 }
1143
1144
1145 std::istream & colvarmodule::read_restart(std::istream &is)
1146 {
1147   bool warn_total_forces = false;
1148
1149   {
1150     // read global restart information
1151     std::string restart_conf;
1152     if (is >> colvarparse::read_block("configuration", restart_conf)) {
1153       if (it_restart_from_state_file) {
1154         parse->get_keyval(restart_conf, "step",
1155                           it_restart, (size_t) 0,
1156                           colvarparse::parse_silent);
1157         it = it_restart;
1158       }
1159       std::string restart_version;
1160       parse->get_keyval(restart_conf, "version",
1161                         restart_version, std::string(""),
1162                         colvarparse::parse_silent);
1163       if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {
1164         cvm::log("This state file was generated with version "+restart_version+"\n");
1165       }
1166       if ((restart_version.size() == 0) || (restart_version.compare(std::string(COLVARS_VERSION)) < 0)) {
1167         // check for total force change
1168         if (proxy->total_forces_enabled()) {
1169           warn_total_forces = true;
1170         }
1171       }
1172     }
1173     is.clear();
1174     parse->clear_keyword_registry();
1175   }
1176
1177   // colvars restart
1178   cvm::increase_depth();
1179   for (std::vector<colvar *>::iterator cvi = colvars.begin();
1180        cvi != colvars.end();
1181        cvi++) {
1182     if ( !((*cvi)->read_restart(is)) ) {
1183       cvm::error("Error: in reading restart configuration for collective variable \""+
1184                  (*cvi)->name+"\".\n",
1185                  INPUT_ERROR);
1186     }
1187   }
1188
1189   // biases restart
1190   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1191        bi != biases.end();
1192        bi++) {
1193     if (!((*bi)->read_state(is))) {
1194       cvm::error("Error: in reading restart configuration for bias \""+
1195                  (*bi)->name+"\".\n",
1196                  INPUT_ERROR);
1197     }
1198   }
1199   cvm::decrease_depth();
1200
1201   if (warn_total_forces) {
1202     cvm::log(cvm::line_marker);
1203     cvm::log("WARNING:\n");
1204     std::string const warning("### CHANGES IN THE DEFINITION OF SYSTEM FORCES (NOW TOTAL FORCES)\n\
1205 \n\
1206 Starting from the version 2016-08-10 of the Colvars module, \n\
1207 the role of system forces has been replaced by total forces.\n\
1208 \n\
1209 These include *all* forces acting on a collective variable, whether they\n\
1210 come from the force field potential or from external terms\n\
1211 (e.g. restraints), including forces applied by Colvars itself.\n\
1212 \n\
1213 In NAMD, forces applied by Colvars, IMD, SMD, TMD, symmetry\n\
1214 restraints and tclForces are now all counted in the total force.\n\
1215 \n\
1216 In LAMMPS, forces applied by Colvars itself are now counted in the total\n\
1217 force (all forces from other fixes were being counted already).\n\
1218 \n\
1219 \n\
1220 ### WHEN IS THIS CHANGE RELEVANT\n\
1221 \n\
1222 This change affects results *only* when (1) outputSystemForce is\n\
1223 requested or (2) the ABF bias is used.  All other usage cases are\n\
1224 *unaffected* (colvar restraints, metadynamics, etc).\n\
1225 \n\
1226 When system forces are reported (flag: outputSystemForce), their values\n\
1227 in the output may change, but the physical trajectory is never affected.\n\
1228 The physical results of ABF calculations may be affected in some cases.\n\
1229 \n\
1230 \n\
1231 ### CHANGES TO ABF CALCULATIONS\n\
1232 \n\
1233 Compared to previous Colvars versions, the ABF method will now attempt\n\
1234 to cancel external forces (for example, boundary walls) and it may be\n\
1235 not possible to resume through a state file a simulation that was\n\
1236 performed with a previous version.\n\
1237 \n\
1238 There are three possible scenarios:\n\
1239 \n\
1240 1. No external forces are applied to the atoms used by ABF: results are\n\
1241 unchanged.\n\
1242 \n\
1243 2. Some of the atoms used by ABF experience external forces, but these\n\
1244 forces are not applied directly to the variables used by ABF\n\
1245 (e.g. another colvar that uses the same atoms, tclForces, etc): in this\n\
1246 case, we recommend beginning a new simulation.\n\
1247 \n\
1248 3. External forces are applied to one or more of the colvars used by\n\
1249 ABF, but no other forces are applied to their atoms: you may use the\n\
1250 subtractAppliedForce keyword inside the corresponding colvars to\n\
1251 continue the previous simulation.\n\n");
1252     cvm::log(warning);
1253     cvm::log(cvm::line_marker);
1254
1255     // update this ahead of time in this special case
1256     output_prefix() = proxy->input_prefix();
1257     cvm::log("All output files will now be saved with the prefix \""+output_prefix()+".tmp.*\".\n");
1258     cvm::log(cvm::line_marker);
1259     cvm::log("Please review the important warning above. After that, you may rename:\n\
1260 \""+output_prefix()+".tmp.colvars.state\"\n\
1261 to:\n\
1262 \""+ proxy->input_prefix()+".colvars.state\"\n");
1263     output_prefix() = output_prefix()+".tmp";
1264     write_output_files();
1265     cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
1266   }
1267
1268   return is;
1269 }
1270
1271
1272 int colvarmodule::backup_file(char const *filename)
1273 {
1274   return proxy->backup_file(filename);
1275 }
1276
1277
1278 int colvarmodule::write_output_files()
1279 {
1280   // if this is a simulation run (i.e. not a postprocessing), output data
1281   // must be written to be able to restart the simulation
1282   std::string const out_name =
1283     (output_prefix().size() ?
1284      std::string(output_prefix()+".colvars.state") :
1285      std::string("colvars.state"));
1286   cvm::log("Saving collective variables state to \""+out_name+"\".\n");
1287
1288   std::ostream * os = proxy->output_stream(out_name);
1289   os->setf(std::ios::scientific, std::ios::floatfield);
1290   this->write_restart(*os);
1291   proxy->close_output_stream(out_name);
1292
1293   cvm::increase_depth();
1294   for (std::vector<colvar *>::iterator cvi = colvars.begin();
1295        cvi != colvars.end();
1296        cvi++) {
1297     (*cvi)->write_output_files();
1298   }
1299   cvm::decrease_depth();
1300
1301   cvm::increase_depth();
1302   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1303        bi != biases.end();
1304        bi++) {
1305     (*bi)->write_output_files();
1306     (*bi)->write_state_to_replicas();
1307   }
1308   cvm::decrease_depth();
1309
1310   if (cv_traj_os != NULL) {
1311     // do not close, there may be another run command
1312     proxy->flush_output_stream(cv_traj_os);
1313   }
1314
1315   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1316 }
1317
1318
1319
1320 int colvarmodule::read_traj(char const *traj_filename,
1321                             long        traj_read_begin,
1322                             long        traj_read_end)
1323 {
1324   cvm::log("Opening trajectory file \""+
1325            std::string(traj_filename)+"\".\n");
1326   std::ifstream traj_is(traj_filename);
1327
1328   while (true) {
1329     while (true) {
1330
1331       std::string line("");
1332
1333       do {
1334         if (!colvarparse::getline_nocomments(traj_is, line)) {
1335           cvm::log("End of file \""+std::string(traj_filename)+
1336                    "\" reached, or corrupted file.\n");
1337           traj_is.close();
1338           return false;
1339         }
1340       } while (line.find_first_not_of(colvarparse::white_space) == std::string::npos);
1341
1342       std::istringstream is(line);
1343
1344       if (!(is >> it)) return false;
1345
1346       if ( (it < traj_read_begin) ) {
1347
1348         if ((it % 1000) == 0)
1349           std::cerr << "Skipping trajectory step " << it
1350                     << "                    \r";
1351
1352         continue;
1353
1354       } else {
1355
1356         if ((it % 1000) == 0)
1357           std::cerr << "Reading from trajectory, step = " << it
1358                     << "                    \r";
1359
1360         if ( (traj_read_end > traj_read_begin) &&
1361              (it > traj_read_end) ) {
1362           std::cerr << "\n";
1363           cvm::error("Reached the end of the trajectory, "
1364                      "read_end = "+cvm::to_str(traj_read_end)+"\n",
1365                      FILE_ERROR);
1366           return COLVARS_ERROR;
1367         }
1368
1369         for (std::vector<colvar *>::iterator cvi = colvars.begin();
1370              cvi != colvars.end();
1371              cvi++) {
1372           if (!(*cvi)->read_traj(is)) {
1373             cvm::error("Error: in reading colvar \""+(*cvi)->name+
1374                        "\" from trajectory file \""+
1375                        std::string(traj_filename)+"\".\n",
1376                        FILE_ERROR);
1377             return COLVARS_ERROR;
1378           }
1379         }
1380
1381         break;
1382       }
1383     }
1384   }
1385   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1386 }
1387
1388
1389 std::ostream & colvarmodule::write_restart(std::ostream &os)
1390 {
1391   os.setf(std::ios::scientific, std::ios::floatfield);
1392   os << "configuration {\n"
1393      << "  step " << std::setw(it_width)
1394      << it << "\n"
1395      << "  dt " << dt() << "\n"
1396      << "  version " << std::string(COLVARS_VERSION) << "\n"
1397      << "}\n\n";
1398
1399   int error_code = COLVARS_OK;
1400
1401   cvm::increase_depth();
1402   for (std::vector<colvar *>::iterator cvi = colvars.begin();
1403        cvi != colvars.end();
1404        cvi++) {
1405     (*cvi)->write_restart(os);
1406     error_code |= (*cvi)->write_output_files();
1407   }
1408
1409   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1410        bi != biases.end();
1411        bi++) {
1412     (*bi)->write_state(os);
1413     error_code |= (*bi)->write_state_to_replicas();
1414     error_code |= (*bi)->write_output_files();
1415   }
1416   cvm::decrease_depth();
1417
1418   if (error_code != COLVARS_OK) {
1419     // TODO make this function return an int instead
1420     os.setstate(std::ios::failbit);
1421   }
1422
1423   return os;
1424 }
1425
1426
1427 int colvarmodule::open_traj_file(std::string const &file_name)
1428 {
1429   if (cv_traj_os != NULL) {
1430     return COLVARS_OK;
1431   }
1432
1433   // (re)open trajectory file
1434   if (cv_traj_append) {
1435     cvm::log("Appending to colvar trajectory file \""+file_name+
1436              "\".\n");
1437     cv_traj_os = (cvm::proxy)->output_stream(file_name, std::ios::app);
1438   } else {
1439     cvm::log("Writing to colvar trajectory file \""+file_name+
1440              "\".\n");
1441     proxy->backup_file(file_name.c_str());
1442     cv_traj_os = (cvm::proxy)->output_stream(file_name);
1443   }
1444
1445   if (cv_traj_os == NULL) {
1446     cvm::error("Error: cannot write to file \""+file_name+"\".\n",
1447                FILE_ERROR);
1448   }
1449
1450   return cvm::get_error();
1451 }
1452
1453
1454 int colvarmodule::close_traj_file()
1455 {
1456   if (cv_traj_os != NULL) {
1457     proxy->close_output_stream(cv_traj_name);
1458     cv_traj_os = NULL;
1459   }
1460   return cvm::get_error();
1461 }
1462
1463
1464 std::ostream & colvarmodule::write_traj_label(std::ostream &os)
1465 {
1466   os.setf(std::ios::scientific, std::ios::floatfield);
1467
1468   os << "# " << cvm::wrap_string("step", cvm::it_width-2)
1469      << " ";
1470
1471   cvm::increase_depth();
1472   for (std::vector<colvar *>::iterator cvi = colvars.begin();
1473        cvi != colvars.end();
1474        cvi++) {
1475     (*cvi)->write_traj_label(os);
1476   }
1477   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1478        bi != biases.end();
1479        bi++) {
1480     (*bi)->write_traj_label(os);
1481   }
1482   os << "\n";
1483
1484   if (cvm::debug()) {
1485     proxy->flush_output_stream(&os);
1486   }
1487
1488   cvm::decrease_depth();
1489   return os;
1490 }
1491
1492
1493 std::ostream & colvarmodule::write_traj(std::ostream &os)
1494 {
1495   os.setf(std::ios::scientific, std::ios::floatfield);
1496
1497   os << std::setw(cvm::it_width) << it
1498      << " ";
1499
1500   cvm::increase_depth();
1501   for (std::vector<colvar *>::iterator cvi = colvars.begin();
1502        cvi != colvars.end();
1503        cvi++) {
1504     (*cvi)->write_traj(os);
1505   }
1506   for (std::vector<colvarbias *>::iterator bi = biases.begin();
1507        bi != biases.end();
1508        bi++) {
1509     (*bi)->write_traj(os);
1510   }
1511   os << "\n";
1512
1513   if (cvm::debug()) {
1514     proxy->flush_output_stream(&os);
1515   }
1516
1517   cvm::decrease_depth();
1518   return os;
1519 }
1520
1521
1522 void cvm::log(std::string const &message)
1523 {
1524   // allow logging when the module is not fully initialized
1525   size_t const d = (cvm::main() != NULL) ? depth() : 0;
1526   if (d > 0)
1527     proxy->log((std::string(2*d, ' '))+message);
1528   else
1529     proxy->log(message);
1530 }
1531
1532
1533 void cvm::increase_depth()
1534 {
1535   (depth())++;
1536 }
1537
1538
1539 void cvm::decrease_depth()
1540 {
1541   if (depth() > 0) {
1542     (depth())--;
1543   }
1544 }
1545
1546
1547 size_t & cvm::depth()
1548 {
1549   // NOTE: do not call log() or error() here, to avoid recursion
1550   colvarmodule *cv = cvm::main();
1551   if (proxy->smp_enabled() == COLVARS_OK) {
1552     int const nt = proxy->smp_num_threads();
1553     if (int(cv->depth_v.size()) != nt) {
1554       proxy->smp_lock();
1555       // update array of depths
1556       if (cv->depth_v.size() > 0) { cv->depth_s = cv->depth_v[0]; }
1557       cv->depth_v.clear();
1558       cv->depth_v.assign(nt, cv->depth_s);
1559       proxy->smp_unlock();
1560     }
1561     return cv->depth_v[proxy->smp_thread_id()];
1562   }
1563   return cv->depth_s;
1564 }
1565
1566
1567 void colvarmodule::set_error_bits(int code)
1568 {
1569   if (code < 0) {
1570     cvm::fatal_error("Error: set_error_bits() received negative error code.\n");
1571     return;
1572   }
1573   proxy->smp_lock();
1574   errorCode |= code | COLVARS_ERROR;
1575   proxy->smp_unlock();
1576 }
1577
1578 bool colvarmodule::get_error_bit(int code)
1579 {
1580   return bool(errorCode & code);
1581 }
1582
1583 void colvarmodule::clear_error()
1584 {
1585   proxy->smp_lock();
1586   errorCode = COLVARS_OK;
1587   proxy->smp_unlock();
1588 }
1589
1590
1591 int colvarmodule::error(std::string const &message, int code)
1592 {
1593   set_error_bits(code);
1594   proxy->error(message);
1595   return get_error();
1596 }
1597
1598
1599 int colvarmodule::fatal_error(std::string const &message)
1600 {
1601   set_error_bits(FATAL_ERROR);
1602   proxy->fatal_error(message);
1603   return get_error();
1604 }
1605
1606
1607 int cvm::read_index_file(char const *filename)
1608 {
1609   std::ifstream is(filename, std::ios::binary);
1610   if (!is.good()) {
1611     cvm::error("Error: in opening index file \""+
1612                std::string(filename)+"\".\n",
1613                FILE_ERROR);
1614   }
1615
1616   while (is.good()) {
1617     char open, close;
1618     std::string group_name;
1619     if ( (is >> open) && (open == '[') &&
1620          (is >> group_name) &&
1621          (is >> close) && (close == ']') ) {
1622       for (std::list<std::string>::iterator names_i = index_group_names.begin();
1623            names_i != index_group_names.end();
1624            names_i++) {
1625         if (*names_i == group_name) {
1626           cvm::error("Error: the group name \""+group_name+
1627                      "\" appears in multiple index files.\n",
1628                      FILE_ERROR);
1629         }
1630       }
1631       index_group_names.push_back(group_name);
1632       index_groups.push_back(std::vector<int>());
1633     } else {
1634       cvm::error("Error: in parsing index file \""+
1635                  std::string(filename)+"\".\n",
1636                  INPUT_ERROR);
1637     }
1638
1639     int atom_number = 1;
1640     size_t pos = is.tellg();
1641     while ( (is >> atom_number) && (atom_number > 0) ) {
1642       (index_groups.back()).push_back(atom_number);
1643       pos = is.tellg();
1644     }
1645     is.clear();
1646     is.seekg(pos, std::ios::beg);
1647     std::string delim;
1648     if ( (is >> delim) && (delim == "[") ) {
1649       // new group
1650       is.clear();
1651       is.seekg(pos, std::ios::beg);
1652     } else {
1653       break;
1654     }
1655   }
1656
1657   cvm::log("The following index groups were read from the index file \""+
1658            std::string(filename)+"\":\n");
1659   std::list<std::string>::iterator names_i = index_group_names.begin();
1660   std::list<std::vector<int> >::iterator lists_i = index_groups.begin();
1661   for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) {
1662     cvm::log("  "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n");
1663   }
1664   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1665 }
1666
1667 int cvm::load_atoms(char const *file_name,
1668                     cvm::atom_group &atoms,
1669                     std::string const &pdb_field,
1670                     double const pdb_field_value)
1671 {
1672   return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
1673 }
1674
1675 int cvm::load_coords(char const *file_name,
1676                      std::vector<cvm::atom_pos> &pos,
1677                      const std::vector<int> &indices,
1678                      std::string const &pdb_field,
1679                      double const pdb_field_value)
1680 {
1681   // Differentiate between PDB and XYZ files
1682   // for XYZ files, use CVM internal parser
1683   // otherwise call proxy function for PDB
1684
1685   std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
1686   if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) {
1687     if ( pdb_field.size() > 0 ) {
1688       cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR);
1689       return COLVARS_ERROR;
1690     }
1691     return cvm::load_coords_xyz(file_name, pos, indices);
1692   } else {
1693     return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value);
1694   }
1695 }
1696
1697
1698 int cvm::load_coords_xyz(char const *filename,
1699                          std::vector<atom_pos> &pos,
1700                          const std::vector<int> &indices)
1701 {
1702   std::ifstream xyz_is(filename);
1703   unsigned int natoms;
1704   char symbol[256];
1705   std::string line;
1706
1707   if ( ! (xyz_is >> natoms) ) {
1708     cvm::error("Error: cannot parse XYZ file "
1709                + std::string(filename) + ".\n", INPUT_ERROR);
1710   }
1711   // skip comment line
1712   cvm::getline(xyz_is, line);
1713   cvm::getline(xyz_is, line);
1714   xyz_is.width(255);
1715   std::vector<atom_pos>::iterator pos_i = pos.begin();
1716
1717   if (pos.size() != natoms) { // Use specified indices
1718     int next = 0; // indices are zero-based
1719     std::vector<int>::const_iterator index = indices.begin();
1720     for ( ; pos_i != pos.end() ; pos_i++, index++) {
1721
1722       while ( next < *index ) {
1723         cvm::getline(xyz_is, line);
1724         next++;
1725       }
1726       xyz_is >> symbol;
1727       xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
1728     }
1729   } else {          // Use all positions
1730     for ( ; pos_i != pos.end() ; pos_i++) {
1731       xyz_is >> symbol;
1732       xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
1733     }
1734   }
1735   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
1736 }
1737
1738
1739 // shared pointer to the proxy object
1740 colvarproxy *colvarmodule::proxy = NULL;
1741
1742 // static runtime data
1743 cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
1744 int       colvarmodule::errorCode = 0;
1745 long      colvarmodule::it = 0;
1746 long      colvarmodule::it_restart = 0;
1747 size_t    colvarmodule::restart_out_freq = 0;
1748 size_t    colvarmodule::cv_traj_freq = 0;
1749 bool      colvarmodule::b_analysis = false;
1750 bool      colvarmodule::use_scripted_forces = false;
1751 bool      colvarmodule::scripting_after_biases = true;
1752
1753 // i/o constants
1754 size_t const colvarmodule::it_width = 12;
1755 size_t const colvarmodule::cv_prec  = 14;
1756 size_t const colvarmodule::cv_width = 21;
1757 size_t const colvarmodule::en_prec  = 14;
1758 size_t const colvarmodule::en_width = 21;
1759 const char * const colvarmodule::line_marker = (const char *)
1760   "----------------------------------------------------------------------\n";