Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvarbias.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 "colvarmodule.h"
11 #include "colvarproxy.h"
12 #include "colvarvalue.h"
13 #include "colvarbias.h"
14 #include "colvargrid.h"
15
16
17 colvarbias::colvarbias(char const *key)
18   : bias_type(to_lower_cppstr(key))
19 {
20   init_cvb_requires();
21
22   rank = 1;
23
24   has_data = false;
25   b_output_energy = false;
26   reset();
27   state_file_step = 0;
28   description = "uninitialized " + cvm::to_str(key) + " bias";
29 }
30
31
32 int colvarbias::init(std::string const &conf)
33 {
34   colvarparse::init(conf);
35
36   size_t i = 0;
37
38   if (name.size() == 0) {
39
40     // first initialization
41
42     cvm::log("Initializing a new \""+bias_type+"\" instance.\n");
43     rank = cvm::main()->num_biases_type(bias_type);
44     get_keyval(conf, "name", name, bias_type+cvm::to_str(rank));
45
46     {
47       colvarbias *bias_with_name = cvm::bias_by_name(this->name);
48       if (bias_with_name != NULL) {
49         if ((bias_with_name->rank != this->rank) ||
50             (bias_with_name->bias_type != this->bias_type)) {
51           cvm::error("Error: this bias cannot have the same name, \""+this->name+
52                      "\", as another bias.\n", INPUT_ERROR);
53           return INPUT_ERROR;
54         }
55       }
56     }
57
58     description = "bias " + name;
59
60     {
61       // lookup the associated colvars
62       std::vector<std::string> colvar_names;
63       if (get_keyval(conf, "colvars", colvar_names)) {
64         if (num_variables()) {
65           cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n",
66                      INPUT_ERROR);
67           return INPUT_ERROR;
68         }
69         for (i = 0; i < colvar_names.size(); i++) {
70           add_colvar(colvar_names[i]);
71         }
72       }
73     }
74
75     if (!num_variables()) {
76       cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
77       return INPUT_ERROR;
78     }
79   } else {
80     cvm::log("Reinitializing bias \""+name+"\".\n");
81   }
82
83   output_prefix = cvm::output_prefix();
84
85   get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
86
87   get_keyval(conf, "timeStepFactor", time_step_factor, 1);
88   if (time_step_factor < 1) {
89     cvm::error("Error: timeStepFactor must be 1 or greater.\n");
90     return COLVARS_ERROR;
91   }
92
93   // Now that children are defined, we can solve dependencies
94   enable(f_cvb_active);
95   if (cvm::debug()) print_state();
96
97   return COLVARS_OK;
98 }
99
100
101 int colvarbias::reset()
102 {
103   bias_energy = 0.0;
104   for (size_t i = 0; i < num_variables(); i++) {
105     colvar_forces[i].reset();
106   }
107   return COLVARS_OK;
108 }
109
110
111 colvarbias::colvarbias()
112   : colvarparse(), has_data(false)
113 {}
114
115
116 colvarbias::~colvarbias()
117 {
118   colvarbias::clear();
119 }
120
121
122 int colvarbias::clear()
123 {
124   free_children_deps();
125
126   // Remove references to this bias from colvars
127   for (std::vector<colvar *>::iterator cvi = colvars.begin();
128        cvi != colvars.end();
129        ++cvi) {
130     for (std::vector<colvarbias *>::iterator bi = (*cvi)->biases.begin();
131          bi != (*cvi)->biases.end();
132          ++bi) {
133       if ( *bi == this) {
134         (*cvi)->biases.erase(bi);
135         break;
136       }
137     }
138   }
139
140   colvarmodule *cv = cvm::main();
141   // ...and from the colvars module
142   for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
143        bi != cv->biases.end();
144        ++bi) {
145     if ( *bi == this) {
146       cv->biases.erase(bi);
147       break;
148     }
149   }
150
151   return COLVARS_OK;
152 }
153
154
155 int colvarbias::clear_state_data()
156 {
157   // no mutable content to delete for base class
158   return COLVARS_OK;
159 }
160
161
162 int colvarbias::add_colvar(std::string const &cv_name)
163 {
164   if (colvar *cv = cvm::colvar_by_name(cv_name)) {
165
166     if (cvm::debug()) {
167       cvm::log("Applying this bias to collective variable \""+
168                cv->name+"\".\n");
169     }
170
171     colvars.push_back(cv);
172
173     colvar_forces.push_back(colvarvalue());
174     colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
175     colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
176     colvar_forces.back().reset();
177
178     previous_colvar_forces.push_back(colvar_forces.back());
179
180     cv->biases.push_back(this); // add back-reference to this bias to colvar
181
182     if (is_enabled(f_cvb_apply_force)) {
183       cv->enable(f_cv_gradient);
184     }
185
186     // Add dependency link.
187     // All biases need at least the value of each colvar
188     // although possibly not at all timesteps
189     add_child(cv);
190
191   } else {
192     cvm::error("Error: cannot find a colvar named \""+
193                cv_name+"\".\n", INPUT_ERROR);
194     return INPUT_ERROR;
195   }
196
197   return COLVARS_OK;
198 }
199
200
201 int colvarbias::update()
202 {
203   if (cvm::debug()) {
204     cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n");
205   }
206
207   has_data = true;
208
209   bias_energy = 0.0;
210   for (size_t ir = 0; ir < num_variables(); ir++) {
211     colvar_forces[ir].reset();
212   }
213
214   return COLVARS_OK;
215 }
216
217
218 void colvarbias::communicate_forces()
219 {
220   size_t i = 0;
221   for (i = 0; i < num_variables(); i++) {
222     if (cvm::debug()) {
223       cvm::log("Communicating a force to colvar \""+
224                variables(i)->name+"\".\n");
225     }
226     // Impulse-style multiple timestep
227     // Note that biases with different values of time_step_factor
228     // may send forces to the same colvar
229     // which is why rescaling has to happen now: the colvar is not
230     // aware of this bias' time_step_factor
231     variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]);
232   }
233   for (i = 0; i < num_variables(); i++) {
234     previous_colvar_forces[i] = colvar_forces[i];
235   }
236 }
237
238
239 int colvarbias::end_of_step()
240 {
241   return COLVARS_OK;
242 }
243
244
245 int colvarbias::change_configuration(std::string const &conf)
246 {
247   cvm::error("Error: change_configuration() not implemented.\n",
248              COLVARS_NOT_IMPLEMENTED);
249   return COLVARS_NOT_IMPLEMENTED;
250 }
251
252
253 cvm::real colvarbias::energy_difference(std::string const &conf)
254 {
255   cvm::error("Error: energy_difference() not implemented.\n",
256              COLVARS_NOT_IMPLEMENTED);
257   return 0.0;
258 }
259
260
261 // So far, these are only implemented in colvarbias_abf
262 int colvarbias::bin_num()
263 {
264   cvm::error("Error: bin_num() not implemented.\n");
265   return COLVARS_NOT_IMPLEMENTED;
266 }
267 int colvarbias::current_bin()
268 {
269   cvm::error("Error: current_bin() not implemented.\n");
270   return COLVARS_NOT_IMPLEMENTED;
271 }
272 int colvarbias::bin_count(int bin_index)
273 {
274   cvm::error("Error: bin_count() not implemented.\n");
275   return COLVARS_NOT_IMPLEMENTED;
276 }
277 int colvarbias::replica_share()
278 {
279   cvm::error("Error: replica_share() not implemented.\n");
280   return COLVARS_NOT_IMPLEMENTED;
281 }
282
283
284 std::string const colvarbias::get_state_params() const
285 {
286   std::ostringstream os;
287   os << "step " << cvm::step_absolute() << "\n"
288      << "name " << this->name << "\n";
289   return os.str();
290 }
291
292
293 int colvarbias::set_state_params(std::string const &conf)
294 {
295   std::string new_name = "";
296   if (colvarparse::get_keyval(conf, "name", new_name,
297                               std::string(""), colvarparse::parse_silent) &&
298       (new_name != this->name)) {
299     cvm::error("Error: in the state file, the "
300                "\""+bias_type+"\" block has a different name, \""+new_name+
301                "\": different system?\n", INPUT_ERROR);
302   }
303
304   if (name.size() == 0) {
305     cvm::error("Error: \""+bias_type+"\" block within the restart file "
306                "has no identifiers.\n", INPUT_ERROR);
307   }
308
309   colvarparse::get_keyval(conf, "step", state_file_step,
310                           cvm::step_absolute(), colvarparse::parse_silent);
311
312   return COLVARS_OK;
313 }
314
315
316 std::ostream & colvarbias::write_state(std::ostream &os)
317 {
318   if (cvm::debug()) {
319     cvm::log("Writing state file for bias \""+name+"\"\n");
320   }
321   os.setf(std::ios::scientific, std::ios::floatfield);
322   os.precision(cvm::cv_prec);
323   os << bias_type << " {\n"
324      << "  configuration {\n";
325   std::istringstream is(get_state_params());
326   std::string line;
327   while (std::getline(is, line)) {
328     os << "    " << line << "\n";
329   }
330   os << "  }\n";
331   write_state_data(os);
332   os << "}\n\n";
333   return os;
334 }
335
336
337 std::istream & colvarbias::read_state(std::istream &is)
338 {
339   size_t const start_pos = is.tellg();
340
341   std::string key, brace, conf;
342   if ( !(is >> key)   || !(key == bias_type) ||
343        !(is >> brace) || !(brace == "{") ||
344        !(is >> colvarparse::read_block("configuration", conf)) ||
345        (set_state_params(conf) != COLVARS_OK) ) {
346     cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
347                this->name+"\" at position "+
348                cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
349     is.clear();
350     is.seekg(start_pos, std::ios::beg);
351     is.setstate(std::ios::failbit);
352     return is;
353   }
354
355   if (!read_state_data(is)) {
356     cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
357                this->name+"\" at position "+
358                cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
359     is.clear();
360     is.seekg(start_pos, std::ios::beg);
361     is.setstate(std::ios::failbit);
362   }
363
364   is >> brace;
365   if (brace != "}") {
366     cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
367                this->name+"\": no matching brace at position "+
368                cvm::to_str(is.tellg())+" in stream.\n");
369     is.setstate(std::ios::failbit);
370   }
371
372   return is;
373 }
374
375
376 std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
377 {
378   size_t const start_pos = is.tellg();
379   std::string key_in;
380   if ( !(is >> key_in) ||
381        !(key_in == to_lower_cppstr(std::string(key))) ) {
382     cvm::error("Error: in reading restart configuration for "+
383                bias_type+" bias \""+this->name+"\" at position "+
384                cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
385     is.clear();
386     is.seekg(start_pos, std::ios::beg);
387     is.setstate(std::ios::failbit);
388     return is;
389   }
390   return is;
391 }
392
393
394
395 std::ostream & colvarbias::write_traj_label(std::ostream &os)
396 {
397   os << " ";
398   if (b_output_energy)
399     os << " E_"
400        << cvm::wrap_string(this->name, cvm::en_width-2);
401   return os;
402 }
403
404
405 std::ostream & colvarbias::write_traj(std::ostream &os)
406 {
407   os << " ";
408   if (b_output_energy)
409     os << " "
410        << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
411        << bias_energy;
412   return os;
413 }
414
415
416
417 colvarbias_ti::colvarbias_ti(char const *key)
418   : colvarbias(key)
419 {
420   provide(f_cvb_calc_ti_samples);
421   ti_avg_forces = NULL;
422   ti_count = NULL;
423 }
424
425
426 colvarbias_ti::~colvarbias_ti()
427 {
428   colvarbias_ti::clear_state_data();
429 }
430
431
432 int colvarbias_ti::clear_state_data()
433 {
434   if (ti_avg_forces != NULL) {
435     delete ti_avg_forces;
436     ti_avg_forces = NULL;
437   }
438   if (ti_count != NULL) {
439     delete ti_count;
440     ti_count = NULL;
441   }
442   return COLVARS_OK;
443 }
444
445
446 int colvarbias_ti::init(std::string const &conf)
447 {
448   int error_code = COLVARS_OK;
449
450   get_keyval_feature(this, conf, "writeTISamples",
451                      f_cvb_write_ti_samples,
452                      is_enabled(f_cvb_write_ti_samples));
453
454   get_keyval_feature(this, conf, "writeTIPMF",
455                      f_cvb_write_ti_pmf,
456                      is_enabled(f_cvb_write_ti_pmf));
457
458   if ((num_variables() > 1) && is_enabled(f_cvb_write_ti_pmf)) {
459     return cvm::error("Error: only 1-dimensional PMFs can be written "
460                       "on the fly.\n"
461                       "Consider using writeTISamples instead and "
462                       "post-processing the sampled free-energy gradients.\n",
463                       COLVARS_NOT_IMPLEMENTED);
464   } else {
465     error_code |= init_grids();
466   }
467
468   if (is_enabled(f_cvb_write_ti_pmf)) {
469     enable(f_cvb_write_ti_samples);
470   }
471
472   if (is_enabled(f_cvb_calc_ti_samples)) {
473     std::vector<std::string> const time_biases =
474       cvm::main()->time_dependent_biases();
475     if (time_biases.size() > 0) {
476       if ((time_biases.size() > 1) || (time_biases[0] != this->name)) {
477         for (size_t i = 0; i < num_variables(); i++) {
478           if (! variables(i)->is_enabled(f_cv_subtract_applied_force)) {
479             return cvm::error("Error: cannot collect TI samples while other "
480                               "time-dependent biases are active and not all "
481                               "variables have subtractAppliedForces on.\n",
482                               INPUT_ERROR);
483           }
484         }
485       }
486     }
487   }
488
489   return error_code;
490 }
491
492
493 int colvarbias_ti::init_grids()
494 {
495   if (is_enabled(f_cvb_calc_ti_samples)) {
496     if (ti_avg_forces == NULL) {
497       ti_bin.resize(num_variables());
498       ti_system_forces.resize(num_variables());
499       for (size_t icv = 0; icv < num_variables(); icv++) {
500         ti_system_forces[icv].type(variables(icv)->value());
501         ti_system_forces[icv].is_derivative();
502         ti_system_forces[icv].reset();
503       }
504       ti_avg_forces = new colvar_grid_gradient(colvars);
505       ti_count = new colvar_grid_count(colvars);
506       ti_avg_forces->samples = ti_count;
507       ti_count->has_parent_data = true;
508     }
509   }
510
511   return COLVARS_OK;
512 }
513
514
515 int colvarbias_ti::update()
516 {
517   return update_system_forces(NULL);
518 }
519
520
521 int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
522                                         *subtract_forces)
523 {
524   if (! is_enabled(f_cvb_calc_ti_samples)) {
525     return COLVARS_OK;
526   }
527
528   has_data = true;
529
530   if (cvm::debug()) {
531     cvm::log("Updating system forces for bias "+this->name+"\n");
532   }
533
534   colvarproxy *proxy = cvm::main()->proxy;
535
536   size_t i;
537
538   if (proxy->total_forces_same_step()) {
539     for (i = 0; i < num_variables(); i++) {
540       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
541     }
542   }
543
544   // Collect total colvar forces
545   if ((cvm::step_relative() > 0) || proxy->total_forces_same_step()) {
546     if (ti_avg_forces->index_ok(ti_bin)) {
547       for (i = 0; i < num_variables(); i++) {
548         if (variables(i)->is_enabled(f_cv_subtract_applied_force)) {
549           // this colvar is already subtracting all applied forces
550           ti_system_forces[i] = variables(i)->total_force();
551         } else {
552           ti_system_forces[i] = variables(i)->total_force() -
553             ((subtract_forces != NULL) ?
554              (*subtract_forces)[i] : previous_colvar_forces[i]);
555         }
556       }
557       ti_avg_forces->acc_value(ti_bin, ti_system_forces);
558     }
559   }
560
561   if (!proxy->total_forces_same_step()) {
562     // Set the index for use in the next iteration, when total forces come in
563     for (i = 0; i < num_variables(); i++) {
564       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
565     }
566   }
567
568   return COLVARS_OK;
569 }
570
571
572 std::string const colvarbias_ti::get_state_params() const
573 {
574   return std::string("");
575 }
576
577
578 int colvarbias_ti::set_state_params(std::string const &state_conf)
579 {
580   return COLVARS_OK;
581 }
582
583
584 std::ostream & colvarbias_ti::write_state_data(std::ostream &os)
585 {
586   if (! is_enabled(f_cvb_calc_ti_samples)) {
587     return os;
588   }
589   os << "\nhistogram\n";
590   ti_count->write_raw(os);
591   os << "\nsystem_forces\n";
592   ti_avg_forces->write_raw(os);
593   return os;
594 }
595
596
597 std::istream & colvarbias_ti::read_state_data(std::istream &is)
598 {
599   if (! is_enabled(f_cvb_calc_ti_samples)) {
600     return is;
601   }
602   if (cvm::debug()) {
603     cvm::log("Reading state data for the TI estimator.\n");
604   }
605   if (! read_state_data_key(is, "histogram")) {
606     return is;
607   }
608   if (! ti_count->read_raw(is)) {
609     return is;
610   }
611   if (! read_state_data_key(is, "system_forces")) {
612     return is;
613   }
614   if (! ti_avg_forces->read_raw(is)) {
615     return is;
616   }
617   if (cvm::debug()) {
618     cvm::log("Done reading state data for the TI estimator.\n");
619   }
620   return is;
621 }
622
623
624 int colvarbias_ti::write_output_files()
625 {
626   if (!has_data) {
627     // nothing to write
628     return COLVARS_OK;
629   }
630
631   std::string const ti_output_prefix = cvm::output_prefix()+"."+this->name;
632
633   std::ostream *os = NULL;
634
635   if (is_enabled(f_cvb_write_ti_samples)) {
636     std::string const ti_count_file_name(ti_output_prefix+".ti.count");
637     os = cvm::proxy->output_stream(ti_count_file_name);
638     if (os) {
639       ti_count->write_multicol(*os);
640       cvm::proxy->close_output_stream(ti_count_file_name);
641     }
642
643     std::string const ti_grad_file_name(ti_output_prefix+".ti.grad");
644     os = cvm::proxy->output_stream(ti_grad_file_name);
645     if (os) {
646       ti_avg_forces->write_multicol(*os);
647       cvm::proxy->close_output_stream(ti_grad_file_name);
648     }
649   }
650
651   if (is_enabled(f_cvb_write_ti_pmf)) {
652     std::string const pmf_file_name(ti_output_prefix+".ti.pmf");
653     cvm::log("Writing TI PMF to file \""+pmf_file_name+"\".\n");
654     os = cvm::proxy->output_stream(pmf_file_name);
655     if (os) {
656       // get the FE gradient
657       ti_avg_forces->multiply_constant(-1.0);
658       ti_avg_forces->write_1D_integral(*os);
659       ti_avg_forces->multiply_constant(-1.0);
660       cvm::proxy->close_output_stream(pmf_file_name);
661     }
662   }
663
664   return COLVARS_OK;
665 }
666
667
668 // Static members
669
670 std::vector<colvardeps::feature *> colvarbias::cvb_features;