LOOS  v2.3.2
AtomicGroup.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 
24 #include <ios>
25 #include <sstream>
26 #include <iomanip>
27 
28 #include <assert.h>
29 #include <vector>
30 #include <list>
31 #include <algorithm>
32 
33 #include <boost/random.hpp>
34 
35 #include <AtomicGroup.hpp>
36 #include <AtomicNumberDeducer.hpp>
37 #include <Selectors.hpp>
38 
39 #include <boost/unordered_map.hpp>
40 
41 namespace loos {
42 
43 
44  typedef boost::unordered_map<int,int> IMap;
45 
46 
47 
48  const double AtomicGroup::superposition_zero_singular_value = 1e-10;
49 
50 
52  return(new AtomicGroup(*this));
53  }
54 
55 
57  const_iterator i;
58  AtomicGroup res;
59 
60  for (i = atoms.begin(); i != atoms.end(); i++) {
61  pAtom pa(new Atom(**i));
62  res.append(pa);
63  }
64  res._sorted = _sorted;
65  res.box = box.copy();
66 
67  return(res);
68  }
69 
70 
71  // Internal: verify the index into the atom array...
72  int AtomicGroup::rangeCheck(int i) const {
73  if (i < 0)
74  i = atoms.size() + i;
75  if ((unsigned int)i >= atoms.size())
76  throw(std::out_of_range("Bad index for an atom"));
77 
78  return(i);
79  }
80 
81 
82  // Returns the ith atom...
83  pAtom AtomicGroup::getAtom(const int i) const {
84  int j = rangeCheck(i);
85 
86  return(atoms[j]);
87  }
88 
89  // Should these invalidate sort status?
90  pAtom& AtomicGroup::operator[](const int i) {
91  int j = rangeCheck(i);
92  return(atoms[j]);
93  }
94 
95  // For const objects...
96  const pAtom& AtomicGroup::operator[](const int i) const {
97  int j = rangeCheck(i);
98  return(atoms[j]);
99  }
100 
101  // Internal: removes an atom from this group based on the address of the shared pointer...
102  void AtomicGroup::deleteAtom(pAtom pa) {
103  std::vector<pAtom>::iterator iter;
104 
105  iter = find(atoms.begin(), atoms.end(), pa);
106  if (iter == atoms.end())
107  throw(LOOSError(*pa, "Attempting to delete an atom that is not in the passed AtomicGroup"));
108 
109  atoms.erase(iter);
110  _sorted = false;
111  }
112 
113 
114  // Append each atom from the passed vector onto this group...
115  AtomicGroup& AtomicGroup::append(std::vector<pAtom> pas) {
116  std::vector<pAtom>::iterator i;
117 
118  for (i=pas.begin(); i != pas.end(); i++)
119  atoms.push_back(*i);
120 
121  _sorted = false;
122  return(*this);
123  }
124 
125 
126  // Append all atoms from the passed group onto this one
128  std::vector<pAtom>::const_iterator i;
129 
130  for (i=grp.atoms.begin(); i != grp.atoms.end(); i++)
131  addAtom(*i);
132 
133  _sorted = false;
134  return(*this);
135  }
136 
137 
138  // Remove all atoms in the passed vector
139  AtomicGroup& AtomicGroup::remove(std::vector<pAtom> pas) {
140  std::vector<pAtom>::iterator i;
141 
142  for (i=pas.begin(); i != pas.end(); i++)
143  deleteAtom(*i);
144 
145  _sorted = false;
146  return(*this);
147  }
148 
149 
150  // Removes all atoms contained in the passed group from this one...
152 
153 
154  if (&grp == this)
155  atoms.clear(); // Assume caller meant to clean out AtomicGroup
156  else {
157  std::vector<pAtom>::const_iterator i;
158 
159  for (i=grp.atoms.begin(); i != grp.atoms.end(); i++)
160  deleteAtom(*i);
161 
162  _sorted = false;
163  return(*this);
164  }
165 
166  return(*this);
167  }
168 
169 
170  // Concatenation operations...
171  AtomicGroup& AtomicGroup::operator+=(const AtomicGroup& rhs) {
172  append(rhs);
173  return(*this);
174  }
175 
176  AtomicGroup& AtomicGroup::operator+=(const pAtom& rhs) {
177  atoms.push_back(rhs);
178  _sorted = false;
179  return(*this);
180  }
181 
182  AtomicGroup AtomicGroup::operator+(const AtomicGroup& rhs) {
183  AtomicGroup res(*this);
184  res += rhs;
185  return(res);
186  }
187 
188  AtomicGroup AtomicGroup::operator+(const pAtom& rhs) {
189  AtomicGroup res(*this);
190  res += rhs;
191  return(res);
192  }
193 
194 
195  bool AtomicGroup::allHaveProperty(const Atom::bits& property) const {
196  for (const_iterator ci = atoms.begin(); ci != atoms.end(); ++ci)
197  if (!(*ci)->checkProperty(property))
198  return(false);
199 
200  return(true);
201  }
202 
203  bool AtomicGroup::anyHaveProperty(const Atom::bits& property) const {
204  for (const_iterator ci = atoms.begin(); ci != atoms.end(); ++ci)
205  if ((*ci)->checkProperty(property))
206  return(true);
207 
208  return(false);
209  }
210 
211  bool AtomicGroup::hasBonds(void) const {
212  const_iterator ci;
213 
214  for (ci = atoms.begin(); ci != atoms.end(); ++ci)
215  if ((*ci)->checkProperty(Atom::bondsbit))
216  return(true);
217 
218  return(false);
219  }
220 
221 
223  const_iterator ci;
224 
225  for (ci = atoms.begin(); ci != atoms.end(); ++ci)
226  (*ci)->clearBonds();
227  }
228 
229 
230  bool AtomicGroup::hasCoords(void) const {
231  const_iterator ci;
232 
233  for (ci = atoms.begin(); ci != atoms.end(); ++ci)
234  if (!((*ci)->checkProperty(Atom::coordsbit)))
235  return(false);
236 
237  return(true);
238  }
239 
240 
241  void AtomicGroup::sort(void) {
242  CmpById comp;
243 
244  if (! _sorted)
245  std::sort(atoms.begin(), atoms.end(), comp);
246 
247  _sorted = true;
248  }
249 
250 
251 
252  // Internal: calculates the start and stop iterators given offset and len args
253  // as in PERL's substr()...
254 
255  boost::tuple<AtomicGroup::iterator, AtomicGroup::iterator> AtomicGroup::calcSubsetIterators(const int offset, const int len) {
256  unsigned int a, b;
257 
258  if (offset < 0) {
259  b = atoms.size() + offset + 1;
260  a = (len == 0) ? 0 : b - len;
261  } else {
262  a = offset;
263  b = (len == 0) ? atoms.size() : a + len;
264  }
265 
266  if (b-a >= atoms.size())
267  throw(std::range_error("Indices out of bounds for subsetting"));
268 
269  boost::tuple<iterator, iterator> res(atoms.begin() + a, atoms.begin() + b);
270 
271  return(res);
272  }
273 
274 
275 
276  AtomicGroup AtomicGroup::subset(const int offset, const int len) {
277  AtomicGroup res;
278 
279  boost::tuple<iterator, iterator> iters = calcSubsetIterators(offset, len);
280  res.atoms.insert(res.atoms.begin(), boost::get<0>(iters), boost::get<1>(iters));
281 
282  res.box = box;
283  return(res);
284  }
285 
286 
287  AtomicGroup AtomicGroup::excise(const int offset, const int len) {
288  AtomicGroup res;
289 
290  boost::tuple<iterator, iterator> iters = calcSubsetIterators(offset, len);
291 
292  res.atoms.insert(res.atoms.begin(), boost::get<0>(iters), boost::get<1>(iters));
293  atoms.erase(boost::get<0>(iters), boost::get<1>(iters));
294 
295  _sorted = false;
296 
297  res.box = box;
298  return(res);
299  }
300 
301  // Select atoms from the current group, adding them to a new group
302  // based on the sel functor/predicate...
303 
305  AtomicGroup res;
306 
307  std::vector<pAtom>::const_iterator i;
308  for (i=atoms.begin(); i != atoms.end(); i++)
309  if (sel(*i))
310  res.addAtom(*i);
311 
312  res.box = box;
313  return(res);
314  }
315 
316 
317  // Split up a group into a vector of groups based on unique segids...
318  std::vector<AtomicGroup> AtomicGroup::splitByUniqueSegid(void) const {
319  const_iterator i;
320  std::list<std::string> segids;
321 
322  for (i = atoms.begin(); i != atoms.end(); i++)
323  if (find(segids.begin(), segids.end(), (*i)->segid()) == segids.end())
324  segids.push_front((*i)->segid());
325 
326  // Reverse the order so that the chunks are in the same order they appear
327  // in the file...
328  uint n = segids.size();
329  std::vector<AtomicGroup> results(n);
330  uint j = n-1;
331  for (std::list<std::string>::const_iterator i = segids.begin(); i != segids.end(); ++i) {
332  SegidSelector segid(*i);
333  results[j--] = select(segid);
334  }
335 
336  std::vector<AtomicGroup>::iterator g;
337  for (g=results.begin(); g!=results.end(); g++) {
338  g->box = box;
339  }
340 
341  return(results);
342  }
343 
344  std::map<std::string, AtomicGroup> AtomicGroup::splitByName(void) const {
345  const_iterator i;
346  std::map<std::string, AtomicGroup> groups;
347 
348  // Loop over atoms, adding them to groups based on their name, creating the
349  // map entry for each new name as we find it
350  std::map<std::string, AtomicGroup>::iterator g;
351  for (i = atoms.begin(); i != atoms.end(); ++i) {
352  g = groups.find((*i)->name());
353  if (g == groups.end()) { // not found, need to create a new AG
354  AtomicGroup ag;
355  ag.append(*i);
356  ag.box = box; // copy the current groups periodic box
357  groups[(*i)->name()] = ag;
358  } else { // found group for that atom name,
359  // so add the atom to it
360  g->second.append(*i);
361  }
362 
363  }
364 
365  return(groups);
366  }
367 
388  std::vector<AtomicGroup> AtomicGroup::sortingSplitByMolecule(void) {
389  std::vector<AtomicGroup> molecules;
390 
391  // If no connectivity, just return the entire group...
392  if (!hasBonds()) {
393  sort();
394  molecules.push_back(*this);
395  } else {
396  HashInt seen; // Track what atoms we've already processed
397  AtomicGroup current; // The molecule we're currently building...
398  AtomicGroup working(*this); // Copy, so we can sort without mucking up original order
399  working.sort();
400 
401  int n = size();
402  for (int i=0; i<n; i++) {
403  HashInt::iterator it = seen.find(working.atoms[i]->id());
404  if (it != seen.end())
405  continue;
406 
407  walkBonds(current, seen, working, working.atoms[i]);
408  if (current.size() != 0) { // Just in case...
409  current.sort();
410  molecules.push_back(current);
411  current = AtomicGroup();
412  }
413 
414  }
415  }
416 
417  // copy the box over
418  std::vector<AtomicGroup>::iterator m;
419  for (m=molecules.begin(); m!=molecules.end(); m++) {
420  m->box = box;
421  }
422  return(molecules);
423  }
424 
425 
426  void AtomicGroup::walkBonds(AtomicGroup& current, HashInt& seen, AtomicGroup& working, pAtom& moi) {
427  int myid = moi->id();
428  HashInt::iterator it = seen.find(myid);
429 
430  // If we've touched this atom before, stop recursing and return.
431  if (it != seen.end())
432  return;
433 
434  // This is a novel atom, so append it into the group we're currently building.
435  seen.insert(myid);
436  current.addAtom(moi);
437 
438  // Just in case it's a solo-atom... This probably should indicate
439  // some kind of error...?
440  if (!(moi->hasBonds()))
441  return;
442 
443  // Now find atoms that are bound to the current atom and recurse
444  // through them...
445 
446  std::vector<int> bonds = moi->getBonds();
447  std::vector<int>::const_iterator citer;
448  for (citer = bonds.begin(); citer != bonds.end(); citer++) {
449  pAtom toi = working.findById(*citer);
450  if (toi != 0)
451  walkBonds(current, seen, working, toi);
452  }
453  }
454 
455 
461  std::vector<AtomicGroup> AtomicGroup::splitByResidue(void) const {
462  std::vector<AtomicGroup> residues;
463 
464  int curr_resid = atoms[0]->resid();
465  std::string curr_segid = atoms[0]->segid();
466 
467  AtomicGroup residue;
468  AtomicGroup::const_iterator ci;
469  for (ci = atoms.begin(); ci != atoms.end(); ++ci) {
470  if (curr_resid != (*ci)->resid() || (*ci)->segid() != curr_segid) {
471  residues.push_back(residue);
472  residue = AtomicGroup();
473  curr_resid = (*ci)->resid();
474  curr_segid = (*ci)->segid();
475  }
476  residue.append(*ci);
477  }
478 
479  if (residue.size() != 0)
480  residues.push_back(residue);
481 
482 
483  // Copy the box information
484  for (std::vector<AtomicGroup>::iterator i = residues.begin(); i != residues.end(); ++i)
485  i->box = box;
486 
487  return(residues);
488  }
489 
490 
491  // Find an atom based on atomid
492  // Returns 0 (null shared_ptr) if not found...
493  pAtom AtomicGroup::findById_binarySearch(const int id) const {
494  int bottom = 0, top = size()-1, middle;
495 
496  while (top > bottom) {
497  middle = bottom + (top - bottom) / 2;
498  if (atoms[middle]->id() < id)
499  bottom = middle + 1;
500  else
501  top = middle;
502  }
503 
504  if (atoms[bottom]->id() == id)
505  return(atoms[bottom]);
506 
507  return(pAtom());
508  }
509 
510 
511  pAtom AtomicGroup::findById_linearSearch(const int id) const {
512  for (AtomicGroup::const_iterator i = begin(); i != end(); ++i)
513  if ((*i)->id() == id)
514  return(*i);
515 
516  return(pAtom());
517  }
518 
519 
520  pAtom AtomicGroup::findById(const int id) const {
521  if (sorted())
522  return(findById_binarySearch(id));
523 
524  return(findById_linearSearch(id));
525  }
526 
527 
539  AtomicGroup AtomicGroup::groupFromID(const std::vector<int> &id_list) const {
540  AtomicGroup result;
541 
542  result.box = box;
543 
544  for (unsigned int i=0; i<id_list.size(); i++) {
545  pAtom pa = findById(id_list[i]);
546  if (pa)
547  result.addAtom(pa);
548  }
549  return(result);
550  }
551 
552  // Get all atoms associated with the residue that contains the
553  // passed atom... The returned atoms will not be in order. If
554  // you want that, then explicitly sort the group.
556  iterator i;
557  AtomicGroup result;
558 
559  result.box = box;
560  i = find(atoms.begin(), atoms.end(), res);
561  if (i == atoms.end())
562  return(result);
563 
564  iterator j = i;
565 
566  while (j >= atoms.begin()) {
567  if ((*j)->resid() == res->resid() && (*j)->segid() == res->segid())
568  result.addAtom(*j);
569  else
570  break;
571  --j;
572  }
573 
574  j = i+1;
575  while (j < atoms.end()) {
576  if ((*j)->resid() == res->resid() && (*j)->segid() == res->segid())
577  result.addAtom(*j);
578  else
579  break;
580 
581  ++j;
582  }
583 
584  return(result);
585  }
586 
587 
588  namespace {
589  // renumber the atomids of the group...
590  void renumberWithoutBonds(AtomicGroup& grp, const int start, const int stride) {
591  AtomicGroup::iterator i;
592  int id = start;
593 
594  for (i=grp.begin(); i != grp.end(); i++, id += stride)
595  (*i)->id(id);
596  }
597 
598 
599  void renumberWithBonds(AtomicGroup& grp, const int start, const int stride) {
600  IMap bond_map;
601 
602  int id = start;
603  for (AtomicGroup::iterator i = grp.begin(); i != grp.end(); ++i, id += stride) {
604  bond_map[(*i)->id()] = id;
605  (*i)->id(id);
606  }
607 
608  // Now transform bond list
609  for (AtomicGroup::iterator i = grp.begin(); i != grp.end(); ++i, id += stride) {
610  if ((*i)->hasBonds()) {
611  std::vector<int> bonds = (*i)->getBonds();
612  for (std::vector<int>::iterator j = bonds.begin(); j != bonds.end(); ++j) {
613  IMap::iterator k = bond_map.find(*j);
614  if (k != bond_map.end())
615  *j = k->second;
616  }
617  (*i)->setBonds(bonds);
618  }
619  }
620  }
621  };
622 
623  void AtomicGroup::renumber(const int start, const int stride) {
624 
625  if (hasBonds())
626  renumberWithBonds(*this, start, stride);
627  else
628  renumberWithoutBonds(*this, start, stride);
629  }
630 
631  // Get the min and max atomid's...
632  int AtomicGroup::minId(void) const {
633  const_iterator i;
634 
635  if (atoms.size() == 0)
636  return(-1);
637 
638  int min = atoms[0]->id();
639  for (i = atoms.begin()+1; i != atoms.end(); i++)
640  if ((*i)->id() < min)
641  min = (*i)->id();
642 
643  return(min);
644  }
645 
646 
647  int AtomicGroup::maxId(void) const {
648  const_iterator i;
649 
650  if (atoms.size() == 0)
651  return(-1);
652  int max = atoms[0]->id();
653  for (i=atoms.begin()+1; i != atoms.end(); i++)
654  if ((*i)->id() > max)
655  max = (*i)->id();
656 
657  return(max);
658  }
659 
660 
661  int AtomicGroup::minResid(void) const {
662  const_iterator i;
663 
664  if (atoms.size() == 0)
665  return(-1);
666 
667  int min = atoms[0]->resid();
668  for (i = atoms.begin()+1; i != atoms.end(); i++)
669  if ((*i)->resid() < min)
670  min = (*i)->resid();
671 
672  return(min);
673  }
674 
675  int AtomicGroup::maxResid(void) const {
676  const_iterator i;
677 
678  if (atoms.size() == 0)
679  return(-1);
680 
681  int max = atoms[0]->resid();
682  for (i = atoms.begin()+1; i != atoms.end(); i++)
683  if ((*i)->resid() > max)
684  max = (*i)->resid();
685 
686  return(max);
687  }
688 
689 
690  // Count the number of higher structural elements...
691  int AtomicGroup::numberOfResidues(void) const {
692 
693  if (atoms.size() == 0)
694  return(0);
695 
696  const_iterator i;
697  int n = 1;
698  int curr_resid = atoms[0]->resid();
699  std::string curr_segid = atoms[0]->segid();
700 
701  for (i=atoms.begin()+1; i !=atoms.end(); i++)
702  if (((*i)->resid() != curr_resid) ||
703  ((*i)->segid() != curr_segid)) {
704  ++n;
705  curr_resid = (*i)->resid();
706  curr_segid = (*i)->segid();
707  }
708 
709  return(n);
710  }
711 
712 
713  int AtomicGroup::numberOfSegids(void) const {
714 
715  if (atoms.size() == 0)
716  return(0);
717 
718  const_iterator i;
719  int n = 1;
720  std::string curr_segid = atoms[0]->segid();
721 
722  for (i=atoms.begin()+1; i !=atoms.end(); i++)
723  if ((*i)->segid() != curr_segid) {
724  ++n;
725  curr_segid = (*i)->segid();
726  }
727 
728  return(n);
729  }
730 
731 
732 
733  AtomicGroup operator+(const pAtom& lhs, const pAtom& rhs) {
734  AtomicGroup res;
735  res.append(lhs);
736  res.append(rhs);
737  return(res);
738  }
739 
740 
741  AtomicGroup operator+(const pAtom& lhs, const AtomicGroup& rhs) {
742  AtomicGroup res(rhs);
743  res += lhs;
744  return(res);
745  }
746 
747 
748 
750 
751  if (size() != rhs.size())
752  return(false);
753 
754  if (this == &rhs)
755  return(true);
756 
757  sort();
758  rhs.sort();
759 
760  int n = size();
761  for (int i = 0; i < n; i++)
762  if (atoms[i] != rhs.atoms[i])
763  return(false);
764 
765  return(true);
766  }
767 
768 
769  bool AtomicGroup::operator==(const AtomicGroup& rhs) const {
770  if (size() != rhs.size())
771  return(false);
772 
773  if (this == &rhs)
774  return(true);
775 
776  const std::vector<pAtom> *lp;
777  const std::vector<pAtom> *rp;
778  std::vector<pAtom> lhs_atoms, rhs_atoms;
779  CmpById comp;
780  if (!sorted()) {
781  lhs_atoms = atoms;
782  std::sort(lhs_atoms.begin(), lhs_atoms.end(), comp);
783  lp = &lhs_atoms;
784  } else
785  lp = &atoms;
786 
787  if (!rhs.sorted()) {
788  rhs_atoms = rhs.atoms;
789  std::sort(rhs_atoms.begin(), rhs_atoms.end(), comp);
790  rp = &rhs_atoms;
791  } else
792  rp = &rhs.atoms;
793 
794  int n = size();
795  for (int i = 0; i<n; i++)
796  if ((*lp)[i] != (*rp)[i])
797  return(false);
798 
799  return(true);
800  }
801 
802 
804  if (!(isPeriodic()))
805  throw(LOOSError("trying to reimage a non-periodic group"));
806  GCoord com = centroid();
807  GCoord reimaged = com;
808  reimaged.reimage(periodicBox());
809  GCoord trans = reimaged - com;
810  const_iterator a;
811  for (a=atoms.begin(); a!=atoms.end(); a++) {
812  (*a)->coords() += trans;
813  }
814  }
815 
817  if (!(isPeriodic()))
818  throw(LOOSError("trying to reimage a non-periodic group"));
819  const_iterator a;
820  GCoord box = periodicBox();
821  for (a=atoms.begin(); a!=atoms.end(); a++) {
822  (*a)->coords().reimage(box);
823  }
824  }
825 
835  void AtomicGroup::mergeImage(pAtom &p ) {
836  GCoord ref = p->coords();
837 
838  translate(-ref);
839  reimageByAtom();
840  translate(ref);
841  }
842 
847  mergeImage(atoms[0]);
848  }
849 
850 
851 
852 
853 
854  void AtomicGroup::findBonds(const double dist) {
855  AtomicGroup::iterator ij;
856  double dist2 = dist * dist;
857 
858  for (ij = begin(); ij != end() - 1; ++ij) {
859  AtomicGroup::iterator ii;
860  GCoord u = (*ij)->coords();
861 
862  for (ii = ij + 1; ii != end(); ++ii) {
863  if (u.distance2((*ii)->coords()) < dist2) {
864  (*ij)->addBond(*ii);
865  (*ii)->addBond(*ij);
866  }
867  }
868  }
869  }
870 
871 
872 
887 
888  for (AtomicGroup::iterator j = begin(); j != end(); ++j)
889  if ((*j)->hasBonds()) {
890  std::vector<int> bonds = (*j)->getBonds();
891  std::vector<int> pruned_bonds;
892  for (std::vector<int>::const_iterator i = bonds.begin(); i != bonds.end(); ++i)
893  if (findById(*i) != 0)
894  pruned_bonds.push_back(*i);
895  (*j)->setBonds(pruned_bonds);
896  }
897  }
898 
899 
910  {
911  for (uint i=0; i<size(); ++i)
912  atoms[i]->index(i);
913  }
914 
915 
924  uint n = 0;
925 
926  for (AtomicGroup::iterator i = begin(); i != end(); ++i)
927  if ((*i)->checkProperty(Atom::massbit)) {
928  uint an = loos::deduceAtomicNumberFromMass((*i)->mass(), tol);
929  if (an) {
930  (*i)->atomic_number(an);
931  ++n;
932  }
933  }
934 
935  return(n);
936  }
937 
938 
939  void AtomicGroup::copyCoordinatesFrom(const AtomicGroup& g, const uint offset, const uint length) {
940  uint n = (length == 0 || length > g.size()) ? g.size() : length;
941 
942 
943  for (uint i=0; i<n && i+offset<atoms.size(); ++i)
944  atoms[i+offset]->coords(g[i]->coords());
945  }
946 
947 
948  std::vector<uint> AtomicGroup::atomOrderMapFrom(const AtomicGroup& g) {
949  if (g.size() != size())
950  throw(LOOSError("Cannot map atom order between groups of different sizes"));
951 
952  std::vector<AtomicGroup> other = g.splitByResidue();
953  std::vector<AtomicGroup> self = splitByResidue();
954 
955  std::vector<uint> order;
956  uint idx = 0;
957  for (uint k=0; k<self.size(); ++k) {
958  if (self[k][0]->resname() != other[k][0]->resname())
959  throw(LOOSError("Mismatched residues while trying to map atom order between groups"));
960  for (uint j=0; j<self[k].size(); ++j) {
961  uint i;
962  for (i=0; i<other[k].size(); ++i)
963  if (self[k][j]->name() == other[k][i]->name())
964  break;
965  if (i < other[k].size())
966  order.push_back(idx + i);
967  else
968  throw(LOOSError(*(self[k][j]), "Cannot find match while constructing atom order map"));
969  }
970  idx += other[k].size();
971  }
972 
973  return(order);
974  }
975 
976  void AtomicGroup::copyMappedCoordinatesFrom(const AtomicGroup& g, const std::vector<uint>& map) {
977  if (g.size() != map.size())
978  throw(LOOSError("Atom order map is of incorrect size to copy coordinates"));
979  if (g.size() != size())
980  throw(LOOSError("Cannot copy coordinates (with atom ordering) from an AtomicGroup of a different size"));
981 
982  for (uint i=0; i<map.size(); ++i)
983  atoms[i]->coords(g[map[i]]->coords());
984  }
985 
987  std::vector<uint> map = atomOrderMapFrom(g);
989  }
990 
991 
992  // XMLish output...
993  std::ostream& operator<<(std::ostream& os, const AtomicGroup& grp) {
994  AtomicGroup::const_iterator i;
995  if (grp.isPeriodic())
996  os << "<GROUP PERIODIC='" << grp.box.box() << "'>\n";
997  else
998  os << "<GROUP>\n";
999  for (i=grp.atoms.begin(); i != grp.atoms.end(); i++)
1000  os << " " << **i << std::endl;
1001  os << "</GROUP>";
1002 
1003  return(os);
1004  }
1005 
1006 
1007  void AtomicGroup::setGroupConnectivity() {
1008  for (AtomicGroup::iterator i = atoms.begin(); i != atoms.end(); ++i)
1009  (*i)->setProperty(Atom::bondsbit);
1010  }
1011 
1012 
1013 
1014  void AtomicGroup::setCoords(double* seq, int m, int n) {
1015  if (n != 3 || static_cast<uint>(m) != size())
1016  throw(LOOSError("Invalid dimensions in AtomicGroup::setCoords()"));
1017 
1018  for (int j=0; j<m; ++j)
1019  for (int i=0; i<n; ++i)
1020  atoms[j]->coords()[i] = seq[j*n+i];
1021  }
1022 
1023 
1024  void AtomicGroup::getCoords(double** outseq, int* m, int* n) {
1025  double* dp = static_cast<double*>(malloc(size() * 3 * sizeof(double)));
1026  for (uint j=0; j<size(); ++j)
1027  for (uint i=0; i<3; ++i)
1028  dp[j*3+i] = atoms[j]->coords()[i];
1029 
1030  *m = size();
1031  *n = 3;
1032 
1033  *outseq = dp;
1034  }
1035 
1037  std::vector<AtomicGroup> mols = splitByMolecule();
1038  AtomicGroup centers;
1039  for (std::vector<AtomicGroup>::const_iterator i = mols.begin(); i != mols.end(); ++i) {
1040  pAtom orig = (*i)[0];
1041  pAtom atom(new Atom(*orig));
1042  atom->name("CEN");
1043  atom->coords((*i).centerOfMass());
1044  centers.append(atom);
1045  }
1046  centers.box = box;
1047  return(centers);
1048  }
1049 
1051  std::vector<AtomicGroup> residues = splitByResidue();
1052  AtomicGroup centers;
1053  for (std::vector<AtomicGroup>::const_iterator i = residues.begin(); i != residues.end(); ++i) {
1054  pAtom orig = (*i)[0];
1055  pAtom atom(new Atom(*orig));
1056  atom->name("CEN");
1057  atom->coords((*i).centerOfMass());
1058  centers.append(atom);
1059  }
1060  centers.box = box;
1061  return(centers);
1062  }
1063 
1064 
1065 
1066 }
AtomicGroup subset(const int offset, const int len=0)
subset() and excise() args are patterned after perl's substr...
pAtom getAtom(const int i) const
Get the ith atom from this group.
Definition: AtomicGroup.cpp:83
bool allHaveProperty(const Atom::bits &property) const
True if all atoms in the group have the passed property(ies)
AtomicGroup select(const AtomSelector &sel) const
Return a group consisting of atoms for which sel predicate returns true...
void reimageByAtom()
Reimage atoms individually into the primary cell.
bool operator==(AtomicGroup &rhs)
Equality test for two groups.
void findBonds(const double dist=1.65)
Distance-based search for bonds.
std::vector< AtomicGroup > splitByResidue(void) const
Returns a vector of AtomicGroups, each comprising a single residue.
AtomicGroup & append(pAtom pa)
Append the atom onto the group.
std::map< std::string, AtomicGroup > splitByName(void) const
Returns a vector of AtomicGroups, each containing atoms with the same name.
bool hasBonds(void) const
Does any atom in the group have bond information???
pAtom findById(const int id) const
Find a contained atom by its atomid.
AtomicGroup groupFromID(const std::vector< int > &id_list) const
Create a new group from a vector of atomids.
double distance2(const Coord< T > &o) const
Distance squared between two coordinates.
Definition: Coord.hpp:429
bool isPeriodic(void) const
Test whether or not periodic boundary conditions are set.
void renumber(const int start=1, const int stride=1)
Renumber the atomid's of the contained atoms...
AtomicGroup getResidue(pAtom res)
void copyCoordinatesFrom(const AtomicGroup &g, const uint offset=0, const uint length=0)
Copy coordinates from g into current group.
Basic Atom class for handling atom properties.
Definition: Atom.hpp:50
std::vector< AtomicGroup > splitByMolecule(void) const
Returns a vector of AtomicGroups split based on bond connectivity.
AtomicGroup centrifyByResidue() const
Replace a group with the cente of masses of contained residues (see centrifyByMolecule()) ...
Virtual base-class for selecting atoms from a group.
Definition: AtomicGroup.hpp:56
void translate(const GCoord &v)
Translate an atomic group by vector v.
void getCoords(double **outseq, int *m, int *n)
bits
Bits in the bitmask that flag what properties have actually been set.
Definition: Atom.hpp:65
Generic LOOS exception.
Definition: exceptions.hpp:40
bool sorted(void) const
Is the array of atoms already sorted???
AtomicGroup & remove(pAtom pa)
Delete a single atom.
void sort(void)
Sort based on atomid.
pAtom & operator[](const int i)
Same as getAtom(i)
Definition: AtomicGroup.cpp:90
void copyMappedCoordinatesFrom(const AtomicGroup &g, const std::vector< uint > &order)
Given a mapping of atom order, copy the coordinates into the current group.
Predicate for selecting atoms based on the passed segid string.
Definition: Selectors.hpp:66
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
uint deduceAtomicNumberFromMass(const double tol=0.1)
Deduce atomic number from mass (if present), returning number of atoms assigned.
std::vector< uint > atomOrderMapFrom(const AtomicGroup &g)
Map the order of atoms in AtomicGroup g into the current group.
void clearBonds(void)
Remove any bonding information present in contained atoms.
AtomicGroup excise(const int offset, const int len=0)
excise returns the excised atoms as a group...
virtual AtomicGroup * clone(void) const
Creates a lightweight clone of this group (for polymorphism)
Definition: AtomicGroup.cpp:51
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
void pruneBonds()
Attempt to prune connectivity (only retain bonds to atoms within this AtomicGroup) ...
unsigned int deduceAtomicNumberFromMass(const double mass, const double tolerance)
Deduce an atomic number from the mass.
void mergeImage()
Takes a group that's split across a periodic boundary and reimages it so it's all together...
Namespace for most things not already encapsulated within a class.
bool anyHaveProperty(const Atom::bits &property) const
True if any atom in the group have the passed property(ies)
bool hasCoords(void) const
Do all the atoms in the group have coordinates?
void resetAtomIndices()
Reset the atom indices (used for interfacing with trajectories)
std::vector< AtomicGroup > splitByUniqueSegid(void) const
Returns a vector of AtomicGroups split from the current group based on segid.
void setCoords(double *seq, int m, int n)
void reimage(const Coord< T > &box)
Handle coordinates with periodic boundary conditions.
Definition: Coord.hpp:404
AtomicGroup centrifyByMolecule() const
Replace a group with the center of masses of contained molecules.
GCoord centroid(void) const
Centroid of atoms (ignores mass, operates in group coordinates)
GCoord periodicBox(void) const
Fetch the periodic boundary conditions.