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