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