Test for energy estimation of interfaces within Modeller (https://salilab.org/modeller)

modeller_test_energy.py

import modeller
import Bio
from Bio.PDB import PDBList

# modeller stuff
modeller.log.none()
env = modeller.environ()
env.edat.nonbonded_sel_atoms = 2
# Read customized topology file with phosphoserines (or standard one)
env.libs.topology.read(file='$(LIB)/top_heav.lib')
# Read customized CHARMM parameter library with phosphoserines (or standard one)
env.libs.parameters.read(file='$(LIB)/par.lib')

# Selecting structures from PDB
pdbl = PDBList()
PDBlist2=['2a3d']
for i in PDBlist2:
    pdbl.retrieve_pdb_file(i,pdir='PDB')

# Loading as Modeller model
mdl = modeller.model(env, file='PDB/pdb2a3d.ent')

# Creating selections for helix 1 and helix 2, and both helix 1 and 2
s0 = modeller.selection(mdl.residue_range('%s:A'%str(2),'%s:A'%str(20)))
s1 = modeller.selection(mdl.residue_range('%s:A'%str(24),'%s:A'%str(43)))
s0s1 = modeller.selection(mdl.residue_range('%s:A'%str(2),'%s:A'%str(20)), \
                    mdl.residue_range('%s:A'%str(24),'%s:A'%str(43)))

# Interface energy as difference between s0s1 energy and s0 + s1 energies
s0s1.assess_dope() - (s0.assess_dope() + s1.assess_dope())