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

modeller_test_energy2.py

import modeller

modeller.log.none()

env = modeller.environ()

env.io.hetatm = True
#soft sphere potential
env.edat.dynamic_sphere=False
env.edat.nonbonded_sel_atoms = 2
#lennard-jones potential (more accurate)
env.edat.dynamic_lennard=True
env.edat.contact_shell = 4.0
env.edat.update_dynamic = 0.39
# 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')

import Bio
from Bio.PDB import PDBList
'''Selecting structures from PDB'''
pdbl = PDBList()
PDBlist2=['2a3d']
for i in PDBlist2:
    pdbl.retrieve_pdb_file(i,pdir='PDB')
    
import Bio.Seq
import Bio.SeqIO
import Bio.SeqUtils
parser = Bio.PDB.PDBParser()
for i in PDBlist2:
    struc = parser.get_structure(i,"PDB/pdb" + i +".ent")
    seq = Bio.Seq.Seq(''.join([Bio.SeqUtils.seq1(x.get_resname()) for x in \
                struc.get_residues()]))
    Bio.SeqIO.write(Bio.SeqRecord.SeqRecord(seq, id="PDB/pdb" + i +".ent"), \
                    open("fasta/%s.fasta"%i, "w"), "fasta")
    
mdl = modeller.model(env, file='PDB/pdb2a3d.ent')
aln = modeller.alignment(env)
aln.append_model(mdl, align_codes='fasta/2a3d.fasta', atom_files="PDB/pdb" + i +".ent")
aln.append_model(mdl, align_codes='fasta/2a3d.fasta', atom_files="PDB/pdb" + i +".ent")

mdl.clear_topology()
mdl.generate_topology(aln['fasta/2a3d.fasta'])
mdl.transfer_xyz(aln)
mdl.build(initialize_xyz=False, build_method='INTERNAL_COORDINATES')

s0 = modeller.selection(mdl.residue_range('%s'%str(2),'%s'%str(20)))
s1 = modeller.selection(mdl.residue_range('%s'%str(24),'%s'%str(43)))
s0s1 = modeller.selection(mdl.residue_range('%s'%str(2),'%s'%str(20)), \
                    mdl.residue_range('%s'%str(24),'%s'%str(43)))

interface_energy = s0s1.energy()[0] - (s0.energy()[0] + s1.energy()[0])