IMP  2.2.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  CHARMMTopology *create_topology(Hierarchy hierarchy) const;
172 
173  //! Get bond parameters for the bond between the two given CHARMM atom types.
174  /** The atom types may match in any order.
175  \throws IndexException if no parameters are present.
176  */
177  const CHARMMBondParameters &get_bond_parameters(std::string type1,
178  std::string type2) const {
179  internal::CHARMMBondNames types = internal::CHARMMBondNames(type1, type2);
180  if (bond_parameters_.find(types) != bond_parameters_.end()) {
181  return bond_parameters_.find(types)->second;
182  } else {
183  IMP_THROW("No CHARMM parameters found for bond " << type1 << "-" << type2,
184  IndexException);
185  }
186  }
187 
188  //! Get parameters for the angle between the three given CHARMM atom types.
189  /** The atom types may match in either forward or reverse order.
190  \throws IndexException if no parameters are present.
191  */
192  const CHARMMBondParameters &get_angle_parameters(std::string type1,
193  std::string type2,
194  std::string type3) const {
195  internal::CHARMMAngleNames types =
196  internal::CHARMMAngleNames(type1, type2, type3);
197  if (angle_parameters_.find(types) != angle_parameters_.end()) {
198  return angle_parameters_.find(types)->second;
199  } else {
200  IMP_THROW("No CHARMM parameters found for angle " << type1 << "-" << type2
201  << "-" << type3,
202  IndexException);
203  }
204  }
205 
206  //! Get parameters for the dihedral between the four given CHARMM atom types.
207  /** The atom types may match in either forward or reverse order. When
208  looking for a match in the library, wildcards are considered; an atom
209  type of X in the library will match any atom type. The most specific
210  match from the library is returned.
211 
212  Multiple sets of parameters can be specified for the same combination
213  of atom types in the library, in which case all of them are returned.
214 
215  \throws IndexException if no parameters are present.
216  */
218  std::string type1, std::string type2, std::string type3,
219  std::string type4) const {
221  internal::CHARMMDihedralNames types =
222  internal::CHARMMDihedralNames(type1, type2, type3, type4);
223  // Get the first match, using wildcards
224  DihedralParameters::const_iterator match = find_dihedral(
225  dihedral_parameters_.begin(), dihedral_parameters_.end(), types, true);
226  if (match != dihedral_parameters_.end()) {
227  // If it matched, look for duplicate dihedral terms (this time the
228  // match must be exactly the same as the first match)
229  param.push_back(match->second);
230  while ((match = find_dihedral(match + 1, dihedral_parameters_.end(),
231  match->first, false)) !=
232  dihedral_parameters_.end()) {
233  param.push_back(match->second);
234  }
235  }
236  if (param.size() == 0) {
237  IMP_THROW("No CHARMM parameters found for dihedral "
238  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
239  IndexException);
240  } else {
241  return param;
242  }
243  }
244 
245  //! Get parameters for the improper between the four given CHARMM atom types.
246  /** The atom types may match in either forward or reverse order. When
247  looking for a match in the library, wildcards are considered; an atom
248  type of X in the library will match any atom type. The most specific
249  match from the library is returned.
250 
251  \throws IndexException if no parameters are present.
252  */
254  std::string type1, std::string type2, std::string type3,
255  std::string type4) const {
256  internal::CHARMMDihedralNames types =
257  internal::CHARMMDihedralNames(type1, type2, type3, type4);
258  // Return just the first match; wildcards are OK
259  DihedralParameters::const_iterator it = find_dihedral(
260  improper_parameters_.begin(), improper_parameters_.end(), types, true);
261  if (it != improper_parameters_.end()) {
262  return it->second;
263  } else {
264  IMP_THROW("No CHARMM parameters found for improper "
265  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
266  IndexException);
267  }
268  }
269 
270  //! Auto-generate Angle particles from the passed list of Bond particles.
271  /** The angles consist of all unique pairs of bonds which share an
272  endpoint. If no parameters are found for an angle, it is simply
273  created without those parameters.
274 
275  The list of newly-created Angle particles can be passed to a
276  StereochemistryPairFilter to exclude 1-3 interactions from the
277  nonbonded list, or to an AngleSingletonScore to score each angle.
278 
279  \return a list of the newly-created Angle particles.
280 
281  \see CHARMMTopology::add_bonds().
282  */
283  kernel::Particles create_angles(kernel::Particles bonds) const;
284 
285  //! Auto-generate Dihedral particles from the passed list of Bond particles.
286  /** The dihedrals consist of all unique triples of bonds which form
287  dihedrals. If no parameters are found for a dihedral, it is simply
288  created without those parameters; if multiple sets of parameters are
289  found, multiple copies of the dihedral are created, each with one set
290  of parameters.
291 
292  The list of newly-created Dihedral particles can be passed to a
293  StereochemistryPairFilter to exclude 1-4 interactions from the
294  nonbonded list, or to a DihedralSingletonScore to score each dihedral.
295 
296  If dihedrals are explicitly listed in the CHARMM topology file, they
297  can be created if desired by calling CHARMMTopology::add_dihedrals()
298  rather than this function.
299 
300  \return a list of the newly-created Dihedral particles.
301 
302  \see CHARMMTopology::add_bonds().
303  */
304  kernel::Particles create_dihedrals(kernel::Particles bonds) const;
305 
307 
308  private:
309  virtual String get_force_field_atom_type(Atom atom) const;
310 
311  void read_parameter_file(base::TextInput input_file);
312  // read topology file
313  void read_topology_file(base::TextInput input_file,
314  bool translate_names_to_pdb);
315 
316  void add_angle(kernel::Particle *p1, kernel::Particle *p2,
317  kernel::Particle *p3, kernel::Particles &ps) const;
318  void add_dihedral(kernel::Particle *p1, kernel::Particle *p2,
320  kernel::Particles &ps) const;
321 
322  ResidueType parse_residue_line(const String &line,
323  bool translate_names_to_pdb);
324  void parse_atom_line(const String &line, ResidueType curr_res_type,
325  CHARMMResidueTopologyBase *residue,
326  bool translate_names_to_pdb);
327  void parse_bond_line(const String &line, ResidueType curr_res_type,
328  CHARMMResidueTopologyBase *residue,
329  bool translate_names_to_pdb);
330 
331  void parse_nonbonded_parameters_line(String line);
332  void parse_bonds_parameters_line(String line);
333  void parse_angles_parameters_line(String line);
334  void parse_dihedrals_parameters_line(String line, DihedralParameters &param);
335  base::WarningContext warn_context_;
336 };
337 
339 
340 /** The default CHARMM parameters support normal amino acid
341  and nucleic acid residues and the atoms found in them.
342  To use CHARMM with heterogens or non-standard residues,
343  a different CHARMM parameters file must be used.
344 
345  No hydrogen parameters are read.
346 
347  \see get_all_atom_CHARMM_parameters();
348 */
350 
351 /** The default CHARMM parameters support normal amino acid
352  and nucleic acid residues and the atoms found in them.
353  To use CHARMM with heterogens or non-standard residues,
354  a different CHARMM parameters file must be used.
355 
356  \see get_heavy_atom_CHARMM_parameters()
357 */
359 
360 IMPATOM_END_NAMESPACE
361 
362 #endif /* IMPATOM_CHARMM_PARAMETERS_H */
The ideal topology of a single residue.
CHARMMParameters * get_heavy_atom_CHARMM_parameters()
Storage and access to force field.
void set_was_used(bool tf) const
The parameters for a CHARMM bond or angle.
The parameters for a CHARMM dihedral or improper.
#define IMP_VALUES(Name, PluralName)
Define the type for storing sets of values.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
force field base class
A smart pointer to a reference counted object.
Definition: base/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.
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.
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
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Handling of file input/output.
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: base/types.h:44