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