Update Colvars to version 2018-12-18
[namd.git] / src / colvarproxy_namd.C
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 <errno.h>
11
12 #include "common.h"
13 #include "fstream_namd.h"
14 #include "BackEnd.h"
15 #include "InfoStream.h"
16 #include "Node.h"
17 #include "Molecule.h"
18 #include "PDB.h"
19 #include "PDBData.h"
20 #include "ReductionMgr.h"
21 #include "ScriptTcl.h"
22
23 #ifdef NAMD_TCL
24 #include <tcl.h>
25 #endif
26
27 #include "colvarmodule.h"
28 #include "colvaratoms.h"
29 #include "colvarproxy.h"
30 #include "colvarproxy_namd.h"
31 #include "colvarscript.h"
32
33
34 colvarproxy_namd::colvarproxy_namd()
35 {
36   version_int = get_version_from_string(COLVARPROXY_VERSION);
37
38   first_timestep = true;
39   total_force_requested = false;
40   requestTotalForce(total_force_requested);
41
42   // initialize pointers to NAMD configuration data
43   simparams = Node::Object()->simParameters;
44
45   if (cvm::debug())
46     iout << "Info: initializing the colvars proxy object.\n" << endi;
47
48   // find the configuration file, if provided
49   StringList *config = Node::Object()->configList->find("colvarsConfig");
50
51   // find the input state file
52   StringList *input_restart = Node::Object()->configList->find("colvarsInput");
53   input_prefix_str = std::string(input_restart ? input_restart->data : "");
54   if (input_prefix_str.rfind(".colvars.state") != std::string::npos) {
55     // strip the extension, if present
56     input_prefix_str.erase(input_prefix_str.rfind(".colvars.state"),
57                            std::string(".colvars.state").size());
58   }
59
60   // get the thermostat temperature
61   if (simparams->rescaleFreq > 0)
62     thermostat_temperature = simparams->rescaleTemp;
63   else if (simparams->reassignFreq > 0)
64     thermostat_temperature = simparams->reassignTemp;
65   else if (simparams->langevinOn)
66     thermostat_temperature = simparams->langevinTemp;
67   else if (simparams->tCoupleOn)
68     thermostat_temperature = simparams->tCoupleTemp;
69   //else if (simparams->loweAndersenOn)
70   //  thermostat_temperature = simparams->loweAndersenTemp;
71   else
72     thermostat_temperature = 0.0;
73
74   random = Random(simparams->randomSeed);
75
76   // take the output prefixes from the namd input
77   output_prefix_str = std::string(simparams->outputFilename);
78   restart_output_prefix_str = std::string(simparams->restartFilename);
79   restart_frequency_s = simparams->restartFrequency;
80
81   // check if it is possible to save output configuration
82   if ((!output_prefix_str.size()) && (!restart_output_prefix_str.size())) {
83     fatal_error("Error: neither the final output state file or "
84                 "the output restart file could be defined, exiting.\n");
85   }
86
87
88 #ifdef NAMD_TCL
89   have_scripts = true;
90
91   init_tcl_pointers();
92
93   // See is user-scripted forces are defined
94   if (Tcl_FindCommand(reinterpret_cast<Tcl_Interp *>(tcl_interp_),
95                       "calc_colvar_forces", NULL, 0) == NULL) {
96     force_script_defined = false;
97   } else {
98     force_script_defined = true;
99   }
100 #else
101   force_script_defined = false;
102   have_scripts = false;
103 #endif
104
105
106   // initiate module: this object will be the communication proxy
107   colvars = new colvarmodule(this);
108   log("Using NAMD interface, version "+
109       cvm::to_str(COLVARPROXY_VERSION)+".\n");
110
111   if (config) {
112     colvars->read_config_file(config->data);
113   }
114
115   colvars->setup();
116   colvars->setup_input();
117   colvars->setup_output();
118
119   // save to Node for Tcl script access
120   Node::Object()->colvars = colvars;
121
122 #ifdef NAMD_TCL
123   // Construct instance of colvars scripting interface
124   script = new colvarscript(this);
125 #endif
126
127   if (simparams->firstTimestep != 0) {
128     log("Initializing step number as firstTimestep.\n");
129     colvars->it = colvars->it_restart = simparams->firstTimestep;
130   }
131
132   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
133
134   if (cvm::debug())
135     iout << "Info: done initializing the colvars proxy object.\n" << endi;
136 }
137
138
139 colvarproxy_namd::~colvarproxy_namd()
140 {
141   delete reduction;
142   if (script != NULL) {
143     delete script;
144     script = NULL;
145   }
146   if (colvars != NULL) {
147     delete colvars;
148     colvars = NULL;
149   }
150 }
151
152
153 int colvarproxy_namd::update_atoms_map(AtomIDList::const_iterator begin,
154                                        AtomIDList::const_iterator end)
155 {
156   for (AtomIDList::const_iterator a_i = begin; a_i != end; a_i++) {
157
158     if (cvm::debug()) {
159       cvm::log("Updating atoms_map for atom ID "+cvm::to_str(*a_i)+"\n");
160     }
161
162     if (atoms_map[*a_i] >= 0) continue;
163
164     for (size_t i = 0; i < atoms_ids.size(); i++) {
165       if (atoms_ids[i] == *a_i) {
166         atoms_map[*a_i] = i;
167         break;
168       }
169     }
170
171     if (atoms_map[*a_i] < 0) {
172       // this atom is probably managed by another GlobalMaster:
173       // add it here anyway to avoid having to test for array boundaries at each step
174       int const index = add_atom_slot(*a_i);
175       atoms_map[*a_i] = index;
176       update_atom_properties(index);
177     }
178   }
179
180   if (cvm::debug()) {
181     log("atoms_map = "+cvm::to_str(atoms_map)+".\n");
182   }
183
184   return COLVARS_OK;
185 }
186
187
188 int colvarproxy_namd::setup()
189 {
190   if (colvars->size() == 0) return COLVARS_OK;
191
192   log("Updating NAMD interface:\n");
193
194   if (simparams->wrapAll) {
195     cvm::log("Warning: enabling wrapAll can lead to inconsistent results "
196              "for Colvars calculations: please disable wrapAll, "
197              "as is the default option in NAMD.\n");
198   }
199
200   log("updating atomic data ("+cvm::to_str(atoms_ids.size())+" atoms).\n");
201
202   size_t i;
203   for (i = 0; i < atoms_ids.size(); i++) {
204     update_atom_properties(i);
205
206     // zero out mutable arrays
207     atoms_positions[i] = cvm::rvector(0.0, 0.0, 0.0);
208     atoms_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
209     atoms_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
210   }
211
212   size_t n_group_atoms = 0;
213   for (int ig = 0; ig < modifyRequestedGroups().size(); ig++) {
214     n_group_atoms += modifyRequestedGroups()[ig].size();
215   }
216
217   log("updating group data ("+cvm::to_str(atom_groups_ids.size())+" scalable groups, "+
218       cvm::to_str(n_group_atoms)+" atoms in total).\n");
219
220   // Note: groupMassBegin, groupMassEnd may be used here, but they won't work for charges
221   for (int ig = 0; ig < modifyRequestedGroups().size(); ig++) {
222
223     // update mass and charge
224     update_group_properties(ig);
225
226     atom_groups_coms[ig] = cvm::rvector(0.0, 0.0, 0.0);
227     atom_groups_total_forces[ig] = cvm::rvector(0.0, 0.0, 0.0);
228     atom_groups_new_colvar_forces[ig] = cvm::rvector(0.0, 0.0, 0.0);
229   }
230
231   return COLVARS_OK;
232 }
233
234
235 int colvarproxy_namd::reset()
236 {
237   int error_code = COLVARS_OK;
238
239   // Unrequest all atoms and group from NAMD
240   modifyRequestedAtoms().clear();
241   modifyRequestedGroups().clear();
242
243   atoms_map.clear();
244
245   // Clear internal Proxy records
246   error_code |= colvarproxy::reset();
247
248   return error_code;
249 }
250
251
252 void colvarproxy_namd::calculate()
253 {
254   if (first_timestep) {
255
256     this->setup();
257     colvars->setup();
258     colvars->setup_input();
259     colvars->setup_output();
260
261     first_timestep = false;
262
263   } else {
264     // Use the time step number inherited from GlobalMaster
265     if ( step - previous_NAMD_step == 1 ) {
266       colvars->it++;
267     }
268     // Other cases could mean:
269     // - run 0
270     // - beginning of a new run statement
271     // then the internal counter should not be incremented
272   }
273
274   previous_NAMD_step = step;
275
276   {
277     Vector const a = lattice->a();
278     Vector const b = lattice->b();
279     Vector const c = lattice->c();
280     unit_cell_x.set(a.x, a.y, a.z);
281     unit_cell_y.set(b.x, b.y, c.z);
282     unit_cell_z.set(c.x, c.y, c.z);
283   }
284
285   if (!lattice->a_p() && !lattice->b_p() && !lattice->c_p()) {
286     boundaries_type = boundaries_non_periodic;
287     reset_pbc_lattice();
288   } else if (lattice->a_p() && lattice->b_p() && lattice->c_p()) {
289     if (lattice->orthogonal()) {
290       boundaries_type = boundaries_pbc_ortho;
291     } else {
292       boundaries_type = boundaries_pbc_triclinic;
293     }
294     colvarproxy_system::update_pbc_lattice();
295   } else {
296     boundaries_type = boundaries_unsupported;
297   }
298
299   if (cvm::debug()) {
300     log(std::string(cvm::line_marker)+
301         "colvarproxy_namd, step no. "+cvm::to_str(colvars->it)+"\n"+
302         "Updating atomic data arrays.\n");
303   }
304
305   // must delete the forces applied at the previous step: we can do
306   // that because they have already been used and copied to other
307   // memory locations
308   modifyForcedAtoms().clear();
309   modifyAppliedForces().clear();
310
311   // prepare local arrays
312   for (size_t i = 0; i < atoms_ids.size(); i++) {
313     atoms_positions[i] = cvm::rvector(0.0, 0.0, 0.0);
314     atoms_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
315     atoms_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
316   }
317
318   for (size_t i = 0; i < atom_groups_ids.size(); i++) {
319     atom_groups_total_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
320     atom_groups_new_colvar_forces[i] = cvm::rvector(0.0, 0.0, 0.0);
321   }
322
323   // create the atom map if needed
324   size_t const n_all_atoms = Node::Object()->molecule->numAtoms;
325   if (atoms_map.size() != n_all_atoms) {
326     atoms_map.resize(n_all_atoms);
327     atoms_map.assign(n_all_atoms, -1);
328     update_atoms_map(getAtomIdBegin(), getAtomIdEnd());
329   }
330
331   // if new atomic positions or forces have been communicated by other GlobalMasters, add them to the atom map
332   if ((int(atoms_ids.size()) < (getAtomIdEnd() - getAtomIdBegin())) ||
333       (int(atoms_ids.size()) < (getForceIdEnd() - getForceIdBegin()))) {
334     update_atoms_map(getAtomIdBegin(), getAtomIdEnd());
335     update_atoms_map(getForceIdBegin(), getForceIdEnd());
336   }
337
338   {
339     if (cvm::debug()) {
340       log("Updating positions arrays.\n");
341     }
342     size_t n_positions = 0;
343     AtomIDList::const_iterator a_i = getAtomIdBegin();
344     AtomIDList::const_iterator a_e = getAtomIdEnd();
345     PositionList::const_iterator p_i = getAtomPositionBegin();
346
347     for ( ; a_i != a_e; ++a_i, ++p_i ) {
348       atoms_positions[atoms_map[*a_i]] = cvm::rvector((*p_i).x, (*p_i).y, (*p_i).z);
349       n_positions++;
350     }
351
352     // The following had to be relaxed because some atoms may be forced without their position being requested
353     // if (n_positions < atoms_ids.size()) {
354     //   cvm::error("Error: did not receive the positions of all atoms.\n", BUG_ERROR);
355     // }
356   }
357
358   if (total_force_requested && cvm::step_relative() > 0) {
359
360     // sort the force arrays the previous step
361     // (but only do so if there *is* a previous step!)
362
363     {
364       if (cvm::debug()) {
365         log("Updating total forces arrays.\n");
366       }
367       size_t n_total_forces = 0;
368       AtomIDList::const_iterator a_i = getForceIdBegin();
369       AtomIDList::const_iterator a_e = getForceIdEnd();
370       ForceList::const_iterator f_i = getTotalForce();
371
372       for ( ; a_i != a_e; ++a_i, ++f_i ) {
373         atoms_total_forces[atoms_map[*a_i]] = cvm::rvector((*f_i).x, (*f_i).y, (*f_i).z);
374         n_total_forces++;
375       }
376
377       if (n_total_forces < atoms_ids.size()) {
378         cvm::error("Error: total forces were requested, but total forces "
379                    "were not received for all atoms.\n"
380                    "The most probable cause is combination of energy "
381                    "minimization with a biasing method that requires MD (e.g. ABF).\n"
382                    "Always run minimization and ABF separately.", INPUT_ERROR);
383       }
384     }
385
386     {
387       if (cvm::debug()) {
388         log("Updating group total forces arrays.\n");
389       }
390       ForceList::const_iterator f_i = getGroupTotalForceBegin();
391       ForceList::const_iterator f_e = getGroupTotalForceEnd();
392       size_t i = 0;
393       if ((f_e - f_i) != ((int) atom_groups_ids.size())) {
394         cvm::error("Error: total forces were requested for scalable groups, "
395                    "but they are not in the same number from the number of groups.\n"
396                    "The most probable cause is combination of energy "
397                    "minimization with a biasing method that requires MD (e.g. ABF).\n"
398                    "Always run minimization and ABF separately.", INPUT_ERROR);
399       }
400       for ( ; f_i != f_e; f_i++, i++) {
401         atom_groups_total_forces[i] = cvm::rvector((*f_i).x, (*f_i).y, (*f_i).z);
402       }
403     }
404   }
405
406   {
407     if (cvm::debug()) {
408       log("Updating group positions arrays.\n");
409     }
410     // update group data (only coms available so far)
411     size_t ig;
412     // note: getGroupMassBegin() could be used here, but masses and charges
413     // have already been calculated from the last call to setup()
414     PositionList::const_iterator gp_i = getGroupPositionBegin();
415     for (ig = 0; gp_i != getGroupPositionEnd(); gp_i++, ig++) {
416       atom_groups_coms[ig] = cvm::rvector(gp_i->x, gp_i->y, gp_i->z);
417     }
418   }
419
420   if (cvm::debug()) {
421     log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
422     log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
423     log("atoms_masses = "+cvm::to_str(atoms_masses)+"\n");
424     log("atoms_charges = "+cvm::to_str(atoms_charges)+"\n");
425     log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
426     log("atoms_total_forces = "+cvm::to_str(atoms_total_forces)+"\n");
427     log(cvm::line_marker);
428
429     log("atom_groups_ids = "+cvm::to_str(atom_groups_ids)+"\n");
430     log("atom_groups_ncopies = "+cvm::to_str(atom_groups_ncopies)+"\n");
431     log("atom_groups_masses = "+cvm::to_str(atom_groups_masses)+"\n");
432     log("atom_groups_charges = "+cvm::to_str(atom_groups_charges)+"\n");
433     log("atom_groups_coms = "+cvm::to_str(atom_groups_coms)+"\n");
434     log("atom_groups_total_forces = "+cvm::to_str(atom_groups_total_forces)+"\n");
435     log(cvm::line_marker);
436   }
437
438   // call the collective variable module
439   if (colvars->calc() != COLVARS_OK) {
440     cvm::error("Error in the collective variables module.\n", COLVARS_ERROR);
441   }
442
443   if (cvm::debug()) {
444     log(cvm::line_marker);
445     log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
446     log(cvm::line_marker);
447     log("atom_groups_new_colvar_forces = "+cvm::to_str(atom_groups_new_colvar_forces)+"\n");
448     log(cvm::line_marker);
449   }
450
451   // communicate all forces to the MD integrator
452   for (size_t i = 0; i < atoms_ids.size(); i++) {
453     cvm::rvector const &f = atoms_new_colvar_forces[i];
454     modifyForcedAtoms().add(atoms_ids[i]);
455     modifyAppliedForces().add(Vector(f.x, f.y, f.z));
456   }
457
458   {
459     // zero out the applied forces on each group
460     modifyGroupForces().resize(modifyRequestedGroups().size());
461     ForceList::iterator gf_i = modifyGroupForces().begin();
462     for (int ig = 0; gf_i != modifyGroupForces().end(); gf_i++, ig++) {
463       cvm::rvector const &f = atom_groups_new_colvar_forces[ig];
464       *gf_i = Vector(f.x, f.y, f.z);
465     }
466   }
467
468   // send MISC energy
469   reduction->submit();
470
471   // NAMD does not destruct GlobalMaster objects, so we must remember
472   // to write all output files at the end of a run
473   if (step == simparams->N) {
474     colvars->write_restart_file(cvm::output_prefix()+".colvars.state");
475     colvars->write_output_files();
476   }
477 }
478
479
480 // Callback functions
481
482 void colvarproxy_namd::init_tcl_pointers()
483 {
484 #ifdef NAMD_TCL
485   // Store pointer to NAMD's Tcl interpreter
486   tcl_interp_ = reinterpret_cast<void *>(Node::Object()->getScript()->interp);
487 #endif
488 }
489
490 int colvarproxy_namd::run_force_callback()
491 {
492   return colvarproxy::tcl_run_force_callback();
493 }
494
495 int colvarproxy_namd::run_colvar_callback(
496                           std::string const &name,
497                           std::vector<const colvarvalue *> const &cvc_values,
498                           colvarvalue &value)
499 {
500   return colvarproxy::tcl_run_colvar_callback(name, cvc_values, value);
501 }
502
503 int colvarproxy_namd::run_colvar_gradient_callback(
504                           std::string const &name,
505                           std::vector<const colvarvalue *> const &cvc_values,
506                           std::vector<cvm::matrix2d<cvm::real> > &gradient)
507 {
508   return colvarproxy::tcl_run_colvar_gradient_callback(name, cvc_values,
509                                                        gradient);
510 }
511
512
513 void colvarproxy_namd::add_energy(cvm::real energy)
514 {
515   reduction->item(REDUCTION_MISC_ENERGY) += energy;
516 }
517
518 void colvarproxy_namd::request_total_force(bool yesno)
519 {
520   if (cvm::debug()) {
521     cvm::log("colvarproxy_namd::request_total_force()\n");
522   }
523   total_force_requested = yesno;
524   requestTotalForce(total_force_requested);
525   if (cvm::debug()) {
526     cvm::log("colvarproxy_namd::request_total_force() end\n");
527   }
528 }
529
530 void colvarproxy_namd::log(std::string const &message)
531 {
532   std::istringstream is(message);
533   std::string line;
534   while (std::getline(is, line))
535     iout << "colvars: " << line << "\n";
536   iout << endi;
537 }
538
539 void colvarproxy_namd::error(std::string const &message)
540 {
541   // In NAMD, all errors are fatal
542   fatal_error(message);
543 }
544
545
546 void colvarproxy_namd::fatal_error(std::string const &message)
547 {
548   log(message);
549   if (errno) {
550     // log(strerror(errno));
551     NAMD_err("Error in the collective variables module");
552   } else {
553     NAMD_die("Error in the collective variables module: exiting.\n");
554   }
555 }
556
557
558 void colvarproxy_namd::exit(std::string const &message)
559 {
560   log(message);
561   BackEnd::exit();
562 }
563
564
565 int colvarproxy_namd::check_atom_id(int atom_number)
566 {
567   // NAMD's internal numbering starts from zero
568   int const aid = (atom_number-1);
569
570   if (cvm::debug())
571     log("Adding atom "+cvm::to_str(atom_number)+
572         " for collective variables calculation.\n");
573
574   if ( (aid < 0) || (aid >= Node::Object()->molecule->numAtoms) ) {
575     cvm::error("Error: invalid atom number specified, "+
576                cvm::to_str(atom_number)+"\n", INPUT_ERROR);
577     return INPUT_ERROR;
578   }
579
580   return aid;
581 }
582
583
584 int colvarproxy_namd::init_atom(int atom_number)
585 {
586   // save time by checking first whether this atom has been requested before
587   // (this is more common than a non-valid atom number)
588   int aid = (atom_number-1);
589
590   for (size_t i = 0; i < atoms_ids.size(); i++) {
591     if (atoms_ids[i] == aid) {
592       // this atom id was already recorded
593       atoms_ncopies[i] += 1;
594       return i;
595     }
596   }
597
598   aid = check_atom_id(atom_number);
599
600   if (aid < 0) {
601     return INPUT_ERROR;
602   }
603
604   int const index = add_atom_slot(aid);
605   modifyRequestedAtoms().add(aid);
606   update_atom_properties(index);
607   return index;
608 }
609
610
611 int colvarproxy_namd::check_atom_id(cvm::residue_id const &residue,
612                                     std::string const     &atom_name,
613                                     std::string const     &segment_id)
614 {
615   int const aid =
616     (segment_id.size() ?
617      Node::Object()->molecule->get_atom_from_name(segment_id.c_str(),
618                                                   residue,
619                                                   atom_name.c_str()) :
620      Node::Object()->molecule->get_atom_from_name("MAIN",
621                                                   residue,
622                                                   atom_name.c_str()));
623
624   if (aid < 0) {
625     // get_atom_from_name() has returned an error value
626     cvm::error("Error: could not find atom \""+
627                atom_name+"\" in residue "+
628                cvm::to_str(residue)+
629                ( (segment_id != "MAIN") ?
630                  (", segment \""+segment_id+"\"") :
631                  ("") )+
632                "\n", INPUT_ERROR);
633     return INPUT_ERROR;
634   }
635
636   return aid;
637 }
638
639
640
641 /// For AMBER topologies, the segment id is automatically set to
642 /// "MAIN" (the segment id assigned by NAMD's AMBER topology parser),
643 /// and is therefore optional when an AMBER topology is used
644 int colvarproxy_namd::init_atom(cvm::residue_id const &residue,
645                                 std::string const     &atom_name,
646                                 std::string const     &segment_id)
647 {
648   int const aid = check_atom_id(residue, atom_name, segment_id);
649
650   for (size_t i = 0; i < atoms_ids.size(); i++) {
651     if (atoms_ids[i] == aid) {
652       // this atom id was already recorded
653       atoms_ncopies[i] += 1;
654       return i;
655     }
656   }
657
658   if (cvm::debug())
659     log("Adding atom \""+
660         atom_name+"\" in residue "+
661         cvm::to_str(residue)+
662         " (index "+cvm::to_str(aid)+
663         ") for collective variables calculation.\n");
664
665   int const index = add_atom_slot(aid);
666   modifyRequestedAtoms().add(aid);
667   update_atom_properties(index);
668   return index;
669 }
670
671
672 void colvarproxy_namd::clear_atom(int index)
673 {
674   colvarproxy::clear_atom(index);
675   // TODO remove it from GlobalMaster arrays?
676 }
677
678
679 void colvarproxy_namd::update_atom_properties(int index)
680 {
681   Molecule *mol = Node::Object()->molecule;
682   // update mass
683   double const mass = mol->atommass(atoms_ids[index]);
684   if (mass <= 0.001) {
685     this->log("Warning: near-zero mass for atom "+
686               cvm::to_str(atoms_ids[index]+1)+
687               "; expect unstable dynamics if you apply forces to it.\n");
688   }
689   atoms_masses[index] = mass;
690   // update charge
691   atoms_charges[index] = mol->atomcharge(atoms_ids[index]);
692 }
693
694
695 cvm::rvector colvarproxy_namd::position_distance(cvm::atom_pos const &pos1,
696                                                  cvm::atom_pos const &pos2)
697   const
698 {
699   Position const p1(pos1.x, pos1.y, pos1.z);
700   Position const p2(pos2.x, pos2.y, pos2.z);
701   // return p2 - p1
702   Vector const d = this->lattice->delta(p2, p1);
703   return cvm::rvector(d.x, d.y, d.z);
704 }
705
706
707
708 enum e_pdb_field {
709   e_pdb_none,
710   e_pdb_occ,
711   e_pdb_beta,
712   e_pdb_x,
713   e_pdb_y,
714   e_pdb_z,
715   e_pdb_ntot
716 };
717
718
719 e_pdb_field pdb_field_str2enum(std::string const &pdb_field_str)
720 {
721   e_pdb_field pdb_field = e_pdb_none;
722
723   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
724       colvarparse::to_lower_cppstr("O")) {
725     pdb_field = e_pdb_occ;
726   }
727
728   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
729       colvarparse::to_lower_cppstr("B")) {
730     pdb_field = e_pdb_beta;
731   }
732
733   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
734       colvarparse::to_lower_cppstr("X")) {
735     pdb_field = e_pdb_x;
736   }
737
738   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
739       colvarparse::to_lower_cppstr("Y")) {
740     pdb_field = e_pdb_y;
741   }
742
743   if (colvarparse::to_lower_cppstr(pdb_field_str) ==
744       colvarparse::to_lower_cppstr("Z")) {
745     pdb_field = e_pdb_z;
746   }
747
748   if (pdb_field == e_pdb_none) {
749     cvm::error("Error: unsupported PDB field, \""+
750                pdb_field_str+"\".\n", INPUT_ERROR);
751   }
752
753   return pdb_field;
754 }
755
756
757 int colvarproxy_namd::load_coords(char const *pdb_filename,
758                                   std::vector<cvm::atom_pos> &pos,
759                                   const std::vector<int> &indices,
760                                   std::string const &pdb_field_str,
761                                   double const pdb_field_value)
762 {
763   if (pdb_field_str.size() == 0 && indices.size() == 0) {
764     cvm::error("Bug alert: either PDB field should be defined or list of "
765                "atom IDs should be available when loading atom coordinates!\n", BUG_ERROR);
766   }
767
768   e_pdb_field pdb_field_index;
769   bool const use_pdb_field = (pdb_field_str.size() > 0);
770   if (use_pdb_field) {
771     pdb_field_index = pdb_field_str2enum(pdb_field_str);
772   }
773
774   // next index to be looked up in PDB file (if list is supplied)
775   std::vector<int>::const_iterator current_index = indices.begin();
776
777   PDB *pdb = new PDB(pdb_filename);
778   size_t const pdb_natoms = pdb->num_atoms();
779
780   if (pos.size() != pdb_natoms) {
781
782     bool const pos_allocated = (pos.size() > 0);
783
784     size_t ipos = 0, ipdb = 0;
785     for ( ; ipdb < pdb_natoms; ipdb++) {
786
787       if (use_pdb_field) {
788         // PDB field mode: skip atoms with wrong value in PDB field
789         double atom_pdb_field_value = 0.0;
790
791         switch (pdb_field_index) {
792         case e_pdb_occ:
793           atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
794           break;
795         case e_pdb_beta:
796           atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
797           break;
798         case e_pdb_x:
799           atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
800           break;
801         case e_pdb_y:
802           atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
803           break;
804         case e_pdb_z:
805           atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
806           break;
807         default:
808           break;
809         }
810
811         if ( (pdb_field_value) &&
812              (atom_pdb_field_value != pdb_field_value) ) {
813           continue;
814         } else if (atom_pdb_field_value == 0.0) {
815           continue;
816         }
817
818       } else {
819         // Atom ID mode: use predefined atom IDs from the atom group
820         if (((int) ipdb) != *current_index) {
821           // Skip atoms not in the list
822           continue;
823         } else {
824           current_index++;
825         }
826       }
827
828       if (!pos_allocated) {
829         pos.push_back(cvm::atom_pos(0.0, 0.0, 0.0));
830       } else if (ipos >= pos.size()) {
831         cvm::error("Error: the PDB file \""+
832                    std::string(pdb_filename)+
833                    "\" contains coordinates for "
834                    "more atoms than needed.\n", BUG_ERROR);
835       }
836
837       pos[ipos] = cvm::atom_pos((pdb->atom(ipdb))->xcoor(),
838                                 (pdb->atom(ipdb))->ycoor(),
839                                 (pdb->atom(ipdb))->zcoor());
840       ipos++;
841       if (!use_pdb_field && current_index == indices.end())
842         break;
843     }
844
845     if (ipos < pos.size() || (!use_pdb_field && current_index != indices.end())) {
846       size_t n_requested = use_pdb_field ? pos.size() : indices.size();
847       cvm::error("Error: number of matching records in the PDB file \""+
848                  std::string(pdb_filename)+"\" ("+cvm::to_str(ipos)+
849                  ") does not match the number of requested coordinates ("+
850                  cvm::to_str(n_requested)+").\n", INPUT_ERROR);
851       return COLVARS_ERROR;
852     }
853   } else {
854
855     // when the PDB contains exactly the number of atoms of the array,
856     // ignore the fields and just read coordinates
857     for (size_t ia = 0; ia < pos.size(); ia++) {
858       pos[ia] = cvm::atom_pos((pdb->atom(ia))->xcoor(),
859                               (pdb->atom(ia))->ycoor(),
860                               (pdb->atom(ia))->zcoor());
861     }
862   }
863
864   delete pdb;
865   return COLVARS_OK;
866 }
867
868
869 int colvarproxy_namd::load_atoms(char const *pdb_filename,
870                                  cvm::atom_group &atoms,
871                                  std::string const &pdb_field_str,
872                                  double const pdb_field_value)
873 {
874   if (pdb_field_str.size() == 0)
875     cvm::error("Error: must define which PDB field to use "
876                "in order to define atoms from a PDB file.\n", INPUT_ERROR);
877
878   PDB *pdb = new PDB(pdb_filename);
879   size_t const pdb_natoms = pdb->num_atoms();
880
881   e_pdb_field pdb_field_index = pdb_field_str2enum(pdb_field_str);
882
883   for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
884
885     double atom_pdb_field_value = 0.0;
886
887     switch (pdb_field_index) {
888     case e_pdb_occ:
889       atom_pdb_field_value = (pdb->atom(ipdb))->occupancy();
890       break;
891     case e_pdb_beta:
892       atom_pdb_field_value = (pdb->atom(ipdb))->temperaturefactor();
893       break;
894     case e_pdb_x:
895       atom_pdb_field_value = (pdb->atom(ipdb))->xcoor();
896       break;
897     case e_pdb_y:
898       atom_pdb_field_value = (pdb->atom(ipdb))->ycoor();
899       break;
900     case e_pdb_z:
901       atom_pdb_field_value = (pdb->atom(ipdb))->zcoor();
902       break;
903     default:
904       break;
905     }
906
907     if ( (pdb_field_value) &&
908          (atom_pdb_field_value != pdb_field_value) ) {
909       continue;
910     } else if (atom_pdb_field_value == 0.0) {
911       continue;
912     }
913
914     if (atoms.is_enabled(colvardeps::f_ag_scalable)) {
915       atoms.add_atom_id(ipdb);
916     } else {
917       atoms.add_atom(cvm::atom(ipdb+1));
918     }
919   }
920
921   delete pdb;
922   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
923 }
924
925
926 std::ostream * colvarproxy_namd::output_stream(std::string const &output_name,
927                                                std::ios_base::openmode mode)
928 {
929   if (cvm::debug()) {
930     cvm::log("Using colvarproxy_namd::output_stream()\n");
931   }
932   std::list<std::ostream *>::iterator osi  = output_files.begin();
933   std::list<std::string>::iterator    osni = output_stream_names.begin();
934   for ( ; osi != output_files.end(); osi++, osni++) {
935     if (*osni == output_name) {
936       return *osi;
937     }
938   }
939   if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
940     colvarproxy::backup_file(output_name);
941   }
942   ofstream_namd *os = new ofstream_namd(output_name.c_str(), mode);
943   if (!os->is_open()) {
944     cvm::error("Error: cannot write to file \""+output_name+"\".\n",
945                FILE_ERROR);
946     return NULL;
947   }
948   output_stream_names.push_back(output_name);
949   output_files.push_back(os);
950   return os;
951 }
952
953
954 int colvarproxy_namd::flush_output_stream(std::ostream *os)
955 {
956   std::list<std::ostream *>::iterator osi  = output_files.begin();
957   std::list<std::string>::iterator    osni = output_stream_names.begin();
958   for ( ; osi != output_files.end(); osi++, osni++) {
959     if (*osi == os) {
960       ((ofstream_namd *) *osi)->flush();
961       return COLVARS_OK;
962     }
963   }
964   return COLVARS_ERROR;
965 }
966
967
968 int colvarproxy_namd::close_output_stream(std::string const &output_name)
969 {
970   std::list<std::ostream *>::iterator osi  = output_files.begin();
971   std::list<std::string>::iterator    osni = output_stream_names.begin();
972   for ( ; osi != output_files.end(); osi++, osni++) {
973     if (*osni == output_name) {
974       if (((ofstream_namd *) *osi)->is_open()) {
975         ((ofstream_namd *) *osi)->close();
976       }
977       delete *osi;
978       output_files.erase(osi);
979       output_stream_names.erase(osni);
980       return COLVARS_OK;
981     }
982   }
983   return COLVARS_ERROR;
984 }
985
986
987 int colvarproxy_namd::backup_file(char const *filename)
988 {
989   if (std::string(filename).rfind(std::string(".colvars.state")) != std::string::npos) {
990     NAMD_backup_file(filename, ".old");
991   } else {
992     NAMD_backup_file(filename, ".BAK");
993   }
994   return COLVARS_OK;
995 }
996
997
998 char const *colvarproxy_namd::script_obj_to_str(unsigned char *obj)
999 {
1000   if (cvm::debug()) {
1001     cvm::log("Called colvarproxy_namd::script_obj_to_str().\n");
1002   }
1003 #ifdef NAMD_TCL
1004   return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
1005 #else
1006   // This is most likely not going to be executed
1007   return colvarproxy::script_obj_to_str(obj);
1008 #endif
1009 }
1010
1011
1012 std::vector<std::string> colvarproxy_namd::script_obj_to_str_vector(unsigned char *obj)
1013 {
1014   if (cvm::debug()) {
1015     cvm::log("Called colvarproxy_namd::script_obj_to_str_vector().\n");
1016   }
1017   std::vector<std::string> result;
1018 #ifdef NAMD_TCL
1019   Tcl_Interp *interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
1020   Tcl_Obj *tcl_obj = reinterpret_cast<Tcl_Obj *>(obj);
1021   Tcl_Obj **tcl_list_elems = NULL;
1022   int count = 0;
1023   if (Tcl_ListObjGetElements(interp, tcl_obj, &count, &tcl_list_elems) ==
1024       TCL_OK) {
1025     result.reserve(count);
1026     for (int i = 0; i < count; i++) {
1027       result.push_back(Tcl_GetString(tcl_list_elems[i]));
1028     }
1029   } else {
1030     Tcl_SetResult(interp,
1031                   const_cast<char *>("Cannot parse Tcl list."), TCL_STATIC);
1032   }
1033 #endif
1034   return result;
1035 }
1036
1037
1038 int colvarproxy_namd::init_atom_group(std::vector<int> const &atoms_ids)
1039 {
1040   if (cvm::debug())
1041     log("Reguesting from NAMD a group of size "+cvm::to_str(atoms_ids.size())+
1042         " for collective variables calculation.\n");
1043
1044   // Note: modifyRequestedGroups is supposed to be in sync with the colvarproxy arrays,
1045   // and to stay that way during a simulation
1046
1047   // compare this new group to those already allocated inside GlobalMaster
1048   int ig;
1049   for (ig = 0; ig < modifyRequestedGroups().size(); ig++) {
1050     AtomIDList const &namd_group = modifyRequestedGroups()[ig];
1051     bool b_match = true;
1052
1053     if (namd_group.size() != ((int) atoms_ids.size())) {
1054       b_match = false;
1055     } else {
1056       int ia;
1057       for (ia = 0; ia < namd_group.size(); ia++) {
1058         int const aid = atoms_ids[ia];
1059         if (namd_group[ia] != aid) {
1060           b_match = false;
1061           break;
1062         }
1063       }
1064     }
1065
1066     if (b_match) {
1067       if (cvm::debug())
1068         log("Group was already added.\n");
1069       // this group already exists
1070       atom_groups_ncopies[ig] += 1;
1071       return ig;
1072     }
1073   }
1074
1075   // add this group (note: the argument of add_atom_group_slot() is redundant for NAMD, and provided only for consistency)
1076   size_t const index = add_atom_group_slot(atom_groups_ids.size());
1077   modifyRequestedGroups().resize(atom_groups_ids.size());
1078   // the following is done in calculate()
1079   // modifyGroupForces().resize(atom_groups_ids.size());
1080   AtomIDList &namd_group = modifyRequestedGroups()[index];
1081   namd_group.resize(atoms_ids.size());
1082   int const n_all_atoms = Node::Object()->molecule->numAtoms;
1083   for (size_t ia = 0; ia < atoms_ids.size(); ia++) {
1084     int const aid = atoms_ids[ia];
1085     if (cvm::debug())
1086       log("Adding atom "+cvm::to_str(aid+1)+
1087           " for collective variables calculation.\n");
1088     if ( (aid < 0) || (aid >= n_all_atoms) ) {
1089       cvm::error("Error: invalid atom number specified, "+
1090                  cvm::to_str(aid+1)+"\n", INPUT_ERROR);
1091       return -1;
1092     }
1093     namd_group[ia] = aid;
1094   }
1095
1096   update_group_properties(index);
1097
1098   if (cvm::debug()) {
1099     log("Group has index "+cvm::to_str(index)+"\n");
1100     log("modifyRequestedGroups length = "+cvm::to_str(modifyRequestedGroups().size())+
1101         ", modifyGroupForces length = "+cvm::to_str(modifyGroupForces().size())+"\n");
1102   }
1103
1104   return index;
1105 }
1106
1107
1108 void colvarproxy_namd::clear_atom_group(int index)
1109 {
1110   // do nothing, keep the NAMD arrays in sync with the colvarproxy ones
1111   colvarproxy::clear_atom_group(index);
1112 }
1113
1114
1115 int colvarproxy_namd::update_group_properties(int index)
1116 {
1117   AtomIDList const &namd_group = modifyRequestedGroups()[index];
1118   if (cvm::debug()) {
1119     log("Re-calculating total mass and charge for scalable group no. "+cvm::to_str(index+1)+" ("+
1120         cvm::to_str(namd_group.size())+" atoms).\n");
1121   }
1122
1123   cvm::real total_mass = 0.0;
1124   cvm::real total_charge = 0.0;
1125   for (int i = 0; i < namd_group.size(); i++) {
1126     total_mass += Node::Object()->molecule->atommass(namd_group[i]);
1127     total_charge += Node::Object()->molecule->atomcharge(namd_group[i]);
1128   }
1129   atom_groups_masses[index] = total_mass;
1130   atom_groups_charges[index] = total_charge;
1131
1132   if (cvm::debug()) {
1133     log("total mass = "+cvm::to_str(total_mass)+
1134         ", total charge = "+cvm::to_str(total_charge)+"\n");
1135   }
1136
1137   return COLVARS_OK;
1138 }
1139
1140
1141 #if CMK_SMP && USE_CKLOOP // SMP only
1142
1143 void calc_colvars_items_smp(int first, int last, void *result, int paramNum, void *param)
1144 {
1145   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
1146   colvarmodule *cv = proxy->colvars;
1147
1148   cvm::increase_depth();
1149   for (int i = first; i <= last; i++) {
1150     colvar *x = (*(cv->variables_active_smp()))[i];
1151     int x_item = (*(cv->variables_active_smp_items()))[i];
1152     if (cvm::debug()) {
1153       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
1154                "]: calc_colvars_items_smp(), first = "+cvm::to_str(first)+
1155                ", last = "+cvm::to_str(last)+", cv = "+
1156                x->name+", cvc = "+cvm::to_str(x_item)+"\n");
1157     }
1158     x->calc_cvcs(x_item, 1);
1159   }
1160   cvm::decrease_depth();
1161 }
1162
1163
1164 int colvarproxy_namd::smp_colvars_loop()
1165 {
1166   colvarmodule *cv = this->colvars;
1167   CkLoop_Parallelize(calc_colvars_items_smp, 1, this,
1168                      cv->variables_active_smp()->size(),
1169                      0, cv->variables_active_smp()->size()-1);
1170   return cvm::get_error();
1171 }
1172
1173
1174 void calc_cv_biases_smp(int first, int last, void *result, int paramNum, void *param)
1175 {
1176   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
1177   colvarmodule *cv = proxy->colvars;
1178
1179   cvm::increase_depth();
1180   for (int i = first; i <= last; i++) {
1181     colvarbias *b = (*(cv->biases_active()))[i];
1182     if (cvm::debug()) {
1183       cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
1184                "]: calc_cv_biases_smp(), first = "+cvm::to_str(first)+
1185                ", last = "+cvm::to_str(last)+", bias = "+
1186                b->name+"\n");
1187     }
1188     b->update();
1189   }
1190   cvm::decrease_depth();
1191 }
1192
1193
1194 int colvarproxy_namd::smp_biases_loop()
1195 {
1196   colvarmodule *cv = this->colvars;
1197   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
1198                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1);
1199   return cvm::get_error();
1200 }
1201
1202
1203 void calc_cv_scripted_forces(int paramNum, void *param)
1204 {
1205   colvarproxy_namd *proxy = (colvarproxy_namd *) param;
1206   colvarmodule *cv = proxy->colvars;
1207   if (cvm::debug()) {
1208     cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+
1209              "]: calc_cv_scripted_forces()\n");
1210   }
1211   cv->calc_scripted_forces();
1212 }
1213
1214
1215 int colvarproxy_namd::smp_biases_script_loop()
1216 {
1217   colvarmodule *cv = this->colvars;
1218   CkLoop_Parallelize(calc_cv_biases_smp, 1, this,
1219                      cv->biases_active()->size(), 0, cv->biases_active()->size()-1,
1220                      1, NULL, CKLOOP_NONE,
1221                      calc_cv_scripted_forces, 1, this);
1222   return cvm::get_error();
1223 }
1224
1225 #endif  // #if CMK_SMP && USE_CKLOOP