IMP logo
IMP Reference Guide  2.21.0
The Integrative Modeling Platform
atom/charmm_forcefield_verbose.py

In this example, a PDB file is read in and scored using the CHARMM forcefield. It is similar to the 'charmm_forcefield.py' example, but fully works through each step of the procedure using lower-level IMP classes. This is useful if you want to customize the way in which the forcefield is applied.

1 ## \example atom/charmm_forcefield_verbose.py
2 # In this example, a PDB file is read in and scored using the CHARMM
3 # forcefield. It is similar to the 'charmm_forcefield.py' example, but fully
4 # works through each step of the procedure using lower-level IMP classes.
5 # This is useful if you want to customize the way in which the forcefield
6 # is applied.
7 #
8 
9 from __future__ import print_function
10 import IMP.atom
11 import IMP.container
12 import sys
13 
14 IMP.setup_from_argv(sys.argv, "CHARMM forcefield verbose")
15 
16 # Create an IMP model and add a heavy atom-only protein from a PDB file
17 m = IMP.Model()
18 prot = IMP.atom.read_pdb(IMP.atom.get_example_path("example_protein.pdb"), m,
20 
21 # Read in the CHARMM heavy atom topology and parameter files
23 
24 # Using the CHARMM libraries, determine the ideal topology (atoms and their
25 # connectivity) for the PDB file's primary sequence
26 topology = ff.create_topology(prot)
27 
28 # Typically this modifies the C and N termini of each chain in the protein by
29 # applying the CHARMM CTER and NTER patches. Patches can also be manually
30 # applied at this point, e.g. to add disulfide bridges.
31 topology.apply_default_patches()
32 
33 # Each atom is mapped to its CHARMM type. These are needed to look up bond
34 # lengths, Lennard-Jones radii etc. in the CHARMM parameter file. Atom types
35 # can also be manually assigned at this point using the CHARMMAtom decorator.
36 topology.add_atom_types(prot)
37 
38 # Remove any atoms that are in the PDB file but not in the topology, and add
39 # in any that are in the topology but not the PDB.
41 topology.add_missing_atoms(prot)
42 
43 # Construct Cartesian coordinates for any atoms that were added
44 topology.add_coordinates(prot)
45 
46 # Generate and return lists of bonds, angles, dihedrals and impropers for
47 # the protein. Each is a Particle in the model which defines the 2, 3 or 4
48 # atoms that are bonded, and adds parameters such as ideal bond length
49 # and force constant. Note that bonds and impropers are explicitly listed
50 # in the CHARMM topology file, while angles and dihedrals are generated
51 # automatically from an existing set of bonds. These particles only define the
52 # bonds, but do not score them or exclude them from the nonbonded list.
53 bonds = topology.add_bonds(prot)
54 angles = ff.create_angles(bonds)
55 dihedrals = ff.create_dihedrals(bonds)
56 impropers = topology.add_impropers(prot)
57 
58 # Maintain stereochemistry by scoring bonds, angles, dihedrals and impropers
59 
60 # Score all of the bonds. This is done by combining IMP 'building blocks':
61 # - A ListSingletonContainer simply manages a list of the bond particles.
62 # - A BondSingletonScore, when given a bond particle, scores the bond by
63 # calculating the distance between the two atoms it bonds, subtracting the
64 # ideal value, and weighting the result by the bond's "stiffness", such that
65 # an "ideal" bond scores zero, and bonds away from equilibrium score
66 # non-zero. It then hands off to a UnaryFunction to actually penalize the
67 # value. In this case, a Harmonic UnaryFunction is used with a mean of zero,
68 # so that bond lengths are harmonically restrained.
69 # - A SingletonsRestraint simply goes through each of the bonds in the
70 # container and scores each one in turn.
71 cont = IMP.container.ListSingletonContainer(m, bonds, "bonds")
73 r = IMP.container.SingletonsRestraint(bss, cont, "bonds")
74 rs = [r]
75 
76 # Score angles, dihedrals, and impropers. In the CHARMM forcefield, angles and
77 # impropers are harmonically restrained, so this is the same as for bonds.
78 # Dihedrals are scored internally by a periodic (cosine) function.
79 cont = IMP.container.ListSingletonContainer(m, angles, "angles")
81 r = IMP.container.SingletonsRestraint(bss, cont, "angles")
82 rs.append(r)
83 
84 cont = IMP.container.ListSingletonContainer(m, dihedrals, "dihedrals")
86 r = IMP.container.SingletonsRestraint(bss, cont, "dihedrals")
87 rs.append(r)
88 
89 cont = IMP.container.ListSingletonContainer(m, impropers, "impropers")
91 rs.append(IMP.container.SingletonsRestraint(bss, cont, "improppers"))
92 
93 # Add non-bonded interaction (in this case, Lennard-Jones). This needs to
94 # know the radii and well depths for each atom, so add them from the forcefield
95 # (they can also be assigned manually using the XYZR or LennardJones
96 # decorators):
97 ff.add_radii(prot)
98 ff.add_well_depths(prot)
99 
100 # Get a list of all atoms in the protein, and put it in a container
101 atoms = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
103 
104 # Add a restraint for the Lennard-Jones interaction. Again, this is built from
105 # a collection of building blocks. First, a ClosePairContainer maintains a list
106 # of all pairs of Particles that are close. A StereochemistryPairFilter is used
107 # to exclude atoms from this list that are bonded to each other or are involved
108 # in an angle or dihedral (1-3 or 1-4 interaction). Then, a
109 # LennardJonesPairScore scores a pair of atoms with the
110 # Lennard-Jones potential.
111 # Finally, a PairsRestraint is used which simply applies the
112 # LennardJonesPairScore to each pair in the ClosePairContainer.
113 nbl = IMP.container.ClosePairContainer(cont, 4.0)
115 pair_filter.set_bonds(bonds)
116 pair_filter.set_angles(angles)
117 pair_filter.set_dihedrals(dihedrals)
118 nbl.add_pair_filter(pair_filter)
119 
120 sf = IMP.atom.ForceSwitch(6.0, 7.0)
122 rs.append(IMP.container.PairsRestraint(ps, nbl))
123 
124 score_func = IMP.core.RestraintsScoringFunction(rs)
125 
126 # it gets awfully slow with internal checks
127 IMP.set_check_level(IMP.USAGE)
128 
129 # Finally, evaluate the score of the whole system (without derivatives)
130 print(score_func.evaluate(False))