I am trying to compute the residue-residue pairwise energy for several residues inside a protein structure. I found the method scorefxn.eval_ci_2b() can compute the pairwise energy between residues for a given scorefxn. This works for the LJ, Solvation and electrostatic potentials. Also as epxected, it ignores intra-residue interactions and backbone parameters when computing with eval_ci_2b(). However, I was wondering why the Hydrogen bonded potential is also not included? Is there an easy way to include the hydrogen bond energy in the pairwise energy?
Not sure if this is helpful to even see, but I have a minimum example of what I am trying to do below. Not sure if it would run as it was cut/pasted from various sections of my code, but it should get at the gist of what I was trying to do.
import numpy as np import pyrosetta as pyr from pyrosetta import Pose from pyrosetta.toolbox import get_hbonds from pyrosetta import teaching as pyrt pose = pyr.pose_from_pdb(my_pdb_file) stuff = "(fa_atr 1) (fa_rep 0.55) (fa_sol 1) (fa_intra_rep 0.005) (fa_intra_sol_xover4 1) (lk_ball_wtd 1) (fa_elec 1) (pro_close 1.25) (hbond_sr_bb 1) (hbond_lr_bb 1) (hbond_bb_sc 1) (hbond_sc 1) (dslf_fa13 1.25) (omega 0.4) (fa_dun 0.7) (p_aa_pp 0.6) (yhh_planarity 0.625) (ref 1) (rama_prepro 0.45)" things = stuff.strip().split(")") weights =  order =  for i_thing_count,thing in enumerate(things): if len(thing) > 0: if i_thing_count == 0: useful_stuff = thing[1:].split() else: useful_stuff = thing[2:].split() order.append(eval("pyrt." + useful_stuff)) weights.append(float(useful_stuff)) assert len(weights) == len(order) scorefxn = pyrt.ScoreFunction() for thing, wt in zip(order,weights): scorefxn.set_weight(thing, wt) hbond_set = get_hbonds(pose) pair_E = np.zeros((nresidues,nresidues)) for contact in these_contacts: idx = contact jdx = contact emap = pyrt.EMapVector() scorefxn.eval_ci_2b(pose.residue(idx), pose.residue(jdx), pose, emap) this_E = 0. for thing,wt in zip(order,weights): if idx == 8 and jdx == 20: print thing print emap[thing] this_E += emap[thing] * wt pair_E[idx-1, jdx-1] = this_E