IMP  2.0.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-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/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
35  << "; ideal value: " << ideal;});
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
49  << "; ideal value: " << ideal;});
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, Pointer<CHARMMIdealResidueTopology> >
74  residue_topologies_;
75  std::map<std::string, 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
85  find_dihedral(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(std::make_pair(ResidueType(res->get_type()),
124  res));
125  }
126 
127  CHARMMIdealResidueTopology *get_residue_topology(ResidueType type) const
128  {
129  std::map<ResidueType, Pointer<CHARMMIdealResidueTopology> >::const_iterator
130  it = 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, 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 "
185  << type1 << "-" << type2, 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 = internal::CHARMMAngleNames(type1, type2,
197  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 "
202  << type1 << "-" << type2 << "-" << type3, 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  */
217  CHARMMDihedralParametersList get_dihedral_parameters(
218  std::string type1, std::string type2, std::string type3,
219  std::string type4) const {
221  internal::CHARMMDihedralNames types = internal::CHARMMDihedralNames(
222  type1, type2, type3, type4);
223  // Get the first match, using wildcards
224  DihedralParameters::const_iterator match =
225  find_dihedral(dihedral_parameters_.begin(),
226  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  */
254  const CHARMMDihedralParameters &get_improper_parameters(
255  std::string type1, std::string type2, std::string type3,
256  std::string type4) const {
257  internal::CHARMMDihedralNames types = internal::CHARMMDihedralNames(
258  type1, type2, type3, type4);
259  // Return just the first match; wildcards are OK
260  DihedralParameters::const_iterator it =
261  find_dihedral(improper_parameters_.begin(),
262  improper_parameters_.end(), types, true);
263  if (it != improper_parameters_.end()) {
264  return it->second;
265  } else {
266  IMP_THROW("No CHARMM parameters found for improper "
267  << type1 << "-" << type2 << "-" << type3 << "-" << type4,
268  IndexException);
269  }
270  }
271 
272  //! Auto-generate Angle particles from the passed list of Bond particles.
273  /** The angles consist of all unique pairs of bonds which share an
274  endpoint. If no parameters are found for an angle, it is simply
275  created without those parameters.
276 
277  The list of newly-created Angle particles can be passed to a
278  StereochemistryPairFilter to exclude 1-3 interactions from the
279  nonbonded list, or to an AngleSingletonScore to score each angle.
280 
281  \return a list of the newly-created Angle particles.
282 
283  \see CHARMMTopology::add_bonds().
284  */
285  Particles create_angles(Particles bonds) const;
286 
287  //! Auto-generate Dihedral particles from the passed list of Bond particles.
288  /** The dihedrals consist of all unique triples of bonds which form
289  dihedrals. If no parameters are found for a dihedral, it is simply
290  created without those parameters; if multiple sets of parameters are
291  found, multiple copies of the dihedral are created, each with one set
292  of parameters.
293 
294  The list of newly-created Dihedral particles can be passed to a
295  StereochemistryPairFilter to exclude 1-4 interactions from the
296  nonbonded list, or to a DihedralSingletonScore to score each dihedral.
297 
298  If dihedrals are explicitly listed in the CHARMM topology file, they
299  can be created if desired by calling CHARMMTopology::add_dihedrals()
300  rather than this function.
301 
302  \return a list of the newly-created Dihedral particles.
303 
304  \see CHARMMTopology::add_bonds().
305  */
306  Particles create_dihedrals(Particles bonds) const;
307 
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(Particle *p1, Particle *p2, Particle *p3, Particles &ps) const;
319  void add_dihedral(Particle *p1, Particle *p2, Particle *p3, Particle *p4,
320  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,
335  DihedralParameters &param);
336  WarningContext warn_context_;
337 };
338 
340 
341 
342 
343 /** The default CHARMM parameters support normal amino acid
344  and nucleic acid residues and the atoms found in them.
345  To use CHARMM with heterogens or non-standard residues,
346  a different CHARMM parameters file must be used.
347 
348  No hydrogen parameters are read.
349 
350  \see get_all_atom_CHARMM_parameters();
351 */
353 
354 
355 /** The default CHARMM parameters support normal amino acid
356  and nucleic acid residues and the atoms found in them.
357  To use CHARMM with heterogens or non-standard residues,
358  a different CHARMM parameters file must be used.
359 
360  \see get_heavy_atom_CHARMM_parameters()
361 */
363 
364 IMPATOM_END_NAMESPACE
365 
366 #endif /* IMPATOM_CHARMM_PARAMETERS_H */