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