c8dcfc076619f881a3a217d4e794d52a9838c44e
[namd.git] / colvars / src / colvaratoms.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 <list>
11 #include <vector>
12 #include <algorithm>
13
14 #include "colvarmodule.h"
15 #include "colvarproxy.h"
16 #include "colvarparse.h"
17 #include "colvaratoms.h"
18
19
20 cvm::atom::atom()
21 {
22   index = -1;
23   id = -1;
24   reset_data();
25 }
26
27
28 cvm::atom::atom(int atom_number)
29 {
30   colvarproxy *p = cvm::proxy;
31   index = p->init_atom(atom_number);
32   if (cvm::debug()) {
33     cvm::log("The index of this atom in the colvarproxy arrays is "+
34              cvm::to_str(index)+".\n");
35   }
36   id = p->get_atom_id(index);
37   update_mass();
38   reset_data();
39 }
40
41
42 cvm::atom::atom(cvm::residue_id const &residue,
43                 std::string const     &atom_name,
44                 std::string const     &segment_id)
45 {
46   colvarproxy *p = cvm::proxy;
47   index = p->init_atom(residue, atom_name, segment_id);
48   if (cvm::debug()) {
49     cvm::log("The index of this atom in the colvarproxy_namd arrays is "+
50              cvm::to_str(index)+".\n");
51   }
52   id = p->get_atom_id(index);
53   update_mass();
54   reset_data();
55 }
56
57
58 cvm::atom::atom(atom const &a)
59   : index(a.index)
60 {
61   id = (cvm::proxy)->get_atom_id(index);
62   update_mass();
63   reset_data();
64 }
65
66
67 cvm::atom::~atom()
68 {
69   if (index >= 0) {
70     (cvm::proxy)->clear_atom(index);
71   }
72 }
73
74
75
76 cvm::atom_group::atom_group()
77 {
78   init();
79 }
80
81
82 cvm::atom_group::atom_group(char const *key_in)
83 {
84   key = key_in;
85   init();
86 }
87
88
89 cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms_in)
90 {
91   init();
92   atoms = atoms_in;
93   setup();
94 }
95
96
97 cvm::atom_group::~atom_group()
98 {
99   if (is_enabled(f_ag_scalable) && !b_dummy) {
100     (cvm::proxy)->clear_atom_group(index);
101     index = -1;
102   }
103
104   if (fitting_group) {
105     delete fitting_group;
106     fitting_group = NULL;
107   }
108 }
109
110
111 int cvm::atom_group::add_atom(cvm::atom const &a)
112 {
113   if (a.id < 0) {
114     return COLVARS_ERROR;
115   }
116
117   for (size_t i = 0; i < atoms_ids.size(); i++) {
118     if (atoms_ids[i] == a.id) {
119       if (cvm::debug())
120         cvm::log("Discarding doubly counted atom with number "+
121                  cvm::to_str(a.id+1)+".\n");
122       return COLVARS_OK;
123     }
124   }
125
126   // for consistency with add_atom_id(), we update the list as well
127   atoms_ids.push_back(a.id);
128   atoms.push_back(a);
129   total_mass += a.mass;
130   total_charge += a.charge;
131
132   return COLVARS_OK;
133 }
134
135
136 int cvm::atom_group::add_atom_id(int aid)
137 {
138   if (aid < 0) {
139     return COLVARS_ERROR;
140   }
141
142   for (size_t i = 0; i < atoms_ids.size(); i++) {
143     if (atoms_ids[i] == aid) {
144       if (cvm::debug())
145         cvm::log("Discarding doubly counted atom with number "+
146                  cvm::to_str(aid+1)+".\n");
147       return COLVARS_OK;
148     }
149   }
150
151   atoms_ids.push_back(aid);
152   return COLVARS_OK;
153 }
154
155
156 int cvm::atom_group::remove_atom(cvm::atom_iter ai)
157 {
158   if (is_enabled(f_ag_scalable)) {
159     cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR);
160     return COLVARS_ERROR;
161   }
162
163   if (!this->size()) {
164     cvm::error("Error: trying to remove an atom from an empty group.\n", INPUT_ERROR);
165     return COLVARS_ERROR;
166   } else {
167     total_mass -= ai->mass;
168     total_charge -= ai->charge;
169     atoms_ids.erase(atoms_ids.begin() + (ai - atoms.begin()));
170     atoms.erase(ai);
171   }
172
173   return COLVARS_OK;
174 }
175
176
177 int cvm::atom_group::init()
178 {
179   if (!key.size()) key = "unnamed";
180   description = "atom group " + key;
181   // These may be overwritten by parse(), if a name is provided
182
183   atoms.clear();
184
185   // TODO: check with proxy whether atom forces etc are available
186   init_ag_requires();
187
188   index = -1;
189
190   b_dummy = false;
191   b_center = false;
192   b_rotate = false;
193   b_user_defined_fit = false;
194   fitting_group = NULL;
195
196   noforce = false;
197
198   total_mass = 0.0;
199   total_charge = 0.0;
200
201   cog.reset();
202   com.reset();
203
204   return COLVARS_OK;
205 }
206
207
208 int cvm::atom_group::setup()
209 {
210   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
211     ai->update_mass();
212     ai->update_charge();
213   }
214   update_total_mass();
215   update_total_charge();
216   return COLVARS_OK;
217 }
218
219
220 void cvm::atom_group::update_total_mass()
221 {
222   if (b_dummy) {
223     total_mass = 1.0;
224     return;
225   }
226
227   if (is_enabled(f_ag_scalable)) {
228     total_mass = (cvm::proxy)->get_atom_group_mass(index);
229   } else {
230     total_mass = 0.0;
231     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
232       total_mass += ai->mass;
233     }
234   }
235 }
236
237
238 void cvm::atom_group::reset_mass(std::string &name, int i, int j)
239 {
240   update_total_mass();
241   cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
242            cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
243            " atoms: total mass = "+cvm::to_str(total_mass)+".\n");
244 }
245
246
247 void cvm::atom_group::update_total_charge()
248 {
249   if (b_dummy) {
250     total_charge = 0.0;
251     return;
252   }
253
254   if (is_enabled(f_ag_scalable)) {
255     total_charge = (cvm::proxy)->get_atom_group_charge(index);
256   } else {
257     total_charge = 0.0;
258     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
259       total_charge += ai->charge;
260     }
261   }
262 }
263
264
265 int cvm::atom_group::parse(std::string const &group_conf)
266 {
267   cvm::log("Initializing atom group \""+key+"\".\n");
268
269   // whether or not to include messages in the log
270   // colvarparse::Parse_Mode mode = parse_silent;
271   // {
272   //   bool b_verbose;
273   //   get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
274   //   if (b_verbose) mode = parse_normal;
275   // }
276   // colvarparse::Parse_Mode mode = parse_normal;
277
278   int parse_error = COLVARS_OK;
279
280   // Optional group name will let other groups reuse atom definition
281   if (get_keyval(group_conf, "name", name)) {
282     if ((cvm::atom_group_by_name(this->name) != NULL) &&
283         (cvm::atom_group_by_name(this->name) != this)) {
284       cvm::error("Error: this atom group cannot have the same name, \""+this->name+
285                         "\", as another atom group.\n",
286                 INPUT_ERROR);
287       return INPUT_ERROR;
288     }
289     cvm::main()->register_named_atom_group(this);
290     description = "atom group " + name;
291   }
292
293   // We need to know about fitting to decide whether the group is scalable
294   // and we need to know about scalability before adding atoms
295   bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
296   bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
297   // is the user setting explicit options?
298   b_user_defined_fit = b_defined_center || b_defined_rotate;
299
300   if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) {
301     enable(f_ag_scalable_com);
302     enable(f_ag_scalable);
303   }
304
305   {
306     std::string atoms_of = "";
307     if (get_keyval(group_conf, "atomsOfGroup", atoms_of)) {
308       atom_group * ag = atom_group_by_name(atoms_of);
309       if (ag == NULL) {
310         cvm::error("Error: cannot find atom group with name " + atoms_of + ".\n");
311         return COLVARS_ERROR;
312       }
313       parse_error |= add_atoms_of_group(ag);
314     }
315   }
316
317 //   if (get_keyval(group_conf, "copyOfGroup", source)) {
318 //     // Goal: Initialize this as a full copy
319 //     // for this we'll need an atom_group copy constructor
320 //     return COLVARS_OK;
321 //   }
322
323   {
324     std::string numbers_conf = "";
325     size_t pos = 0;
326     while (key_lookup(group_conf, "atomNumbers", &numbers_conf, &pos)) {
327       parse_error |= add_atom_numbers(numbers_conf);
328       numbers_conf = "";
329     }
330   }
331
332   {
333     std::string index_group_name;
334     if (get_keyval(group_conf, "indexGroup", index_group_name)) {
335       // use an index group from the index file read globally
336       parse_error |= add_index_group(index_group_name);
337     }
338   }
339
340   {
341     std::string range_conf = "";
342     size_t pos = 0;
343     while (key_lookup(group_conf, "atomNumbersRange",
344                       &range_conf, &pos)) {
345       parse_error |= add_atom_numbers_range(range_conf);
346       range_conf = "";
347     }
348   }
349
350   {
351     std::vector<std::string> psf_segids;
352     get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string>());
353     std::vector<std::string>::iterator psii;
354     for (psii = psf_segids.begin(); psii < psf_segids.end(); ++psii) {
355       if ( (psii->size() == 0) || (psii->size() > 4) ) {
356         cvm::error("Error: invalid PSF segment identifier provided, \""+
357                    (*psii)+"\".\n", INPUT_ERROR);
358       }
359     }
360
361     std::string range_conf = "";
362     size_t pos = 0;
363     size_t range_count = 0;
364     psii = psf_segids.begin();
365     while (key_lookup(group_conf, "atomNameResidueRange",
366                       &range_conf, &pos)) {
367       range_count++;
368       if (psf_segids.size() && (range_count > psf_segids.size())) {
369         cvm::error("Error: more instances of \"atomNameResidueRange\" than "
370                    "values of \"psfSegID\".\n", INPUT_ERROR);
371       } else {
372         parse_error |= add_atom_name_residue_range(psf_segids.size() ?
373           *psii : std::string(""), range_conf);
374         if (psf_segids.size()) psii++;
375       }
376       range_conf = "";
377     }
378   }
379
380   {
381     // read the atoms from a file
382     std::string atoms_file_name;
383     if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""))) {
384
385       std::string atoms_col;
386       if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""))) {
387         cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n",
388                    INPUT_ERROR);
389       }
390
391       double atoms_col_value;
392       bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0);
393       if (atoms_col_value_defined && (!atoms_col_value)) {
394         cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
395       }
396
397       // NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
398       cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
399     }
400   }
401
402   // Catch any errors from all the initialization steps above
403   if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
404
405   // checks of doubly-counted atoms have been handled by add_atom() already
406
407   if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos())) {
408     b_dummy = true;
409     // note: atoms_ids.size() is used here in lieu of atoms.size(),
410     // which can be empty for scalable groups
411     if (atoms_ids.size()) {
412       cvm::error("Error: cannot set up group \""+
413                  key+"\" as a dummy atom "
414                  "and provide it with atom definitions.\n", INPUT_ERROR);
415     }
416   } else {
417     b_dummy = false;
418
419     if (!(atoms_ids.size())) {
420       cvm::error("Error: no atoms defined for atom group \""+
421                  key+"\".\n", INPUT_ERROR);
422     }
423
424     // whether these atoms will ever receive forces or not
425     bool enable_forces = true;
426     // disableForces is deprecated
427     if (get_keyval(group_conf, "enableForces", enable_forces, true)) {
428       noforce = !enable_forces;
429     } else {
430       get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
431     }
432   }
433
434   // Now that atoms are defined we can parse the detailed fitting options
435   parse_error |= parse_fitting_options(group_conf);
436
437   if (is_enabled(f_ag_scalable) && !b_dummy) {
438     cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
439     index = (cvm::proxy)->init_atom_group(atoms_ids);
440   }
441
442   bool b_print_atom_ids = false;
443   get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false);
444
445   // Calculate all required properties (such as total mass)
446   setup();
447
448   if (cvm::debug())
449     cvm::log("Done initializing atom group \""+key+"\".\n");
450
451   cvm::log("Atom group \""+key+"\" defined, "+
452             cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+
453             cvm::to_str(total_mass)+", total charge = "+
454             cvm::to_str(total_charge)+".\n");
455
456   if (b_print_atom_ids) {
457     cvm::log("Internal definition of the atom group:\n");
458     cvm::log(print_atom_ids());
459   }
460
461   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
462 }
463
464
465 int cvm::atom_group::add_atoms_of_group(atom_group const * ag)
466 {
467   std::vector<int> const &source_ids = ag->atoms_ids;
468
469   if (source_ids.size()) {
470     atoms_ids.reserve(atoms_ids.size()+source_ids.size());
471
472     if (is_enabled(f_ag_scalable)) {
473       for (size_t i = 0; i < source_ids.size(); i++) {
474         add_atom_id(source_ids[i]);
475       }
476     } else {
477       atoms.reserve(atoms.size()+source_ids.size());
478       for (size_t i = 0; i < source_ids.size(); i++) {
479         // We could use the atom copy constructor, but only if the source
480         // group is not scalable - whereas this works in both cases
481         // atom constructor expects 1-based atom number
482         add_atom(cvm::atom(source_ids[i] + 1));
483       }
484     }
485
486     if (cvm::get_error()) return COLVARS_ERROR;
487   } else {
488     cvm::error("Error: source atom group contains no atoms\".\n", INPUT_ERROR);
489     return COLVARS_ERROR;
490   }
491
492   return COLVARS_OK;
493 }
494
495
496 int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
497 {
498   std::vector<int> atom_indexes;
499
500   if (numbers_conf.size()) {
501     std::istringstream is(numbers_conf);
502     int ia;
503     while (is >> ia) {
504       atom_indexes.push_back(ia);
505     }
506   }
507
508   if (atom_indexes.size()) {
509     atoms_ids.reserve(atoms_ids.size()+atom_indexes.size());
510
511     if (is_enabled(f_ag_scalable)) {
512       for (size_t i = 0; i < atom_indexes.size(); i++) {
513         add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i]));
514       }
515     } else {
516       // if we are handling the group on rank 0, better allocate the vector in one shot
517       atoms.reserve(atoms.size()+atom_indexes.size());
518       for (size_t i = 0; i < atom_indexes.size(); i++) {
519         add_atom(cvm::atom(atom_indexes[i]));
520       }
521     }
522
523     if (cvm::get_error()) return COLVARS_ERROR;
524   } else {
525     cvm::error("Error: no numbers provided for \""
526                "atomNumbers\".\n", INPUT_ERROR);
527     return COLVARS_ERROR;
528   }
529
530   return COLVARS_OK;
531 }
532
533
534 int cvm::atom_group::add_index_group(std::string const &index_group_name)
535 {
536   colvarmodule *cv = cvm::main();
537
538   std::list<std::string>::iterator names_i = cv->index_group_names.begin();
539   std::list<std::vector<int> >::iterator index_groups_i = cv->index_groups.begin();
540   for ( ; names_i != cv->index_group_names.end() ; ++names_i, ++index_groups_i) {
541     if (*names_i == index_group_name)
542       break;
543   }
544
545   if (names_i == cv->index_group_names.end()) {
546     cvm::error("Error: could not find index group "+
547                index_group_name+" among those provided by the index file.\n",
548                INPUT_ERROR);
549     return COLVARS_ERROR;
550   }
551
552   atoms_ids.reserve(atoms_ids.size()+index_groups_i->size());
553
554   if (is_enabled(f_ag_scalable)) {
555     for (size_t i = 0; i < index_groups_i->size(); i++) {
556       add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i]));
557     }
558   } else {
559     atoms.reserve(atoms.size()+index_groups_i->size());
560     for (size_t i = 0; i < index_groups_i->size(); i++) {
561       add_atom(cvm::atom((*index_groups_i)[i]));
562     }
563   }
564
565   if (cvm::get_error())
566     return COLVARS_ERROR;
567
568   return COLVARS_OK;
569 }
570
571
572 int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf)
573 {
574   if (range_conf.size()) {
575     std::istringstream is(range_conf);
576     int initial, final;
577     char dash;
578     if ( (is >> initial) && (initial > 0) &&
579          (is >> dash) && (dash == '-') &&
580          (is >> final) && (final > 0) ) {
581
582       atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
583
584       if (is_enabled(f_ag_scalable)) {
585         for (int anum = initial; anum <= final; anum++) {
586           add_atom_id((cvm::proxy)->check_atom_id(anum));
587         }
588       } else {
589         atoms.reserve(atoms.size() + (final - initial + 1));
590         for (int anum = initial; anum <= final; anum++) {
591           add_atom(cvm::atom(anum));
592         }
593       }
594
595     }
596     if (cvm::get_error()) return COLVARS_ERROR;
597   } else {
598     cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
599                range_conf+"\".\n", INPUT_ERROR);
600     return COLVARS_ERROR;
601   }
602
603   return COLVARS_OK;
604 }
605
606
607 int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
608                                                  std::string const &range_conf)
609 {
610   if (range_conf.size()) {
611     std::istringstream is(range_conf);
612     std::string atom_name;
613     int initial, final;
614     char dash;
615     if ( (is >> atom_name) && (atom_name.size())  &&
616          (is >> initial) && (initial > 0) &&
617          (is >> dash) && (dash == '-') &&
618          (is >> final) && (final > 0) ) {
619
620       atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
621
622       if (is_enabled(f_ag_scalable)) {
623         for (int resid = initial; resid <= final; resid++) {
624           add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid));
625         }
626       } else {
627         atoms.reserve(atoms.size() + (final - initial + 1));
628         for (int resid = initial; resid <= final; resid++) {
629           add_atom(cvm::atom(resid, atom_name, psf_segid));
630         }
631       }
632
633       if (cvm::get_error()) return COLVARS_ERROR;
634     } else {
635       cvm::error("Error: cannot parse definition for \""
636                  "atomNameResidueRange\", \""+
637                  range_conf+"\".\n");
638       return COLVARS_ERROR;
639     }
640   } else {
641     cvm::error("Error: atomNameResidueRange with empty definition.\n");
642     return COLVARS_ERROR;
643   }
644   return COLVARS_OK;
645 }
646
647
648 std::string const cvm::atom_group::print_atom_ids() const
649 {
650   size_t line_count = 0;
651   std::ostringstream os("");
652   for (size_t i = 0; i < atoms_ids.size(); i++) {
653     os << " " << std::setw(9) << atoms_ids[i];
654     if (++line_count == 7) {
655       os << "\n";
656       line_count = 0;
657     }
658   }
659   return os.str();
660 }
661
662
663 int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
664 {
665   if (b_center || b_rotate) {
666
667     if (b_dummy)
668       cvm::error("Error: centerReference or rotateReference "
669                  "cannot be defined for a dummy atom.\n");
670
671     bool b_ref_pos_group = false;
672     std::string fitting_group_conf;
673     if (key_lookup(group_conf, "refPositionsGroup", &fitting_group_conf)) {
674       b_ref_pos_group = true;
675       cvm::log("Warning: keyword \"refPositionsGroup\" is deprecated: please use \"fittingGroup\" instead.\n");
676     }
677
678     if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup", &fitting_group_conf)) {
679       // instead of this group, define another group to compute the fit
680       if (fitting_group) {
681         cvm::error("Error: the atom group \""+
682                    key+"\" has already a reference group "
683                    "for the rototranslational fit, which was communicated by the "
684                    "colvar component.  You should not use fittingGroup "
685                    "in this case.\n", INPUT_ERROR);
686         return INPUT_ERROR;
687       }
688       cvm::log("Within atom group \""+key+"\":\n");
689       fitting_group = new atom_group("fittingGroup");
690       if (fitting_group->parse(fitting_group_conf) == COLVARS_OK) {
691         fitting_group->check_keywords(fitting_group_conf, "fittingGroup");
692         if (cvm::get_error()) {
693           cvm::error("Error setting up atom group \"fittingGroup\".", INPUT_ERROR);
694           return INPUT_ERROR;
695         }
696       }
697     }
698
699     atom_group *group_for_fit = fitting_group ? fitting_group : this;
700
701     get_keyval(group_conf, "refPositions", ref_pos, ref_pos);
702
703     std::string ref_pos_file;
704     if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""))) {
705
706       if (ref_pos.size()) {
707         cvm::error("Error: cannot specify \"refPositionsFile\" and "
708                    "\"refPositions\" at the same time.\n");
709       }
710
711       std::string ref_pos_col;
712       double ref_pos_col_value=0.0;
713
714       if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""))) {
715         // if provided, use PDB column to select coordinates
716         bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0);
717         if (found && ref_pos_col_value == 0.0) {
718           cvm::error("Error: refPositionsColValue, "
719                      "if provided, must be non-zero.\n", INPUT_ERROR);
720           return COLVARS_ERROR;
721         }
722       }
723
724       ref_pos.resize(group_for_fit->size());
725       cvm::load_coords(ref_pos_file.c_str(), &ref_pos, group_for_fit,
726                        ref_pos_col, ref_pos_col_value);
727     }
728
729     if (ref_pos.size()) {
730
731       if (b_rotate) {
732         if (ref_pos.size() != group_for_fit->size())
733           cvm::error("Error: the number of reference positions provided("+
734                      cvm::to_str(ref_pos.size())+
735                      ") does not match the number of atoms within \""+
736                      key+
737                      "\" ("+cvm::to_str(group_for_fit->size())+
738                      "): to perform a rotational fit, "+
739                      "these numbers should be equal.\n", INPUT_ERROR);
740       }
741
742       // save the center of geometry of ref_pos and subtract it
743       center_ref_pos();
744
745     } else {
746       cvm::error("Error: no reference positions provided.\n", INPUT_ERROR);
747       return COLVARS_ERROR;
748     }
749
750     if (b_rotate && !noforce) {
751       cvm::log("Warning: atom group \""+key+
752                "\" will be aligned to a fixed orientation given by the reference positions provided.  "
753                "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
754                "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "
755                "If that happens, use fittingGroup (or a different definition for it if already defined) "
756                "to align the coordinates.\n");
757       // initialize rot member data
758       rot.request_group1_gradients(group_for_fit->size());
759     }
760   }
761
762   // Enable fit gradient calculation only if necessary, and not disabled by the user
763   // This must happen after fitting group is defined so that side-effects are performed
764   // properly (ie. allocating fitting group gradients)
765   {
766     bool b_fit_gradients;
767     get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
768
769     if (b_fit_gradients && (b_center || b_rotate)) {
770       enable(f_ag_fit_gradients);
771     }
772   }
773
774   return COLVARS_OK;
775 }
776
777
778 void cvm::atom_group::do_feature_side_effects(int id)
779 {
780   // If enabled features are changed upstream, the features below should be refreshed
781   switch (id) {
782     case f_ag_fit_gradients:
783       if (b_center || b_rotate) {
784         atom_group *group_for_fit = fitting_group ? fitting_group : this;
785         group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
786         rot.request_group1_gradients(group_for_fit->size());
787       }
788       break;
789   }
790 }
791
792
793 int cvm::atom_group::create_sorted_ids()
794 {
795   // Only do the work if the vector is not yet populated
796   if (sorted_atoms_ids.size())
797     return COLVARS_OK;
798
799   // Sort the internal IDs
800   std::list<int> sorted_atoms_ids_list;
801   for (size_t i = 0; i < this->size(); i++) {
802     sorted_atoms_ids_list.push_back(atoms_ids[i]);
803   }
804   sorted_atoms_ids_list.sort();
805   sorted_atoms_ids_list.unique();
806   if (sorted_atoms_ids_list.size() != this->size()) {
807     return cvm::error("Error: duplicate atom IDs in atom group? (found " +
808                       cvm::to_str(sorted_atoms_ids_list.size()) +
809                       " unique atom IDs instead of " +
810                       cvm::to_str(this->size()) + ").\n", BUG_ERROR);
811   }
812
813   // Compute map between sorted and unsorted elements
814   sorted_atoms_ids.resize(this->size());
815   sorted_atoms_ids_map.resize(this->size());
816   std::list<int>::iterator lsii = sorted_atoms_ids_list.begin();
817   size_t ii = 0;
818   for ( ; ii < this->size(); lsii++, ii++) {
819     sorted_atoms_ids[ii] = *lsii;
820     size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) -
821       atoms_ids.begin();
822     sorted_atoms_ids_map[ii] = pos;
823   }
824
825   return COLVARS_OK;
826 }
827
828
829 int cvm::atom_group::overlap(const atom_group &g1, const atom_group &g2){
830   for (cvm::atom_const_iter ai1 = g1.begin(); ai1 != g1.end(); ai1++) {
831     for (cvm::atom_const_iter ai2 = g2.begin(); ai2 != g2.end(); ai2++) {
832       if (ai1->id == ai2->id) {
833         return (ai1->id + 1); // 1-based index to allow boolean usage
834       }
835     }
836   }
837   return 0;
838 }
839
840
841 void cvm::atom_group::center_ref_pos()
842 {
843   ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
844   std::vector<cvm::atom_pos>::iterator pi;
845   for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
846     ref_pos_cog += *pi;
847   }
848   ref_pos_cog /= (cvm::real) ref_pos.size();
849   for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
850     (*pi) -= ref_pos_cog;
851   }
852 }
853
854
855 void cvm::atom_group::read_positions()
856 {
857   if (b_dummy) return;
858
859   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
860     ai->read_position();
861   }
862
863   if (fitting_group)
864     fitting_group->read_positions();
865 }
866
867
868 int cvm::atom_group::calc_required_properties()
869 {
870   // TODO check if the com is needed?
871   calc_center_of_mass();
872   calc_center_of_geometry();
873
874   if (!is_enabled(f_ag_scalable)) {
875     if (b_center || b_rotate) {
876       if (fitting_group) {
877         fitting_group->calc_center_of_geometry();
878       }
879
880       calc_apply_roto_translation();
881
882       // update COM and COG after fitting
883       calc_center_of_geometry();
884       calc_center_of_mass();
885       if (fitting_group) {
886         fitting_group->calc_center_of_geometry();
887       }
888     }
889   }
890
891   // TODO calculate elements of scalable cvc's here before reduction
892
893   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
894 }
895
896
897 void cvm::atom_group::calc_apply_roto_translation()
898 {
899   // store the laborarory-frame COGs for when they are needed later
900   cog_orig = this->center_of_geometry();
901   if (fitting_group) {
902     fitting_group->cog_orig = fitting_group->center_of_geometry();
903   }
904
905   if (b_center) {
906     // center on the origin first
907     cvm::atom_pos const rpg_cog = fitting_group ?
908       fitting_group->center_of_geometry() : this->center_of_geometry();
909     apply_translation(-1.0 * rpg_cog);
910     if (fitting_group) {
911       fitting_group->apply_translation(-1.0 * rpg_cog);
912     }
913   }
914
915   if (b_rotate) {
916     // rotate the group (around the center of geometry if b_center is
917     // true, around the origin otherwise)
918     rot.calc_optimal_rotation(fitting_group ?
919                               fitting_group->positions() :
920                               this->positions(),
921                               ref_pos);
922
923     cvm::atom_iter ai;
924     for (ai = this->begin(); ai != this->end(); ai++) {
925       ai->pos = rot.rotate(ai->pos);
926     }
927     if (fitting_group) {
928       for (ai = fitting_group->begin(); ai != fitting_group->end(); ai++) {
929         ai->pos = rot.rotate(ai->pos);
930       }
931     }
932   }
933
934   if (b_center) {
935     // align with the center of geometry of ref_pos
936     apply_translation(ref_pos_cog);
937     if (fitting_group) {
938       fitting_group->apply_translation(ref_pos_cog);
939     }
940   }
941   // update of COM and COG is done from the calling routine
942 }
943
944
945 void cvm::atom_group::apply_translation(cvm::rvector const &t)
946 {
947   if (b_dummy) {
948     cvm::error("Error: cannot translate the coordinates of a dummy atom group.\n", INPUT_ERROR);
949     return;
950   }
951
952   if (is_enabled(f_ag_scalable)) {
953     cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR);
954     return;
955   }
956
957   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
958     ai->pos += t;
959   }
960 }
961
962
963 void cvm::atom_group::read_velocities()
964 {
965   if (b_dummy) return;
966
967   if (b_rotate) {
968
969     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
970       ai->read_velocity();
971       ai->vel = rot.rotate(ai->vel);
972     }
973
974   } else {
975
976     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
977       ai->read_velocity();
978     }
979   }
980 }
981
982
983 // TODO make this a calc function
984 void cvm::atom_group::read_total_forces()
985 {
986   if (b_dummy) return;
987
988   if (b_rotate) {
989
990     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
991       ai->read_total_force();
992       ai->total_force = rot.rotate(ai->total_force);
993     }
994
995   } else {
996
997     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
998       ai->read_total_force();
999     }
1000   }
1001 }
1002
1003
1004 int cvm::atom_group::calc_center_of_geometry()
1005 {
1006   if (b_dummy) {
1007     cog = dummy_atom_pos;
1008   } else {
1009     cog.reset();
1010     for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
1011       cog += ai->pos;
1012     }
1013     cog /= this->size();
1014   }
1015   return COLVARS_OK;
1016 }
1017
1018
1019 int cvm::atom_group::calc_center_of_mass()
1020 {
1021   if (b_dummy) {
1022     com = dummy_atom_pos;
1023     if (cvm::debug()) {
1024       cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n");
1025     }
1026   } else if (is_enabled(f_ag_scalable)) {
1027     com = (cvm::proxy)->get_atom_group_com(index);
1028   } else {
1029     com.reset();
1030     for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
1031       com += ai->mass * ai->pos;
1032     }
1033     com /= total_mass;
1034   }
1035   return COLVARS_OK;
1036 }
1037
1038
1039 int cvm::atom_group::calc_dipole(cvm::atom_pos const &com)
1040 {
1041   if (b_dummy) {
1042     cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR);
1043     return COLVARS_ERROR;
1044   }
1045   dip.reset();
1046   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
1047     dip += ai->charge * (ai->pos - com);
1048   }
1049   return COLVARS_OK;
1050 }
1051
1052
1053 void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
1054 {
1055   if (b_dummy) return;
1056
1057   if (is_enabled(f_ag_scalable)) {
1058     scalar_com_gradient = grad;
1059     return;
1060   }
1061
1062   for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
1063     ai->grad = (ai->mass/total_mass) * grad;
1064   }
1065 }
1066
1067
1068 void cvm::atom_group::calc_fit_gradients()
1069 {
1070   if (b_dummy || ! is_enabled(f_ag_fit_gradients)) return;
1071
1072   if (cvm::debug())
1073     cvm::log("Calculating fit gradients.\n");
1074
1075   cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
1076
1077   if (b_center) {
1078     // add the center of geometry contribution to the gradients
1079     cvm::rvector atom_grad;
1080
1081     for (size_t i = 0; i < this->size(); i++) {
1082       atom_grad += atoms[i].grad;
1083     }
1084     if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
1085     atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
1086
1087     for (size_t j = 0; j < group_for_fit->size(); j++) {
1088       group_for_fit->fit_gradients[j] = atom_grad;
1089     }
1090   }
1091
1092   if (b_rotate) {
1093
1094     // add the rotation matrix contribution to the gradients
1095     cvm::rotation const rot_inv = rot.inverse();
1096
1097     for (size_t i = 0; i < this->size(); i++) {
1098
1099       // compute centered, unrotated position
1100       cvm::atom_pos const pos_orig =
1101         rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
1102
1103       // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
1104       cvm::quaternion const dxdq =
1105         rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
1106
1107       for (size_t j = 0; j < group_for_fit->size(); j++) {
1108         // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients
1109         for (size_t iq = 0; iq < 4; iq++) {
1110           group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
1111         }
1112       }
1113     }
1114   }
1115
1116   if (cvm::debug())
1117     cvm::log("Done calculating fit gradients.\n");
1118 }
1119
1120
1121 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
1122 {
1123   if (b_dummy) {
1124     cvm::error("Error: positions are not available "
1125                "from a dummy atom group.\n", INPUT_ERROR);
1126   }
1127
1128   if (is_enabled(f_ag_scalable)) {
1129     cvm::error("Error: atomic positions are not available "
1130                "from a scalable atom group.\n", INPUT_ERROR);
1131   }
1132
1133   std::vector<cvm::atom_pos> x(this->size(), 0.0);
1134   cvm::atom_const_iter ai = this->begin();
1135   std::vector<cvm::atom_pos>::iterator xi = x.begin();
1136   for ( ; ai != this->end(); ++xi, ++ai) {
1137     *xi = ai->pos;
1138   }
1139   return x;
1140 }
1141
1142 std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted(cvm::rvector const &shift) const
1143 {
1144   if (b_dummy) {
1145     cvm::error("Error: positions are not available "
1146                "from a dummy atom group.\n", INPUT_ERROR);
1147   }
1148
1149   if (is_enabled(f_ag_scalable)) {
1150     cvm::error("Error: atomic positions are not available "
1151                "from a scalable atom group.\n", INPUT_ERROR);
1152   }
1153
1154   std::vector<cvm::atom_pos> x(this->size(), 0.0);
1155   cvm::atom_const_iter ai = this->begin();
1156   std::vector<cvm::atom_pos>::iterator xi = x.begin();
1157   for ( ; ai != this->end(); ++xi, ++ai) {
1158     *xi = (ai->pos + shift);
1159   }
1160   return x;
1161 }
1162
1163 std::vector<cvm::rvector> cvm::atom_group::velocities() const
1164 {
1165   if (b_dummy) {
1166     cvm::error("Error: velocities are not available "
1167                "from a dummy atom group.\n", INPUT_ERROR);
1168   }
1169
1170   if (is_enabled(f_ag_scalable)) {
1171     cvm::error("Error: atomic velocities are not available "
1172                "from a scalable atom group.\n", INPUT_ERROR);
1173   }
1174
1175   std::vector<cvm::rvector> v(this->size(), 0.0);
1176   cvm::atom_const_iter ai = this->begin();
1177   std::vector<cvm::atom_pos>::iterator vi = v.begin();
1178   for ( ; ai != this->end(); vi++, ai++) {
1179     *vi = ai->vel;
1180   }
1181   return v;
1182 }
1183
1184 std::vector<cvm::rvector> cvm::atom_group::total_forces() const
1185 {
1186   if (b_dummy) {
1187     cvm::error("Error: total forces are not available "
1188                "from a dummy atom group.\n", INPUT_ERROR);
1189   }
1190
1191   if (is_enabled(f_ag_scalable)) {
1192     cvm::error("Error: atomic total forces are not available "
1193                "from a scalable atom group.\n", INPUT_ERROR);
1194   }
1195
1196   std::vector<cvm::rvector> f(this->size(), 0.0);
1197   cvm::atom_const_iter ai = this->begin();
1198   std::vector<cvm::atom_pos>::iterator fi = f.begin();
1199   for ( ; ai != this->end(); ++fi, ++ai) {
1200     *fi = ai->total_force;
1201   }
1202   return f;
1203 }
1204
1205
1206 // TODO make this an accessor
1207 cvm::rvector cvm::atom_group::total_force() const
1208 {
1209   if (b_dummy) {
1210     cvm::error("Error: total total forces are not available "
1211                "from a dummy atom group.\n", INPUT_ERROR);
1212   }
1213
1214   if (is_enabled(f_ag_scalable)) {
1215     return (cvm::proxy)->get_atom_group_total_force(index);
1216   }
1217
1218   cvm::rvector f(0.0);
1219   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
1220     f += ai->total_force;
1221   }
1222   return f;
1223 }
1224
1225
1226 void cvm::atom_group::apply_colvar_force(cvm::real const &force)
1227 {
1228   if (cvm::debug()) {
1229     log("Communicating a colvar force from atom group to the MD engine.\n");
1230   }
1231
1232   if (b_dummy) return;
1233
1234   if (noforce) {
1235     cvm::error("Error: sending a force to a group that has "
1236                "\"enableForces\" set to off.\n");
1237     return;
1238   }
1239
1240   if (is_enabled(f_ag_scalable)) {
1241     (cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
1242     return;
1243   }
1244
1245   if (b_rotate) {
1246
1247     // rotate forces back to the original frame
1248     cvm::rotation const rot_inv = rot.inverse();
1249     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
1250       ai->apply_force(rot_inv.rotate(force * ai->grad));
1251     }
1252
1253   } else {
1254
1255     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
1256       ai->apply_force(force * ai->grad);
1257     }
1258   }
1259
1260   if ((b_center || b_rotate) && is_enabled(f_ag_fit_gradients)) {
1261
1262     atom_group *group_for_fit = fitting_group ? fitting_group : this;
1263
1264     // Fit gradients are already calculated in "laboratory" frame
1265     for (size_t j = 0; j < group_for_fit->size(); j++) {
1266       (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
1267     }
1268   }
1269 }
1270
1271
1272 void cvm::atom_group::apply_force(cvm::rvector const &force)
1273 {
1274   if (cvm::debug()) {
1275     log("Communicating a colvar force from atom group to the MD engine.\n");
1276   }
1277
1278   if (b_dummy) return;
1279
1280   if (noforce) {
1281     cvm::error("Error: sending a force to a group that has "
1282                "\"enableForces\" set to off.\n");
1283     return;
1284   }
1285
1286   if (is_enabled(f_ag_scalable)) {
1287     (cvm::proxy)->apply_atom_group_force(index, force);
1288     return;
1289   }
1290
1291   if (b_rotate) {
1292
1293     cvm::rotation const rot_inv = rot.inverse();
1294     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
1295       ai->apply_force(rot_inv.rotate((ai->mass/total_mass) * force));
1296     }
1297
1298   } else {
1299
1300     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
1301       ai->apply_force((ai->mass/total_mass) * force);
1302     }
1303   }
1304 }
1305
1306
1307 // Static members
1308
1309 std::vector<colvardeps::feature *> cvm::atom_group::ag_features;
1310
1311