IMP  2.0.0
The Integrative Modeling Platform
charmm_topology.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/charmm_topology.h
3  * \brief Classes for handling CHARMM-style topology.
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPATOM_CHARMM_TOPOLOGY_H
10 #define IMPATOM_CHARMM_TOPOLOGY_H
11 
12 #include "IMP/Object.h"
13 #include "Hierarchy.h"
14 #include "Atom.h"
15 #include <IMP/atom/atom_config.h>
16 
17 #include <string>
18 #include <vector>
19 
20 IMPATOM_BEGIN_NAMESPACE
21 
22 //! A single atom in a CHARMM topology.
23 /** Each atom has a name (unique to the residue), a CHARMM type (used to
24  look up parameters such as radii and bond lengths in the parameter file)
25  and an electrostatic charge.
26  \see CHARMMAtom
27  */
29  std::string name_;
30  std::string charmm_type_;
31  double charge_;
32 public:
33  CHARMMAtomTopology(std::string name) : name_(name) {};
34 
35  CHARMMAtomTopology(std::string name, const CHARMMAtomTopology &other)
36  : name_(name), charmm_type_(other.charmm_type_),
37  charge_(other.charge_) {};
38 
39  std::string get_name() const { return name_; }
40  std::string get_charmm_type() const { return charmm_type_; }
41  double get_charge() const { return charge_; }
42  void set_charmm_type(std::string charmm_type) { charmm_type_ = charmm_type; }
43  void set_charge(double charge) { charge_ = charge; }
44  IMP_SHOWABLE_INLINE(CHARMMAtomTopology, {out << "name: " << name_
45  << "; CHARMM type: " << charmm_type_
46  << "; charge: " << charge_;});
47 };
48 
50 
52 
53 //! The end of a bond, angle, dihedral, improper, or internal coordinate.
54 /** An endpoint is an atom in a residue; bonds have two endpoints,
55  angles three, and dihedrals, impropers and internal coordinates, four.
56 
57  If residue_ is not nullptr, the endpoint is the named atom in the
58  pointed-to residue. (This is used for bonds that can span multiple
59  residues, perhaps in different chains, and is typically only created
60  by applying a two-residue patch such as DISU or LINK).
61 
62  If residue_ is nullptr, the endpoint is the named atom in the "current"
63  residue, unless the atom name has a +, -, 1: or 2: prefix.
64  Atoms prefixed with + are found in the next residue in the chain.
65  Atoms prefixed with - are found in the previous residue in the chain.
66  1: or 2: prefixes are used by two-residue patches to refer to atoms in
67  each residue patched by the patch.
68  */
69 class IMPATOMEXPORT CHARMMBondEndpoint {
70  std::string atom_name_;
71  Pointer<Object> residue_;
72 public:
73  CHARMMBondEndpoint(std::string atom_name,
74  CHARMMResidueTopology *residue=nullptr);
75 
76  std::string get_atom_name() const { return atom_name_; }
77 
78  //! Map the endpoint to an Atom particle.
79  Atom get_atom(const CHARMMResidueTopology *current_residue,
80  const CHARMMResidueTopology *previous_residue,
81  const CHARMMResidueTopology *next_residue,
82  const std::map<const CHARMMResidueTopology *,
83  Hierarchy> &resmap) const;
84 
85  IMP_SHOWABLE_INLINE(CHARMMBondEndpoint, {out << atom_name_;});
86 };
87 
89 
90 //! A connection (bond, angle, dihedral) between some number of endpoints.
91 template <unsigned int D>
93 {
94 protected:
96 public:
97  CHARMMConnection(const IMP::Strings &atoms) {
98  IMP_INTERNAL_CHECK(atoms.size() == D, "wrong number of bond endpoints");
99  for (Strings::const_iterator it = atoms.begin();
100  it != atoms.end(); ++it) {
101  endpoints_.push_back(CHARMMBondEndpoint(*it));
102  }
103  }
104 
105 #ifndef SWIG
107  : endpoints_(endpoints) {
108  IMP_INTERNAL_CHECK(endpoints.size() == D, "wrong number of bond endpoints");
109  }
110 #endif
111 
112  const IMP::atom::CHARMMBondEndpoint & get_endpoint(unsigned int i) const {
113  return endpoints_[i];
114  }
115 
116  //! \return true if the bond contains the named atom.
117  bool get_contains_atom(std::string name) const {
119  it = endpoints_.begin(); it != endpoints_.end(); ++it) {
120  if (it->get_atom_name() == name) {
121  return true;
122  }
123  }
124  return false;
125  }
126 
127 #ifndef SWIG
128  //! Map the bond to a list of Atom particles.
129  Atoms get_atoms(const CHARMMResidueTopology *current_residue,
130  const CHARMMResidueTopology *previous_residue,
131  const CHARMMResidueTopology *next_residue,
132  const std::map<const CHARMMResidueTopology *,
133  Hierarchy> &resmap) const {
134  Atoms as;
136  it = endpoints_.begin(); it != endpoints_.end(); ++it) {
137  Atom a = it->get_atom(current_residue, previous_residue,
138  next_residue, resmap);
139  if (a) {
140  as.push_back(a);
141  } else {
142  return Atoms();
143  }
144  }
145  return as;
146  }
147 #endif
150  it = endpoints_.begin(); it != endpoints_.end(); ++it) {
151  if (it != endpoints_.begin()) {
152  out << "-";
153  }
154  out << it->get_atom_name();
155  }
156  });
157 };
158 
159 typedef CHARMMConnection<2> CHARMMBond;
160 typedef CHARMMConnection<3> CHARMMAngle;
161 typedef CHARMMConnection<4> CHARMMDihedral;
162 
163 
164 IMP_VALUES(CHARMMBond, CHARMMBonds);
165 
166 IMP_VALUES(CHARMMAngle, CHARMMAngles);
167 
168 IMP_VALUES(CHARMMDihedral, CHARMMDihedrals);
169 
170 #ifdef SWIG
171 // Ugly, but SWIG needs a template instantiation before it is used as
172 // a base class
173 %template(CHARMMDihedral) CHARMMConnection<4>;
174 #endif
175 
176 //! A geometric relationship between four atoms.
177 /** The atoms (denoted i,j,k,l here) are uniquely positioned in 3D space
178  relative to each other by means of two distances, two angles, and
179  a dihedral.
180 
181  A regular internal coordinate stores the distances between ij and
182  kl respectively, and the angles between ijk and jkl.
183 
184  An improper internal coordinate stores the distances between ik and
185  kl respectively, and the angles between ikj and jkl.
186 
187  In both cases the dihedral is the angle between the two planes formed
188  by ijk and jkl.
189  */
191 {
192  float first_distance_, second_distance_, first_angle_, second_angle_,
193  dihedral_;
194  bool improper_;
195 public:
196  CHARMMInternalCoordinate(const IMP::Strings &atoms, float first_distance,
197  float first_angle, float dihedral,
198  float second_angle, float second_distance,
199  bool improper)
200  : CHARMMConnection<4>(atoms), first_distance_(first_distance),
201  second_distance_(second_distance), first_angle_(first_angle),
202  second_angle_(second_angle), dihedral_(dihedral),
203  improper_(improper) {}
205  float first_distance, float first_angle,
206  float dihedral, float second_angle,
207  float second_distance, bool improper)
208  : CHARMMConnection<4>(endpoints), first_distance_(first_distance),
209  second_distance_(second_distance), first_angle_(first_angle),
210  second_angle_(second_angle), dihedral_(dihedral),
211  improper_(improper) {}
212 
213  float get_first_distance() const { return first_distance_; }
214  float get_second_distance() const { return second_distance_; }
215  float get_first_angle() const { return first_angle_; }
216  float get_second_angle() const { return second_angle_; }
217  float get_dihedral() const { return dihedral_; }
218  bool get_improper() const { return improper_; }
219 
222  out << "; distances: " << first_distance_ << ", "
223  << second_distance_ << "; angles: " << first_angle_
224  << ", " << second_angle_ << "; dihedral: "
225  << dihedral_;
226  if (improper_) {
227  out << "; improper";
228  }
229  });
230 };
232 
233 //! Base class for all CHARMM residue-based topology
234 class IMPATOMEXPORT CHARMMResidueTopologyBase : public IMP::base::Object
235 {
236  std::string type_;
237 protected:
239  CHARMMBonds bonds_;
240  CHARMMAngles angles_;
241  CHARMMDihedrals dihedrals_;
242  CHARMMDihedrals impropers_;
243  CHARMMInternalCoordinates internal_coordinates_;
244 
245  CHARMMResidueTopologyBase(std::string type) :
246  Object("CharmmResidueTopoBase%1%"),
247  type_(type) {
248  set_name(std::string("CHARMM residue ") + type);
249  }
250 public:
251  std::string get_type() const { return type_; }
252 
253  unsigned int get_number_of_atoms() const { return atoms_.size(); }
254  const CHARMMAtomTopology &get_atom(unsigned int i) const { return atoms_[i]; }
255 
256  void add_atom(const CHARMMAtomTopology &atom);
257 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
258  CHARMMAtomTopology &get_atom(std::string name);
259 #endif
260  const CHARMMAtomTopology &get_atom(AtomType type) const {
261  return get_atom(type.get_string());
262  }
263 
264  const CHARMMAtomTopology &get_atom(std::string name) const;
265 
266  unsigned int get_number_of_bonds() const { return bonds_.size(); }
267  void add_bond(const CHARMMBond &bond) {
268  bonds_.push_back(bond);
269  }
270 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
271  // cannot return non const ref to swig
272  CHARMMBond &get_bond(unsigned int index) { return bonds_[index]; }
273 #endif
274 
275  unsigned int get_number_of_angles() const { return angles_.size(); }
276  void add_angle(const CHARMMAngle &bond) {
277  angles_.push_back(bond);
278  }
279 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
280  // cannot return non const ref to swig
281  CHARMMAngle &get_angle(unsigned int index) { return angles_[index]; }
282 #endif
283  unsigned int get_number_of_dihedrals() const { return dihedrals_.size(); }
284  void add_dihedral(const CHARMMDihedral &bond) {
285  dihedrals_.push_back(bond);
286  }
287 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
288  // cannot return non const ref to swig
289  CHARMMDihedral &get_dihedral(unsigned int index) { return dihedrals_[index]; }
290 #endif
291  unsigned int get_number_of_impropers() const { return impropers_.size(); }
292  void add_improper(const CHARMMDihedral &bond) {
293  impropers_.push_back(bond);
294  }
295 #if !defined(SWIG) && !defined(IMP_DOXYGEN)
296  CHARMMDihedral &get_improper(unsigned int index) { return impropers_[index]; }
297 #endif
298 
299  unsigned int get_number_of_internal_coordinates() const {
300  return internal_coordinates_.size();
301  }
302  void add_internal_coordinate(const CHARMMInternalCoordinate &ic) {
303  internal_coordinates_.push_back(ic);
304  }
306  &get_internal_coordinate(unsigned int index) const {
307  return internal_coordinates_[index];
308  }
309 
310  const CHARMMBond &get_bond(unsigned int index) const {
311  return bonds_[index];
312  }
313  const CHARMMAngle &get_angle(unsigned int index) const {
314  return angles_[index];
315  }
316  const CHARMMDihedral &get_dihedral(unsigned int index) const {
317  return dihedrals_[index];
318  }
319  const CHARMMDihedral &get_improper(unsigned int index) const {
320  return impropers_[index];
321  }
322 
324 };
325 
327 
328 //! The ideal topology of a single residue
329 /** These residue topologies can be constructed manually (by adding
330  CHARMMAtomTopology objects with add_atom()). However, they are more
331  commonly populated automatically by reading a CHARMM topology file
332  (CHARMMParameters::get_residue_topology()).
333  */
334 class IMPATOMEXPORT CHARMMIdealResidueTopology
335  : public CHARMMResidueTopologyBase {
336  std::string default_first_patch_, default_last_patch_;
337 public:
338  CHARMMIdealResidueTopology(std::string type)
339  : CHARMMResidueTopologyBase(type) {}
341  : CHARMMResidueTopologyBase(type.get_string()) {}
342 
343  //! Remove the named atom
344  /** Any bonds/angles that involve this atom are also removed.
345  */
346  void remove_atom(std::string name);
347 
348  void set_default_first_patch(std::string patch) {
349  default_first_patch_ = patch;
350  }
351  void set_default_last_patch(std::string patch) {
352  default_last_patch_ = patch;
353  }
354  std::string get_default_first_patch() const { return default_first_patch_; }
355  std::string get_default_last_patch() const { return default_last_patch_; }
357 };
358 
360 
362 
363 //! A CHARMM patch residue.
364 /** Patch residues are similar to regular residues, except that they are
365  used to modify an existing residue. Any atoms they contain replace or
366  add to those in the residue; they can also remove atoms.
367  Atom names are as for regular residues, except for patches that affect
368  two residues (e.g. DISU, LINK) in which each atom contains a 1 or 2
369  prefix (plus an optional colon) to identify the residue which will be
370  patched; for example, both 1:CA and 1CA identify the CA atom in the
371  first residue.
372  */
373 class IMPATOMEXPORT CHARMMPatch : public CHARMMResidueTopologyBase {
374  base::Vector<std::string> deleted_atoms_;
375 public:
376  //! Construct a new, empty patch residue.
377  /** To get an existing patch, use CHARMMParameters::get_patch() instead.
378  */
379  CHARMMPatch(std::string type) : CHARMMResidueTopologyBase(type) {
380  set_name(std::string("CHARMM patching residue ") + type);
381  }
382 
383  void add_removed_atom(std::string name) { deleted_atoms_.push_back(name); }
384 
385  unsigned int get_number_of_removed_atoms() const {
386  return deleted_atoms_.size();
387  }
388 
389  std::string get_removed_atom(unsigned int i) const {
390  return deleted_atoms_[i];
391  }
392 
393  //! Apply the patch to the residue, modifying its topology accordingly.
394  /** \note Most CHARMM patches are designed to be applied in isolation;
395  it is usually an error to try to apply two different patches
396  to the same residue. Thus, by default \imp prohibits this.
397  To allow an already-patched residue to be re-patched, first
398  call CHARMMResidueTopology::set_patched(false).
399  */
400  void apply(CHARMMResidueTopology *res) const;
401 
402  //! Apply the patch to the given pair of residues.
403  /** This can only be used for special two-residue patches, such as
404  DISU or LINK. In a two-residue patch, each atom has a 1: or 2: prefix
405  to identify the residue it refers to.
406 
407  \throws ValueException if the patch is not a two-residue patch.
408  */
409  void apply(CHARMMResidueTopology *res1, CHARMMResidueTopology *res2) const;
410  IMP_OBJECT(CHARMMPatch);
411 };
412 
414 
415 
416 //! The topology of a single residue in a model.
417 /** Each CHARMMResidueTopology object can represent an 'unknown' residue
418  (containing no atoms or bonds) or is a copy of an existing
419  CHARMMIdealResidueTopology. Optionally, patches can be applied to
420  residues (see CHARMMPatch::apply()) to add N- or C-termini, disulfide
421  bridges, sidechain modifications, etc.
422  */
424 {
425  bool patched_;
426 public:
427 
428  //! Create an empty topology, containing no atoms or bonds.
430  : CHARMMIdealResidueTopology(type), patched_(false) {}
431 
432  //! Construct residue topology as a copy of an existing topology.
434  : CHARMMIdealResidueTopology(ideal->get_type()), patched_(false) {
435  set_default_first_patch(ideal->get_default_first_patch());
436  set_default_last_patch(ideal->get_default_last_patch());
437 
438  // Add atoms from existing topology
439  for (unsigned int i = 0; i < ideal->get_number_of_atoms(); ++i) {
440  add_atom(ideal->get_atom(i));
441  }
442  // Add angles/bonds/dihedrals/impropers
443  for (unsigned int i = 0; i < ideal->get_number_of_bonds(); ++i) {
444  add_bond(ideal->get_bond(i));
445  }
446  for (unsigned int i = 0; i < ideal->get_number_of_angles(); ++i) {
447  add_angle(ideal->get_angle(i));
448  }
449  for (unsigned int i = 0; i < ideal->get_number_of_dihedrals(); ++i) {
450  add_dihedral(ideal->get_dihedral(i));
451  }
452  for (unsigned int i = 0; i < ideal->get_number_of_impropers(); ++i) {
453  add_improper(ideal->get_improper(i));
454  }
455  for (unsigned int i = 0; i < ideal->get_number_of_internal_coordinates();
456  ++i) {
457  add_internal_coordinate(ideal->get_internal_coordinate(i));
458  }
459  }
460 
461  bool get_patched() const { return patched_; }
462  void set_patched(bool patched) { patched_ = patched; }
463 
464  IMP_OBJECT(CHARMMResidueTopology);
465 };
466 
468 
469 IMPATOM_END_NAMESPACE
470 
471 #endif /* IMPATOM_CHARMM_TOPOLOGY_H */