IMP logo
IMP Reference Guide  2.10.0
The Integrative Modeling Platform
charmm_forcefield.py
1 ## \example atom/charmm_forcefield.py
2 # In this example, a PDB file is read in and scored using the CHARMM
3 # forcefield. For more control over the setup of the forcefield, see
4 # the 'charmm_forcefield_verbose.py' example.
5 
6 from __future__ import print_function
7 import IMP.atom
8 import IMP.container
9 import sys
10 
11 IMP.setup_from_argv(sys.argv, "CHARMM forcefield")
12 
13 # Create an IMP model and add a heavy atom-only protein from a PDB file
14 m = IMP.Model()
15 prot = IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m,
17 
18 # Read in the CHARMM heavy atom topology and parameter files
20 
21 # Using the CHARMM libraries, determine the ideal topology (atoms and their
22 # connectivity) for the PDB file's primary sequence
23 topology = ff.create_topology(prot)
24 
25 # Typically this modifies the C and N termini of each chain in the protein by
26 # applying the CHARMM CTER and NTER patches. Patches can also be manually
27 # applied at this point, e.g. to add disulfide bridges.
28 topology.apply_default_patches()
29 
30 # Make the PDB file conform with the topology; i.e. if it contains extra
31 # atoms that are not in the CHARMM topology file, remove them; if it is
32 # missing atoms (e.g. sidechains, hydrogens) that are in the CHARMM topology,
33 # add them and construct their Cartesian coordinates from internal coordinate
34 # information.
35 topology.setup_hierarchy(prot)
36 
37 # Set up and evaluate the stereochemical part (bonds, angles, dihedrals,
38 # impropers) of the CHARMM forcefield
40 
41 # Add non-bonded interaction (in this case, Lennard-Jones). This needs to
42 # know the radii and well depths for each atom, so add them from the forcefield
43 # (they can also be assigned manually using the XYZR or LennardJones
44 # decorators):
45 ff.add_radii(prot)
46 ff.add_well_depths(prot)
47 
48 # Get a list of all atoms in the protein, and put it in a container
49 atoms = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
51 
52 # Add a restraint for the Lennard-Jones interaction. This is built from
53 # a collection of building blocks. First, a ClosePairContainer maintains a list
54 # of all pairs of Particles that are close. Next, all 1-2, 1-3 and 1-4 pairs
55 # from the stereochemistry created above are filtered out.
56 # Then, a LennardJonesPairScore scores a pair of atoms with the Lennard-Jones
57 # potential. Finally, a PairsRestraint is used which simply applies the
58 # LennardJonesPairScore to each pair in the ClosePairContainer.
59 nbl = IMP.container.ClosePairContainer(cont, 4.0)
60 nbl.add_pair_filter(r.get_pair_filter())
61 
62 sf = IMP.atom.ForceSwitch(6.0, 7.0)
64 restraints = [r, IMP.container.PairsRestraint(ps, nbl)]
65 scoring_function = IMP.core.RestraintsScoringFunction(restraints)
66 
67 # it gets awfully slow with internal checks
68 IMP.set_check_level(IMP.USAGE)
69 # Finally, evaluate the score of the whole system (without derivatives)
70 print(scoring_function.evaluate(False))
Select non water and non hydrogen atoms.
Definition: pdb.h:245
CHARMMParameters * get_heavy_atom_CHARMM_parameters()
Enforce CHARMM stereochemistry on the given Hierarchy.
Strings setup_from_argv(const Strings &argv, std::string description, std::string positional_description, int num_positional)
Various classes to hold sets of particles.
std::string get_example_path(std::string file_name)
Return the full path to one of this module's example files.
Return all close unordered pairs of particles taken from the SingletonContainer.
void read_pdb(TextInput input, int model, Hierarchy h)
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Lennard-Jones score between a pair of particles.
Store a list of ParticleIndexes.
Smooth interaction scores by switching the derivatives (force switch).
Functionality for loading, creating, manipulating and scoring atomic structures.
Applies a PairScore to each Pair in a list.
void set_check_level(CheckLevel tf)
Control runtime checks in the code.
Definition: exception.h:73