wepy_glue_utilities/ppi_scale.py at main · ADicksonLab/wepy_glue_utilities · GitHub

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
"""
    This is an example script on how to scale down the PPI between two proteins.

    Seref Berk Atik (atiksere@msu.edu)
    Dickson Lab
    Department of Computational Mathematics, Science & Engineering
    Michigan State University

    December 2025

    References:
        Jo, S., Kim, T., Iyer, V. G., & Im, W. (2008). CHARMM-GUI: a web-based graphical user interface for CHARMM. Journal of Computational Chemistry, 29(11), 1859–1865. doi:10.1002/jcc.20945
        Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., … Im, W. (2016). CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. Journal of Chemical Theory and Computation, 12(1), 405–413. doi:10.1021/acs.jctc.5b00935
        Lee, J., Hitzenberger, M., Rieger, M., Kern, N. R., Zacharias, M., & Im, W. (2020). CHARMM-GUI supports the Amber force fields. The Journal of Chemical Physics, 153(3), 035103. doi:10.1063/5.0012280
"""

import pickle as pkl

import sys
import mdtraj as mdj

import openmm as omm

if __name__ == "__main__":
    sys_tag = sys.argv[1]
    ppi_scale_fac = sys.argv[2]

    pdb_path = "path_to_pdb"
    system_pkl_path = "path_to_system_pkl"

    with open(system_pkl_path, 'rb') as f:
        system = pkl.load(f)

    pdb = mdj.load_pdb(pdb_path)
    prot1_idxs = pdb.top.select("protein_1 selection string")
    prot2_idxs = pdb.top.select("protein_2 selection string")

    forces = system.getForces()

    old_nb = forces[6]
    assert isinstance(old_nb, omm.NonbondedForce)
    old_cnb = forces[7]
    assert isinstance(old_cnb, omm.CustomNonbondedForce)

    if isinstance(forces[9], omm.CustomNonbondedForce):
        system.removeForce(9)


    # NOTE: Might be subject to change, double check your forcefield strings prior to running and update the script accordingly.
    efunc = "-escale * (138.935456 * q1 * q2 / r + select(step(r-Ron), (cr12*rjunk12 - cr6*rjunk6), (ccnba/r^12-ccnbb/r^6-ccnba/onoff^6+ccnbb/onoff^3))); cr12 = ccnba*ofdif6; cr6  = ccnbb*ofdif3; rjunk12 = (1.0/r^6-1.0/Roff6)^2; rjunk6 = (1.0/r^3-1.0/Roff3)^2; ccnbb = bcoef(type1, type2); ccnba = acoef(type1, type2)^2; ofdif6 = Roff6/(Roff6 - Ron^6); ofdif3 = Roff3/(Roff3 - Ron^3); onoff = Roff * Ron; Roff6 = Roff^6; Roff3 = Roff^3; Ron  = 1.000000; Roff = 1.200000;"

    new_cnb = omm.CustomNonbondedForce(efunc)

    new_cnb.addGlobalParameter('escale',1-float(ppi_scale_fac))

    new_cnb.setCutoffDistance(old_cnb.getCutoffDistance())
    new_cnb.setForceGroup(old_cnb.getForceGroup())
    new_cnb.setNonbondedMethod(old_cnb.getNonbondedMethod())

    for i in range(old_cnb.getNumExclusions()):
        new_cnb.addExclusion(*old_cnb.getExclusionParticles(i))

    for i in range(old_cnb.getNumTabulatedFunctions()):
        tab_name = old_cnb.getTabulatedFunctionName(i)
        new_tf = old_cnb.getTabulatedFunction(i).Copy()
        new_tf.thisown = True
        new_cnb.addTabulatedFunction(tab_name,new_tf)

    new_cnb.addPerParticleParameter('type')
    new_cnb.addPerParticleParameter('q')

    new_cnb.addInteractionGroup(set(prot1_idxs),set(prot2_idxs))

    system.addForce(new_cnb)

    for force in system.getForces():
        print(force)

    pkl.dump("scaled_system_pk_path",'wb')