IMP  2.1.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-2013 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 
92  //! Construction with CHARMM topology (and optionally parameters) file.
93  /** For addition of atom types, the topology file alone is enough;
94  for adding bonds and radii, both files are needed.
95 
96  Atom and residue naming in the topology file differs slightly between
97  CHARMM and PDB. If translate_names_to_pdb is set to true, some simple
98  translations are done to map CHARMM-style names to PDB-style for
99  standard amino acids and some other commonly-used residues and patches.
100  (This translation has already been done for the topology files included
101  with IMP and MODELLER, so it is only needed for topology files taken
102  from CHARMM itself or other sources.) The modifications are as follows:
103  - CHARMM HSD (unprotonated histidine) is mapped to PDB HIS.
104  - CD1 and CD2 atoms in LEU are swapped.
105  - OT1 and OT2 in CTER are mapped to O and OXT.
106  - CHARMM hydrogen names are mapped to PDB equivalents.
107  - CHARMM NTER, GLYP and CTER residues are modified slightly to avoid
108  removing the HN, HN and O atoms respectively, and adding excess bonds
109  to these atoms.
110  */
111  CHARMMParameters(base::TextInput topology_file_name,
112  base::TextInput par_file_name = base::TextInput(),
113  bool translate_names_to_pdb = false);
114 
115  /** \name Residue topology
116 
117  The class stores the topology of each residue type as defined in the
118  topology file, as a set of CHARMMIdealResidueTopology objects.
119  */
120  /**@{*/
121  void add_residue_topology(CHARMMIdealResidueTopology *res) {
122  res->set_was_used(true);
123  residue_topologies_.insert(
124  std::make_pair(ResidueType(res->get_type()), res));
125  }
126 
127  CHARMMIdealResidueTopology *get_residue_topology(ResidueType type) const {
128  std::map<ResidueType,
129  base::Pointer<CHARMMIdealResidueTopology> >::const_iterator it =
130  residue_topologies_.find(type);
131  if (it != residue_topologies_.end()) {
132  return it->second;
133  } else {
134  IMP_THROW("Residue " << type << " does not exist", ValueException);
135  }
136  }
137 
138  /** \name Patches
139 
140  The class stores patches as defined in the topology file, as a set of
141  CHARMMPatch objects.
142  */
143  /**@{*/
144  void add_patch(CHARMMPatch *patch) {
145  patch->set_was_used(true);
146  patches_.insert(std::make_pair(patch->get_type(), patch));
147  }
148 #if 0
149  // return of non const ref values is not allowed
150  CHARMMPatch &get_patch(std::string name) {
151  std::map<std::string, CHARMMPatch>::iterator it = patches_.find(name);
152  if (it != patches_.end()) {
153  return it->second;
154  } else {
155  IMP_THROW("Patch " << name << " does not exist", ValueException);
156  }
157  }
158 #endif
159 
160  CHARMMPatch *get_patch(std::string name) const {
161  std::map<std::string, base::Pointer<CHARMMPatch> >::const_iterator it =
162  patches_.find(name);
163  if (it != patches_.end()) {
164  return it->second;
165  } else {
166  IMP_THROW("Patch " << name << " does not exist", ValueException);
167  }
168  }
169  /**@}*/
170 
171  //! Create topology that corresponds to the primary sequence of the Hierarchy.
172  CHARMMTopology *create_topology(Hierarchy hierarchy) const;
173 
174  //! Get bond parameters for the bond between the two given CHARMM atom types.
175  /** The atom types may match in any order.
176  \throws IndexException if no parameters are present.
177  */
178  const CHARMMBondParameters &get_bond_parameters(std::string type1,
179  std::string type2) const {
180  internal::CHARMMBondNames types = internal::CHARMMBondNames(type1, type2);
181  if (bond_parameters_.find(types) != bond_parameters_.end()) {
182  return bond_parameters_.find(types)->second;
183  } else {
184  IMP_THROW("No CHARMM parameters found for bond " << type1 << "-" << type2,
185  IndexException);
186  }
187  }
188 
189  //! Get parameters for the angle between the three given CHARMM atom types.
190  /** The atom types may match in either forward or reverse order.
191  \throws IndexException if no parameters are present.
192  */
193  const CHARMMBondParameters &get_angle_parameters(std::string type1,
194  std::string type2,
195  std::string type3) const {
196  internal::CHARMMAngleNames types =
197  internal::CHARMMAngleNames(type1, type2, type3);
198  if (angle_parameters_.find(types) != angle_parameters_.end()) {
199  return angle_parameters_.find(types)->second;
200  } else {
201  IMP_THROW("No CHARMM parameters found for angle " << type1 << "-" << type2
202  << "-" << type3,
203  IndexException);
204  }
205  }
206 
207  //! Get parameters for the dihedral between the four given CHARMM atom types.
208  /** The atom types may match in either forward or reverse order. When
209  looking for a match in the library, wildcards are considered; an atom
210  type of X in the library will match any atom type. The most specific
211  match from the library is returned.
212 
213  Multiple sets of parameters can be specified for the same combination
214  of atom types in the library, in which case all of them are returned.
215 
216  \throws IndexException if no parameters are present.
217  */
219  std::string type1, std::string type2, std::string type3,
220  std::string type4) const {
222  internal::CHARMMDihedralNames types =
223  internal::CHARMMDihedralNames(type1, type2, type3, type4);
224  // Get the first match, using wildcards
225  DihedralParameters::const_iterator match = find_dihedral(
226  dihedral_parameters_.begin(), dihedral_parameters_.end(), types, true);
227  if (match != dihedral_parameters_.end()) {
228  // If it matched, look for duplicate dihedral terms (this time the
229  // match must be exactly the same as the first match)
230  param.push_back(match->second);
231  while ((match = find_dihedral(match + 1, dihedral_parameters_.end(),
232  match->first, false)) !=
233  dihedral_parameters_.end()) {
234  param.push_back(match->second);
235  }
236  }
237  if (param.size() == 0) {
238  IMP_THROW("No CHARMM parameters found for dihedral "
239  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
240  IndexException);
241  } else {
242  return param;
243  }
244  }
245 
246  //! Get parameters for the improper between the four given CHARMM atom types.
247  /** The atom types may match in either forward or reverse order. When
248  looking for a match in the library, wildcards are considered; an atom
249  type of X in the library will match any atom type. The most specific
250  match from the library is returned.
251 
252  \throws IndexException if no parameters are present.
253  */
255  std::string type1, std::string type2, std::string type3,
256  std::string type4) const {
257  internal::CHARMMDihedralNames types =
258  internal::CHARMMDihedralNames(type1, type2, type3, type4);
259  // Return just the first match; wildcards are OK
260  DihedralParameters::const_iterator it = find_dihedral(
261  improper_parameters_.begin(), improper_parameters_.end(), types, true);
262  if (it != improper_parameters_.end()) {
263  return it->second;
264  } else {
265  IMP_THROW("No CHARMM parameters found for improper "
266  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
267  IndexException);
268  }
269  }
270 
271  //! Auto-generate Angle particles from the passed list of Bond particles.
272  /** The angles consist of all unique pairs of bonds which share an
273  endpoint. If no parameters are found for an angle, it is simply
274  created without those parameters.
275 
276  The list of newly-created Angle particles can be passed to a
277  StereochemistryPairFilter to exclude 1-3 interactions from the
278  nonbonded list, or to an AngleSingletonScore to score each angle.
279 
280  \return a list of the newly-created Angle particles.
281 
282  \see CHARMMTopology::add_bonds().
283  */
284  kernel::Particles create_angles(kernel::Particles bonds) const;
285 
286  //! Auto-generate Dihedral particles from the passed list of Bond particles.
287  /** The dihedrals consist of all unique triples of bonds which form
288  dihedrals. If no parameters are found for a dihedral, it is simply
289  created without those parameters; if multiple sets of parameters are
290  found, multiple copies of the dihedral are created, each with one set
291  of parameters.
292 
293  The list of newly-created Dihedral particles can be passed to a
294  StereochemistryPairFilter to exclude 1-4 interactions from the
295  nonbonded list, or to a DihedralSingletonScore to score each dihedral.
296 
297  If dihedrals are explicitly listed in the CHARMM topology file, they
298  can be created if desired by calling CHARMMTopology::add_dihedrals()
299  rather than this function.
300 
301  \return a list of the newly-created Dihedral particles.
302 
303  \see CHARMMTopology::add_bonds().
304  */
305  kernel::Particles create_dihedrals(kernel::Particles bonds) const;
306 
308 
309  private:
310 
311  virtual String get_force_field_atom_type(Atom atom) const;
312 
313  void read_parameter_file(base::TextInput input_file);
314  // read topology file
315  void read_topology_file(base::TextInput input_file,
316  bool translate_names_to_pdb);
317 
318  void add_angle(kernel::Particle *p1, kernel::Particle *p2,
319  kernel::Particle *p3, kernel::Particles &ps) const;
320  void add_dihedral(kernel::Particle *p1, kernel::Particle *p2,
322  kernel::Particles &ps) const;
323 
324  ResidueType parse_residue_line(const String &line,
325  bool translate_names_to_pdb);
326  void parse_atom_line(const String &line, ResidueType curr_res_type,
327  CHARMMResidueTopologyBase *residue,
328  bool translate_names_to_pdb);
329  void parse_bond_line(const String &line, ResidueType curr_res_type,
330  CHARMMResidueTopologyBase *residue,
331  bool translate_names_to_pdb);
332 
333  void parse_nonbonded_parameters_line(String line);
334  void parse_bonds_parameters_line(String line);
335  void parse_angles_parameters_line(String line);
336  void parse_dihedrals_parameters_line(String line, DihedralParameters &param);
337  base::WarningContext warn_context_;
338 };
339 
341 
342 /** The default CHARMM parameters support normal amino acid
343  and nucleic acid residues and the atoms found in them.
344  To use CHARMM with heterogens or non-standard residues,
345  a different CHARMM parameters file must be used.
346 
347  No hydrogen parameters are read.
348 
349  \see get_all_atom_CHARMM_parameters();
350 */
352 
353 /** The default CHARMM parameters support normal amino acid
354  and nucleic acid residues and the atoms found in them.
355  To use CHARMM with heterogens or non-standard residues,
356  a different CHARMM parameters file must be used.
357 
358  \see get_heavy_atom_CHARMM_parameters()
359 */
361 
362 IMPATOM_END_NAMESPACE
363 
364 #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