1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
|
from ost import io, geom
from promod3 import loop
# setup system creator
ff_lookup = loop.ForcefieldLookup.GetDefault()
mm_sys = loop.MmSystemCreator(ff_lookup)
# load example (has res. numbering starting at 1)
ent = io.LoadPDB("data/1CRN.pdb")
res_list = ent.residues
num_residues = len(res_list)
# get all atom positions for full protein
all_atoms = loop.AllAtomPositions(res_list)
# here full structure in res_indices but in practice this could
# be just a subset of residues relevant to the loop of interest
res_indices = list(range(num_residues))
# define two loops (indices into res_indices list)
loop_start_indices = [10, 20]
loop_lengths = [6, 4]
# define which of the res_indices is terminal
is_n_ter = [True] + [False] * (num_residues - 1)
is_c_ter = [False] * (num_residues - 1) + [True]
# get disulfid bridges
disulfid_bridges = mm_sys.GetDisulfidBridges(all_atoms, res_indices)
# setup MM system
mm_sys.SetupSystem(all_atoms, res_indices, loop_start_indices,
loop_lengths, is_n_ter, is_c_ter, disulfid_bridges)
# run simulation
sim = mm_sys.GetSimulation()
print("Potential energy before: %g" % sim.GetPotentialEnergy())
sim.ApplySD(0.01, 100)
print("Potential energy after: %g" % sim.GetPotentialEnergy())
# extract new loop positions and store it
mm_sys.ExtractLoopPositions(all_atoms, res_indices)
io.SavePDB(all_atoms.ToEntity(), "mm_sys_output.pdb")
|
