IMP  2.3.1
The Integrative Modeling Platform
CHARMMParameters.h
Go to the documentation of this file.
1 /**
2  * \file IMP/atom/CHARMMParameters.h
3  * \brief access to CHARMM force field parameters
4  *
5  * Copyright 2007-2014 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPATOM_CHARMM_PARAMETERS_H
10 #define IMPATOM_CHARMM_PARAMETERS_H
11 
12 #include "internal/charmm_helpers.h"
13 #include "ForceFieldParameters.h"
14 #include "charmm_topology.h"
15 #include "atom_macros.h"
16 #include <IMP/base/file.h>
17 
18 #include <string>
19 // swig is being dumb
20 #ifdef SWIG
21 IMPKERNEL_BEGIN_NAMESPACE
22 class VersionInfo;
23 IMPKERNEL_END_NAMESPACE
24 #endif
25 
26 IMPATOM_BEGIN_NAMESPACE
27 
28 //! The parameters for a CHARMM bond or angle.
30  double force_constant;
31  double ideal;
34  out << "force constant: " << force_constant << "; ideal value: " << ideal;
35  });
36 };
37 
39 
40 //! The parameters for a CHARMM dihedral or improper.
42  double force_constant;
43  int multiplicity;
44  double ideal;
47  out << "force constant: " << force_constant
48  << "; multiplicity: " << multiplicity << "; ideal value: " << ideal;
49  });
50 };
51 
53 
54 class CHARMMTopology;
55 
56 //! CHARMM force field parameters.
57 /** This class reads in topology and parameter files in CHARMM format and
58  stores the information.
59 
60  It does not actually evaluate the force field itself - there are other
61  classes that use the parameters provided by this class to do that. For
62  example, the LennardJonesPairScore or CoulombPairScore evaluate the
63  nonbond terms of the CHARMM force field, while
64  CHARMMStereochemistryRestraint (or BondSingletonScore, AngleSingletonScore,
65  DihedralSingletonScore and ImproperSingletonScore) cover the bond terms.
66 
67  Typically, the create_topology() method is used to create a new
68  CHARMMTopology object for a given Hierarchy; that object can then be
69  used to assign atomic radii, bonds, etc.
70  */
71 class IMPATOMEXPORT CHARMMParameters : public ForceFieldParameters {
72  std::map<std::string, Element> atom_type_to_element_;
73  std::map<ResidueType, base::Pointer<CHARMMIdealResidueTopology> >
74  residue_topologies_;
75  std::map<std::string, base::Pointer<CHARMMPatch> > patches_;
76  std::map<internal::CHARMMBondNames, CHARMMBondParameters> bond_parameters_;
77  std::map<internal::CHARMMAngleNames, CHARMMBondParameters> angle_parameters_;
78 
79  typedef base::Vector<std::pair<internal::CHARMMDihedralNames,
81  DihedralParameters dihedral_parameters_;
82  DihedralParameters improper_parameters_;
83 
84  DihedralParameters::const_iterator find_dihedral(
85  DihedralParameters::const_iterator begin,
86  DihedralParameters::const_iterator end,
87  const internal::CHARMMDihedralNames &dihedral,
88  bool allow_wildcards) const;
89 
90  public:
91  //! Construction with CHARMM topology (and optionally parameters) file.
92  /** For addition of atom types, the topology file alone is enough;
93  for adding bonds and radii, both files are needed.
94 
95  Atom and residue naming in the topology file differs slightly between
96  CHARMM and PDB. If translate_names_to_pdb is set to true, some simple
97  translations are done to map CHARMM-style names to PDB-style for
98  standard amino acids and some other commonly-used residues and patches.
99  (This translation has already been done for the topology files included
100  with IMP and MODELLER, so it is only needed for topology files taken
101  from CHARMM itself or other sources.) The modifications are as follows:
102  - CHARMM HSD (unprotonated histidine) is mapped to PDB HIS.
103  - CD1 and CD2 atoms in LEU are swapped.
104  - OT1 and OT2 in CTER are mapped to O and OXT.
105  - CHARMM hydrogen names are mapped to PDB equivalents.
106  - CHARMM NTER, GLYP and CTER residues are modified slightly to avoid
107  removing the HN, HN and O atoms respectively, and adding excess bonds
108  to these atoms.
109  */
110  CHARMMParameters(base::TextInput topology_file_name,
111  base::TextInput par_file_name = base::TextInput(),
112  bool translate_names_to_pdb = false);
113 
114  /** \name Residue topology
115 
116  The class stores the topology of each residue type as defined in the
117  topology file, as a set of CHARMMIdealResidueTopology objects.
118  */
119  /**@{*/
120  void add_residue_topology(CHARMMIdealResidueTopology *res) {
121  res->set_was_used(true);
122  residue_topologies_.insert(
123  std::make_pair(ResidueType(res->get_type()), res));
124  }
125 
126  CHARMMIdealResidueTopology *get_residue_topology(ResidueType type) const {
127  std::map<ResidueType,
128  base::Pointer<CHARMMIdealResidueTopology> >::const_iterator it =
129  residue_topologies_.find(type);
130  if (it != residue_topologies_.end()) {
131  return it->second;
132  } else {
133  IMP_THROW("Residue " << type << " does not exist", ValueException);
134  }
135  }
136 
137  /** \name Patches
138 
139  The class stores patches as defined in the topology file, as a set of
140  CHARMMPatch objects.
141  */
142  /**@{*/
143  void add_patch(CHARMMPatch *patch) {
144  patch->set_was_used(true);
145  patches_.insert(std::make_pair(patch->get_type(), patch));
146  }
147 #if 0
148  // return of non const ref values is not allowed
149  CHARMMPatch &get_patch(std::string name) {
150  std::map<std::string, CHARMMPatch>::iterator it = patches_.find(name);
151  if (it != patches_.end()) {
152  return it->second;
153  } else {
154  IMP_THROW("Patch " << name << " does not exist", ValueException);
155  }
156  }
157 #endif
158 
159  CHARMMPatch *get_patch(std::string name) const {
160  std::map<std::string, base::Pointer<CHARMMPatch> >::const_iterator it =
161  patches_.find(name);
162  if (it != patches_.end()) {
163  return it->second;
164  } else {
165  IMP_THROW("Patch " << name << " does not exist", ValueException);
166  }
167  }
168  /**@}*/
169 
170  //! Create topology that corresponds to the primary sequence of the Hierarchy.
171  /** Residues are placed in different segments if they are deemed to be
172  "disconnected", i.e. they have different Chains as parents, they are
173  in different Fragments and the fragments are not consecutive (the last
174  residue index of the first Fragment and the first index of the second
175  Fragment are not consecutive), or their ancestors are different
176  (e.g. one Residue lives in a Chain and another does not have a Chain
177  parent) */
178  CHARMMTopology *create_topology(Hierarchy hierarchy) const;
179 
180  //! Get bond parameters for the bond between the two given CHARMM atom types.
181  /** The atom types may match in any order.
182  \throws IndexException if no parameters are present.
183  */
184  const CHARMMBondParameters &get_bond_parameters(std::string type1,
185  std::string type2) const {
186  internal::CHARMMBondNames types = internal::CHARMMBondNames(type1, type2);
187  if (bond_parameters_.find(types) != bond_parameters_.end()) {
188  return bond_parameters_.find(types)->second;
189  } else {
190  IMP_THROW("No CHARMM parameters found for bond " << type1 << "-" << type2,
191  IndexException);
192  }
193  }
194 
195  //! Get parameters for the angle between the three given CHARMM atom types.
196  /** The atom types may match in either forward or reverse order.
197  \throws IndexException if no parameters are present.
198  */
199  const CHARMMBondParameters &get_angle_parameters(std::string type1,
200  std::string type2,
201  std::string type3) const {
202  internal::CHARMMAngleNames types =
203  internal::CHARMMAngleNames(type1, type2, type3);
204  if (angle_parameters_.find(types) != angle_parameters_.end()) {
205  return angle_parameters_.find(types)->second;
206  } else {
207  IMP_THROW("No CHARMM parameters found for angle " << type1 << "-" << type2
208  << "-" << type3,
209  IndexException);
210  }
211  }
212 
213  //! Get parameters for the dihedral between the four given CHARMM atom types.
214  /** The atom types may match in either forward or reverse order. When
215  looking for a match in the library, wildcards are considered; an atom
216  type of X in the library will match any atom type. The most specific
217  match from the library is returned.
218 
219  Multiple sets of parameters can be specified for the same combination
220  of atom types in the library, in which case all of them are returned.
221 
222  \throws IndexException if no parameters are present.
223  */
225  std::string type1, std::string type2, std::string type3,
226  std::string type4) const {
228  internal::CHARMMDihedralNames types =
229  internal::CHARMMDihedralNames(type1, type2, type3, type4);
230  // Get the first match, using wildcards
231  DihedralParameters::const_iterator match = find_dihedral(
232  dihedral_parameters_.begin(), dihedral_parameters_.end(), types, true);
233  if (match != dihedral_parameters_.end()) {
234  // If it matched, look for duplicate dihedral terms (this time the
235  // match must be exactly the same as the first match)
236  param.push_back(match->second);
237  while ((match = find_dihedral(match + 1, dihedral_parameters_.end(),
238  match->first, false)) !=
239  dihedral_parameters_.end()) {
240  param.push_back(match->second);
241  }
242  }
243  if (param.size() == 0) {
244  IMP_THROW("No CHARMM parameters found for dihedral "
245  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
246  IndexException);
247  } else {
248  return param;
249  }
250  }
251 
252  //! Get parameters for the improper between the four given CHARMM atom types.
253  /** The atom types may match in either forward or reverse order. When
254  looking for a match in the library, wildcards are considered; an atom
255  type of X in the library will match any atom type. The most specific
256  match from the library is returned.
257 
258  \throws IndexException if no parameters are present.
259  */
261  std::string type1, std::string type2, std::string type3,
262  std::string type4) const {
263  internal::CHARMMDihedralNames types =
264  internal::CHARMMDihedralNames(type1, type2, type3, type4);
265  // Return just the first match; wildcards are OK
266  DihedralParameters::const_iterator it = find_dihedral(
267  improper_parameters_.begin(), improper_parameters_.end(), types, true);
268  if (it != improper_parameters_.end()) {
269  return it->second;
270  } else {
271  IMP_THROW("No CHARMM parameters found for improper "
272  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
273  IndexException);
274  }
275  }
276 
277  //! Auto-generate Angle particles from the passed list of Bond particles.
278  /** The angles consist of all unique pairs of bonds which share an
279  endpoint. If no parameters are found for an angle, it is simply
280  created without those parameters.
281 
282  The list of newly-created Angle particles can be passed to a
283  StereochemistryPairFilter to exclude 1-3 interactions from the
284  nonbonded list, or to an AngleSingletonScore to score each angle.
285 
286  \return a list of the newly-created Angle particles.
287 
288  \see CHARMMTopology::add_bonds().
289  */
290  kernel::Particles create_angles(kernel::Particles bonds) const;
291 
292  //! Auto-generate Dihedral particles from the passed list of Bond particles.
293  /** The dihedrals consist of all unique triples of bonds which form
294  dihedrals. If no parameters are found for a dihedral, it is simply
295  created without those parameters; if multiple sets of parameters are
296  found, multiple copies of the dihedral are created, each with one set
297  of parameters.
298 
299  The list of newly-created Dihedral particles can be passed to a
300  StereochemistryPairFilter to exclude 1-4 interactions from the
301  nonbonded list, or to a DihedralSingletonScore to score each dihedral.
302 
303  If dihedrals are explicitly listed in the CHARMM topology file, they
304  can be created if desired by calling CHARMMTopology::add_dihedrals()
305  rather than this function.
306 
307  \return a list of the newly-created Dihedral particles.
308 
309  \see CHARMMTopology::add_bonds().
310  */
311  kernel::Particles create_dihedrals(kernel::Particles bonds) const;
312 
314 
315  private:
316  virtual String get_force_field_atom_type(Atom atom) const;
317 
318  void read_parameter_file(base::TextInput input_file);
319  // read topology file
320  void read_topology_file(base::TextInput input_file,
321  bool translate_names_to_pdb);
322 
323  void add_angle(kernel::Particle *p1, kernel::Particle *p2,
324  kernel::Particle *p3, kernel::Particles &ps) const;
325  void add_dihedral(kernel::Particle *p1, kernel::Particle *p2,
327  kernel::Particles &ps) const;
328 
329  ResidueType parse_residue_line(const String &line,
330  bool translate_names_to_pdb);
331  void parse_atom_line(const String &line, ResidueType curr_res_type,
332  CHARMMResidueTopologyBase *residue,
333  bool translate_names_to_pdb);
334  void parse_bond_line(const String &line, ResidueType curr_res_type,
335  CHARMMResidueTopologyBase *residue,
336  bool translate_names_to_pdb);
337 
338  void parse_nonbonded_parameters_line(String line);
339  void parse_bonds_parameters_line(String line);
340  void parse_angles_parameters_line(String line);
341  void parse_dihedrals_parameters_line(String line, DihedralParameters &param);
342  base::WarningContext warn_context_;
343 };
344 
346 
347 /** The default CHARMM parameters support normal amino acid
348  and nucleic acid residues and the atoms found in them.
349  To use CHARMM with heterogens or non-standard residues,
350  a different CHARMM parameters file must be used.
351 
352  No hydrogen parameters are read.
353 
354  \see get_all_atom_CHARMM_parameters();
355 */
357 
358 /** The default CHARMM parameters support normal amino acid
359  and nucleic acid residues and the atoms found in them.
360  To use CHARMM with heterogens or non-standard residues,
361  a different CHARMM parameters file must be used.
362 
363  \see get_heavy_atom_CHARMM_parameters()
364 */
366 
367 IMPATOM_END_NAMESPACE
368 
369 #endif /* IMPATOM_CHARMM_PARAMETERS_H */
The ideal topology of a single residue.
CHARMMParameters * get_heavy_atom_CHARMM_parameters()
Storage and access to force field.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
void set_was_used(bool tf) const
The parameters for a CHARMM bond or angle.
The parameters for a CHARMM dihedral or improper.
force field base class
Handling of file input/output.
A smart pointer to a reference counted object.
Definition: Pointer.h:87
const CHARMMBondParameters & get_angle_parameters(std::string type1, std::string type2, std::string type3) const
Get parameters for the angle between the three given CHARMM atom types.
const CHARMMBondParameters & get_bond_parameters(std::string type1, std::string type2) const
Get bond parameters for the bond between the two given CHARMM atom types.
CHARMM force field parameters.
Classes for handling CHARMM-style topology.
const CHARMMDihedralParameters & get_improper_parameters(std::string type1, std::string type2, std::string type3, std::string type4) const
Get parameters for the improper between the four given CHARMM atom types.
#define IMP_VALUES(Name, PluralName)
Define the type for storing sets of values.
Definition: value_macros.h:23
The standard decorator for manipulating molecular structures.
Base class for all CHARMM residue-based topology.
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
The type for a residue.
Class to handle individual model particles.
A CHARMM patch residue.
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing sets of objects.
Definition: object_macros.h:52
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
Various important macros for implementing decorators.
#define IMP_FORCE_FIELD_PARAMETERS(Name)
Define the basic things you need for a ForceFieldParameters.
Definition: atom_macros.h:22
CHARMMDihedralParametersList get_dihedral_parameters(std::string type1, std::string type2, std::string type3, std::string type4) const
Get parameters for the dihedral between the four given CHARMM atom types.
The topology of a complete CHARMM model.
CHARMMParameters * get_all_atom_CHARMM_parameters()
std::string String
Basic string value.
Definition: types.h:44