Break propka.calculations into parts.
Broken into three files to simplify import analysis: * `calculations.py` - distance calculations and simple math * `energy.py` - energy calculations * `hydrogens.py` - hydrogen geometry calculations Working towards addressing https://github.com/jensengroup/propka-3.1/issues/49
This commit is contained in:
Nathan Baker
@@ -66,880 +66,3 @@ def get_smallest_distance(atoms1, atoms2):
|
|||||||
res_atom1 = atom1
|
res_atom1 = atom1
|
||||||
res_atom2 = atom2
|
res_atom2 = atom2
|
||||||
return [res_atom1, math.sqrt(res_dist), res_atom2]
|
return [res_atom1, math.sqrt(res_dist), res_atom2]
|
||||||
|
|
||||||
|
|
||||||
# TODO - the next set of functions form a distinct "module" for hydrogen addition
|
|
||||||
|
|
||||||
|
|
||||||
def setup_bonding_and_protonation(parameters, molecular_container):
|
|
||||||
"""Set up bonding and protonation for a molecule.
|
|
||||||
|
|
||||||
NOTE - the unused `parameters` argument is required for compatibility in version.py
|
|
||||||
TODO - figure out why there is a similar function in version.py
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: not used
|
|
||||||
molecular_container: molecule container.
|
|
||||||
"""
|
|
||||||
# make bonds
|
|
||||||
my_bond_maker = setup_bonding(molecular_container)
|
|
||||||
# set up ligand atom names
|
|
||||||
set_ligand_atom_names(molecular_container)
|
|
||||||
# apply information on pi electrons
|
|
||||||
my_bond_maker.add_pi_electron_information(molecular_container)
|
|
||||||
# Protonate atoms
|
|
||||||
if molecular_container.options.protonate_all:
|
|
||||||
PROTONATOR = propka.protonate.Protonate(verbose=False)
|
|
||||||
PROTONATOR.protonate(molecular_container)
|
|
||||||
|
|
||||||
|
|
||||||
def setup_bonding(molecular_container):
|
|
||||||
"""Set up bonding for a molecular container.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
molecular_container: the molecular container in question
|
|
||||||
Returns:
|
|
||||||
BondMaker object
|
|
||||||
"""
|
|
||||||
my_bond_maker = propka.bonds.BondMaker()
|
|
||||||
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
|
|
||||||
return my_bond_maker
|
|
||||||
|
|
||||||
|
|
||||||
def setup_bonding_and_protonation_30_style(parameters, molecular_container):
|
|
||||||
"""Set up bonding for a molecular container.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: parameters for calculation
|
|
||||||
molecular_container: the molecular container in question
|
|
||||||
Returns:
|
|
||||||
BondMaker object
|
|
||||||
"""
|
|
||||||
# Protonate atoms
|
|
||||||
protonate_30_style(molecular_container)
|
|
||||||
# make bonds
|
|
||||||
bond_maker = propka.bonds.BondMaker()
|
|
||||||
bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
|
|
||||||
return bond_maker
|
|
||||||
|
|
||||||
|
|
||||||
def protonate_30_style(molecular_container):
|
|
||||||
"""Protonate the molecule.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
molecular_container: molecule
|
|
||||||
"""
|
|
||||||
for name in molecular_container.conformation_names:
|
|
||||||
info('Now protonating', name)
|
|
||||||
# split atom into residues
|
|
||||||
curres = -1000000
|
|
||||||
residue = []
|
|
||||||
o_atom = None
|
|
||||||
c_atom = None
|
|
||||||
for atom in molecular_container.conformations[name].atoms:
|
|
||||||
if atom.res_num != curres:
|
|
||||||
curres = atom.res_num
|
|
||||||
if len(residue) > 0:
|
|
||||||
#backbone
|
|
||||||
[o_atom, c_atom] = add_backbone_hydrogen(
|
|
||||||
residue, o_atom, c_atom)
|
|
||||||
#arginine
|
|
||||||
if residue[0].res_name == 'ARG':
|
|
||||||
add_arg_hydrogen(residue)
|
|
||||||
#histidine
|
|
||||||
if residue[0].res_name == 'HIS':
|
|
||||||
add_his_hydrogen(residue)
|
|
||||||
#tryptophan
|
|
||||||
if residue[0].res_name == 'TRP':
|
|
||||||
add_trp_hydrogen(residue)
|
|
||||||
#amides
|
|
||||||
if residue[0].res_name in ['GLN', 'ASN']:
|
|
||||||
add_amd_hydrogen(residue)
|
|
||||||
residue = []
|
|
||||||
if atom.type == 'atom':
|
|
||||||
residue.append(atom)
|
|
||||||
|
|
||||||
|
|
||||||
def set_ligand_atom_names(molecular_container):
|
|
||||||
"""Set names for ligands in molecular container.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
molecular_container: molecular container for ligand names
|
|
||||||
"""
|
|
||||||
for name in molecular_container.conformation_names:
|
|
||||||
molecular_container.conformations[name].set_ligand_atom_names()
|
|
||||||
|
|
||||||
|
|
||||||
def add_arg_hydrogen(residue):
|
|
||||||
"""Adds Arg hydrogen atoms to residues according to the 'old way'.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
residue: arginine residue to protonate
|
|
||||||
Returns:
|
|
||||||
list of hydrogen atoms
|
|
||||||
"""
|
|
||||||
#info('Adding arg H',residue)
|
|
||||||
for atom in residue:
|
|
||||||
if atom.name == "CD":
|
|
||||||
cd_atom = atom
|
|
||||||
elif atom.name == "CZ":
|
|
||||||
cz_atom = atom
|
|
||||||
elif atom.name == "NE":
|
|
||||||
ne_atom = atom
|
|
||||||
elif atom.name == "NH1":
|
|
||||||
nh1_atom = atom
|
|
||||||
elif atom.name == "NH2":
|
|
||||||
nh2_atom = atom
|
|
||||||
h1_atom = protonate_sp2(cd_atom, ne_atom, cz_atom)
|
|
||||||
h1_atom.name = "HE"
|
|
||||||
h2_atom = protonate_direction(nh1_atom, ne_atom, cz_atom)
|
|
||||||
h2_atom.name = "HN1"
|
|
||||||
h3_atom = protonate_direction(nh1_atom, ne_atom, cd_atom)
|
|
||||||
h3_atom.name = "HN2"
|
|
||||||
h4_atom = protonate_direction(nh2_atom, ne_atom, cz_atom)
|
|
||||||
h4_atom.name = "HN3"
|
|
||||||
h5_atom = protonate_direction(nh2_atom, ne_atom, h1_atom)
|
|
||||||
h5_atom.name = "HN4"
|
|
||||||
return [h1_atom, h2_atom, h3_atom, h4_atom, h5_atom]
|
|
||||||
|
|
||||||
|
|
||||||
def add_his_hydrogen(residue):
|
|
||||||
"""Adds His hydrogen atoms to residues according to the 'old way'.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
residue: histidine residue to protonate
|
|
||||||
"""
|
|
||||||
for atom in residue:
|
|
||||||
if atom.name == "CG":
|
|
||||||
cg_atom = atom
|
|
||||||
elif atom.name == "ND1":
|
|
||||||
nd_atom = atom
|
|
||||||
elif atom.name == "CD2":
|
|
||||||
cd_atom = atom
|
|
||||||
elif atom.name == "CE1":
|
|
||||||
ce_atom = atom
|
|
||||||
elif atom.name == "NE2":
|
|
||||||
ne_atom = atom
|
|
||||||
hd_atom = protonate_sp2(cg_atom, nd_atom, ce_atom)
|
|
||||||
hd_atom.name = "HND"
|
|
||||||
he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom)
|
|
||||||
he_atom.name = "HNE"
|
|
||||||
|
|
||||||
|
|
||||||
def add_trp_hydrogen(residue):
|
|
||||||
"""Adds Trp hydrogen atoms to residues according to the 'old way'.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
residue: tryptophan residue to protonate
|
|
||||||
"""
|
|
||||||
cd_atom = None
|
|
||||||
ne_atom = None
|
|
||||||
for atom in residue:
|
|
||||||
if atom.name == "CD1":
|
|
||||||
cd_atom = atom
|
|
||||||
elif atom.name == "NE1":
|
|
||||||
ne_atom = atom
|
|
||||||
elif atom.name == "CE2":
|
|
||||||
ce_atom = atom
|
|
||||||
if (cd_atom is None) or (ne_atom is None) or (ce_atom is None):
|
|
||||||
str_ = "Unable to find all atoms for {0:s} {1:s}".format(
|
|
||||||
residue[0].res_name, residue[0].res_num)
|
|
||||||
raise ValueError(str_)
|
|
||||||
he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom)
|
|
||||||
he_atom.name = "HNE"
|
|
||||||
|
|
||||||
|
|
||||||
def add_amd_hydrogen(residue):
|
|
||||||
"""Adds Gln & Asn hydrogen atoms to residues according to the 'old way'.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
residue: glutamine or asparagine residue to protonate
|
|
||||||
"""
|
|
||||||
c_atom = None
|
|
||||||
o_atom = None
|
|
||||||
n_atom = None
|
|
||||||
for atom in residue:
|
|
||||||
if ((atom.res_name == "GLN" and atom.name == "CD")
|
|
||||||
or (atom.res_name == "ASN" and atom.name == "CG")):
|
|
||||||
c_atom = atom
|
|
||||||
elif ((atom.res_name == "GLN" and atom.name == "OE1")
|
|
||||||
or (atom.res_name == "ASN" and atom.name == "OD1")):
|
|
||||||
o_atom = atom
|
|
||||||
elif ((atom.res_name == "GLN" and atom.name == "NE2")
|
|
||||||
or (atom.res_name == "ASN" and atom.name == "ND2")):
|
|
||||||
n_atom = atom
|
|
||||||
if (c_atom is None) or (o_atom is None) or (n_atom is None):
|
|
||||||
str_ = "Unable to find all atoms for {0:s} {1:s}".format(
|
|
||||||
residue[0].res_name, residue[0].res_num)
|
|
||||||
raise ValueError(str_)
|
|
||||||
h1_atom = protonate_direction(n_atom, o_atom, c_atom)
|
|
||||||
h1_atom.name = "HN1"
|
|
||||||
h2_atom = protonate_average_direction(n_atom, c_atom, o_atom)
|
|
||||||
h2_atom.name = "HN2"
|
|
||||||
|
|
||||||
|
|
||||||
def add_backbone_hydrogen(residue, o_atom, c_atom):
|
|
||||||
"""Adds hydrogen backbone atoms to residues according to the old way.
|
|
||||||
|
|
||||||
dR is wrong for the N-terminus (i.e. first residue) but it doesn't affect
|
|
||||||
anything at the moment. Could be improved, but works for now.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
residue: residue to protonate
|
|
||||||
o_atom: backbone oxygen atom
|
|
||||||
c_atom: backbone carbon atom
|
|
||||||
Returns:
|
|
||||||
[new backbone oxygen atom, new backbone carbon atom]
|
|
||||||
"""
|
|
||||||
new_c_atom = None
|
|
||||||
new_o_atom = None
|
|
||||||
n_atom = None
|
|
||||||
for atom in residue:
|
|
||||||
if atom.name == "N":
|
|
||||||
n_atom = atom
|
|
||||||
if atom.name == "C":
|
|
||||||
new_c_atom = atom
|
|
||||||
if atom.name == "O":
|
|
||||||
new_o_atom = atom
|
|
||||||
if None in [c_atom, o_atom, n_atom]:
|
|
||||||
return [new_o_atom, new_c_atom]
|
|
||||||
if n_atom.res_name == "PRO":
|
|
||||||
# PRO doesn't have an H-atom; do nothing
|
|
||||||
pass
|
|
||||||
else:
|
|
||||||
h_atom = protonate_direction(n_atom, o_atom, c_atom)
|
|
||||||
h_atom.name = "H"
|
|
||||||
return [new_o_atom, new_c_atom]
|
|
||||||
|
|
||||||
|
|
||||||
def protonate_direction(x1_atom, x2_atom, x3_atom):
|
|
||||||
"""Protonates an atom, x1_atom, given a direction.
|
|
||||||
|
|
||||||
New direction for x1_atom proton is (x2_atom -> x3_atom).
|
|
||||||
|
|
||||||
Args:
|
|
||||||
x1_atom: atom to be protonated
|
|
||||||
x2_atom: atom for direction
|
|
||||||
x3_atom: other atom for direction
|
|
||||||
Returns:
|
|
||||||
new hydrogen atom
|
|
||||||
"""
|
|
||||||
dx = (x3_atom.x - x2_atom.x)
|
|
||||||
dy = (x3_atom.y - x2_atom.y)
|
|
||||||
dz = (x3_atom.z - x2_atom.z)
|
|
||||||
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
|
||||||
x = x1_atom.x + dx/length
|
|
||||||
y = x1_atom.y + dy/length
|
|
||||||
z = x1_atom.z + dz/length
|
|
||||||
h_atom = make_new_h(x1_atom, x, y, z)
|
|
||||||
h_atom.name = "H"
|
|
||||||
return h_atom
|
|
||||||
|
|
||||||
|
|
||||||
def protonate_average_direction(x1_atom, x2_atom, x3_atom):
|
|
||||||
"""Protonates an atom, x1_atom, given a direction.
|
|
||||||
|
|
||||||
New direction for x1_atom is (x1_atom/x2_atom -> x3_atom).
|
|
||||||
Note, this one uses the average of x1_atom & x2_atom (N & O) unlike
|
|
||||||
the previous N - C = O
|
|
||||||
|
|
||||||
Args:
|
|
||||||
x1_atom: atom to be protonated
|
|
||||||
x2_atom: atom for direction
|
|
||||||
x3_atom: other atom for direction
|
|
||||||
Returns:
|
|
||||||
new hydrogen atom
|
|
||||||
"""
|
|
||||||
dx = (x3_atom.x + x1_atom.x)*0.5 - x2_atom.x
|
|
||||||
dy = (x3_atom.y + x1_atom.y)*0.5 - x2_atom.y
|
|
||||||
dz = (x3_atom.z + x1_atom.z)*0.5 - x2_atom.z
|
|
||||||
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
|
||||||
x = x1_atom.x + dx/length
|
|
||||||
y = x1_atom.y + dy/length
|
|
||||||
z = x1_atom.z + dz/length
|
|
||||||
h_atom = make_new_h(x1_atom, x, y, z)
|
|
||||||
h_atom.name = "H"
|
|
||||||
return h_atom
|
|
||||||
|
|
||||||
|
|
||||||
def protonate_sp2(x1_atom, x2_atom, x3_atom):
|
|
||||||
"""Protonates a SP2 atom, given a list of atoms
|
|
||||||
|
|
||||||
Args:
|
|
||||||
x1_atom: atom to set direction
|
|
||||||
x2_atom: atom to be protonated
|
|
||||||
x3_atom: other atom to set direction
|
|
||||||
Returns:
|
|
||||||
new hydrogen atom
|
|
||||||
"""
|
|
||||||
dx = (x1_atom.x + x3_atom.x)*0.5 - x2_atom.x
|
|
||||||
dy = (x1_atom.y + x3_atom.y)*0.5 - x2_atom.y
|
|
||||||
dz = (x1_atom.z + x3_atom.z)*0.5 - x2_atom.z
|
|
||||||
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
|
||||||
x = x2_atom.x - dx/length
|
|
||||||
y = x2_atom.y - dy/length
|
|
||||||
z = x2_atom.z - dz/length
|
|
||||||
h_atom = make_new_h(x2_atom, x, y, z)
|
|
||||||
h_atom.name = "H"
|
|
||||||
return h_atom
|
|
||||||
|
|
||||||
|
|
||||||
def make_new_h(atom, x, y, z):
|
|
||||||
"""Add a new hydrogen to an atom at the specified position.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
atom: atom to protonate
|
|
||||||
x: x position of hydrogen
|
|
||||||
y: y position of hydrogen
|
|
||||||
z: z position of hydrogen
|
|
||||||
Returns:
|
|
||||||
new hydrogen atom
|
|
||||||
"""
|
|
||||||
new_h = propka.atom.Atom()
|
|
||||||
new_h.set_property(
|
|
||||||
numb=None, name='H{0:s}'.format(atom.name[1:]),
|
|
||||||
res_name=atom.res_name, chain_id=atom.chain_id,
|
|
||||||
res_num=atom.res_num, x=x, y=y, z=z, occ=None, beta=None)
|
|
||||||
new_h.element = 'H'
|
|
||||||
new_h.bonded_atoms = [atom]
|
|
||||||
new_h.charge = 0
|
|
||||||
new_h.steric_number = 0
|
|
||||||
new_h.number_of_lone_pairs = 0
|
|
||||||
new_h.number_of_protons_to_add = 0
|
|
||||||
new_h.num_pi_elec_2_3_bonds = 0
|
|
||||||
atom.bonded_atoms.append(new_h)
|
|
||||||
atom.conformation_container.add_atom(new_h)
|
|
||||||
return new_h
|
|
||||||
|
|
||||||
|
|
||||||
# TODO - the remaining functions form a distinct "module" for desolvation
|
|
||||||
|
|
||||||
|
|
||||||
# TODO - I have no idea what these constants mean so I labeled them "UNK_"
|
|
||||||
UNK_MIN_DISTANCE = 2.75
|
|
||||||
MIN_DISTANCE_4TH = math.pow(UNK_MIN_DISTANCE, 4)
|
|
||||||
UNK_DIELECTRIC1 = 160
|
|
||||||
UNK_DIELECTRIC2 = 30
|
|
||||||
UNK_PKA_SCALING1 = 244.12
|
|
||||||
UNK_BACKBONE_DISTANCE1 = 6.0
|
|
||||||
UNK_BACKBONE_DISTANCE2 = 3.0
|
|
||||||
UNK_FANGLE_MIN = 0.001
|
|
||||||
UNK_PKA_SCALING2 = 0.8
|
|
||||||
COMBINED_NUM_BURIED_MAX = 900
|
|
||||||
SEPARATE_NUM_BURIED_MAX = 400
|
|
||||||
|
|
||||||
|
|
||||||
def radial_volume_desolvation(parameters, group):
|
|
||||||
"""Calculate desolvation terms for group.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: parameters for desolvation calculation
|
|
||||||
group: group of atoms for calculation
|
|
||||||
"""
|
|
||||||
all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms()
|
|
||||||
volume = 0.0
|
|
||||||
# TODO - Nathan really wants to rename the num_volume variable.
|
|
||||||
# He had to re-read the original paper to figure out what it was.
|
|
||||||
# A better name would be num_volume.
|
|
||||||
group.num_volume = 0
|
|
||||||
min_dist_4th = MIN_DISTANCE_4TH
|
|
||||||
for atom in all_atoms:
|
|
||||||
# ignore atoms in the same residue
|
|
||||||
if (atom.res_num == group.atom.res_num
|
|
||||||
and atom.chain_id == group.atom.chain_id):
|
|
||||||
continue
|
|
||||||
sq_dist = squared_distance(group, atom)
|
|
||||||
# desolvation
|
|
||||||
if sq_dist < parameters.desolv_cutoff_squared:
|
|
||||||
# use a default relative volume of 1.0 if the volume of the element
|
|
||||||
# is not found in parameters
|
|
||||||
# insert check for methyl groups
|
|
||||||
if atom.element == 'C' and atom.name not in ['CA', 'C']:
|
|
||||||
dvol = parameters.VanDerWaalsVolume['C4']
|
|
||||||
else:
|
|
||||||
dvol = parameters.VanDerWaalsVolume.get(atom.element, 1.0)
|
|
||||||
dv_inc = dvol/max(min_dist_4th, sq_dist*sq_dist)
|
|
||||||
volume += dv_inc
|
|
||||||
# buried
|
|
||||||
if sq_dist < parameters.buried_cutoff_squared:
|
|
||||||
group.num_volume += 1
|
|
||||||
group.buried = calculate_weight(parameters, group.num_volume)
|
|
||||||
scale_factor = calculate_scale_factor(parameters, group.buried)
|
|
||||||
volume_after_allowance = max(0.00, volume-parameters.desolvationAllowance)
|
|
||||||
group.energy_volume = (
|
|
||||||
group.charge * parameters.desolvationPrefactor
|
|
||||||
* volume_after_allowance * scale_factor)
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_scale_factor(parameters, weight):
|
|
||||||
"""Calculate desolvation scaling factor.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: parameters for desolvation calculation
|
|
||||||
weight: weight for scaling factor
|
|
||||||
Returns:
|
|
||||||
scaling factor
|
|
||||||
"""
|
|
||||||
scale_factor = 1.0 - (1.0 - parameters.desolvationSurfaceScalingFactor)*(1.0 - weight)
|
|
||||||
return scale_factor
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_weight(parameters, num_volume):
|
|
||||||
"""Calculate the atom-based desolvation weight.
|
|
||||||
|
|
||||||
TODO - figure out why a similar function exists in version.py
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: parameters for desolvation calculation
|
|
||||||
num_volume: number of heavy atoms within desolvation calculation volume
|
|
||||||
Returns:
|
|
||||||
desolvation weight
|
|
||||||
"""
|
|
||||||
weight = (
|
|
||||||
float(num_volume - parameters.Nmin)
|
|
||||||
/ float(parameters.Nmax - parameters.Nmin))
|
|
||||||
weight = min(1.0, weight)
|
|
||||||
weight = max(0.0, weight)
|
|
||||||
return weight
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_pair_weight(parameters, num_volume1, num_volume2):
|
|
||||||
"""Calculate the atom-pair based desolvation weight.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
num_volume1: number of heavy atoms within first desolvation volume
|
|
||||||
num_volume2: number of heavy atoms within second desolvation volume
|
|
||||||
Returns:
|
|
||||||
desolvation weight
|
|
||||||
"""
|
|
||||||
num_volume = num_volume1 + num_volume2
|
|
||||||
num_min = 2*parameters.Nmin
|
|
||||||
num_max = 2*parameters.Nmax
|
|
||||||
weight = float(num_volume - num_min)/float(num_max - num_min)
|
|
||||||
weight = min(1.0, weight)
|
|
||||||
weight = max(0.0, weight)
|
|
||||||
return weight
|
|
||||||
|
|
||||||
|
|
||||||
def hydrogen_bond_energy(dist, dpka_max, cutoffs, f_angle=1.0):
|
|
||||||
"""Calculate hydrogen-bond interaction pKa shift.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
dist: distance for hydrogen bond
|
|
||||||
dpka_max: maximum pKa value shift
|
|
||||||
cutoffs: array with max and min distance values
|
|
||||||
f_angle: angle scaling factor
|
|
||||||
Returns:
|
|
||||||
pKa shift value
|
|
||||||
"""
|
|
||||||
if dist < cutoffs[0]:
|
|
||||||
value = 1.00
|
|
||||||
elif dist > cutoffs[1]:
|
|
||||||
value = 0.00
|
|
||||||
else:
|
|
||||||
value = 1.0 - (dist - cutoffs[0])/(cutoffs[1] - cutoffs[0])
|
|
||||||
dpka = dpka_max*value*f_angle
|
|
||||||
return abs(dpka)
|
|
||||||
|
|
||||||
|
|
||||||
def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None):
|
|
||||||
"""Calculate distance and angle factors for three atoms for backbone interactions.
|
|
||||||
|
|
||||||
NOTE - you need to use atom1 to be the e.g. ASP atom if distance is reset at
|
|
||||||
return: [O1 -- H2-N3].
|
|
||||||
|
|
||||||
Also generalized to be able to be used for residue 'centers' for C=O COO interactions.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
atom1: first atom for calculation (could be None)
|
|
||||||
atom2: second atom for calculation
|
|
||||||
atom3: third atom for calculation
|
|
||||||
center: center point between atoms 1 and 2
|
|
||||||
Returns
|
|
||||||
[distance factor between atoms 1 and 2,
|
|
||||||
angle factor,
|
|
||||||
distance factor between atoms 2 and 3]
|
|
||||||
"""
|
|
||||||
# The angle factor
|
|
||||||
#
|
|
||||||
# ---closest_atom1/2
|
|
||||||
# .
|
|
||||||
# .
|
|
||||||
# the_hydrogen---closest_atom2/1---
|
|
||||||
dx_32 = atom2.x - atom3.x
|
|
||||||
dy_32 = atom2.y - atom3.y
|
|
||||||
dz_32 = atom2.z - atom3.z
|
|
||||||
dist_23 = math.sqrt(dx_32 * dx_32 + dy_32 * dy_32 + dz_32 * dz_32)
|
|
||||||
dx_32 = dx_32/dist_23
|
|
||||||
dy_32 = dy_32/dist_23
|
|
||||||
dz_32 = dz_32/dist_23
|
|
||||||
if atom1 is None:
|
|
||||||
dx_21 = center[0] - atom2.x
|
|
||||||
dy_21 = center[1] - atom2.y
|
|
||||||
dz_21 = center[2] - atom2.z
|
|
||||||
else:
|
|
||||||
dx_21 = atom1.x - atom2.x
|
|
||||||
dy_21 = atom1.y - atom2.y
|
|
||||||
dz_21 = atom1.z - atom2.z
|
|
||||||
dist_12 = math.sqrt(dx_21 * dx_21 + dy_21 * dy_21 + dz_21 * dz_21)
|
|
||||||
dx_21 = dx_21/dist_12
|
|
||||||
dy_21 = dy_21/dist_12
|
|
||||||
dz_21 = dz_21/dist_12
|
|
||||||
f_angle = dx_21*dx_32 + dy_21*dy_32 + dz_21*dz_32
|
|
||||||
return dist_12, f_angle, dist_23
|
|
||||||
|
|
||||||
|
|
||||||
def hydrogen_bond_interaction(group1, group2, version):
|
|
||||||
"""Calculate energy for hydrogen bond interactions between two groups.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first interacting group
|
|
||||||
group2: second interacting group
|
|
||||||
version: an object that contains version-specific parameters
|
|
||||||
Returns:
|
|
||||||
hydrogen bond interaction energy
|
|
||||||
"""
|
|
||||||
# find the smallest distance between interacting atoms
|
|
||||||
atoms1 = group1.get_interaction_atoms(group2)
|
|
||||||
atoms2 = group2.get_interaction_atoms(group1)
|
|
||||||
[closest_atom1, dist, closest_atom2] = get_smallest_distance(atoms1, atoms2)
|
|
||||||
if None in [closest_atom1, closest_atom2]:
|
|
||||||
warning(
|
|
||||||
'Side chain interaction failed for {0:s} and {1:s}'.format(
|
|
||||||
group1.label, group2.label))
|
|
||||||
return None
|
|
||||||
# get the parameters
|
|
||||||
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,
|
|
||||||
closest_atom2)
|
|
||||||
if (dpka_max is None) or (None in cutoff):
|
|
||||||
return None
|
|
||||||
# check that the closest atoms are close enough
|
|
||||||
if dist >= cutoff[1]:
|
|
||||||
return None
|
|
||||||
# check that bond distance criteria is met
|
|
||||||
min_hbond_dist = version.parameters.min_bond_distance_for_hydrogen_bonds
|
|
||||||
if group1.atom.is_atom_within_bond_distance(group2.atom, min_hbond_dist, 1):
|
|
||||||
return None
|
|
||||||
# set angle factor
|
|
||||||
f_angle = 1.0
|
|
||||||
if group2.type in version.parameters.angular_dependent_sidechain_interactions:
|
|
||||||
if closest_atom2.element == 'H':
|
|
||||||
heavy_atom = closest_atom2.bonded_atoms[0]
|
|
||||||
hydrogen = closest_atom2
|
|
||||||
dist, f_angle, _ = angle_distance_factors(closest_atom1, hydrogen,
|
|
||||||
heavy_atom)
|
|
||||||
else:
|
|
||||||
# Either the structure is corrupt (no hydrogen), or the heavy atom
|
|
||||||
# is closer to the titratable atom than the hydrogen. In either
|
|
||||||
# case, we set the angle factor to 0
|
|
||||||
f_angle = 0.0
|
|
||||||
elif group1.type in version.parameters.angular_dependent_sidechain_interactions:
|
|
||||||
if closest_atom1.element == 'H':
|
|
||||||
heavy_atom = closest_atom1.bonded_atoms[0]
|
|
||||||
hydrogen = closest_atom1
|
|
||||||
dist, f_angle, _ = angle_distance_factors(closest_atom2, hydrogen,
|
|
||||||
heavy_atom)
|
|
||||||
else:
|
|
||||||
# Either the structure is corrupt (no hydrogen), or the heavy atom
|
|
||||||
# is closer to the titratable atom than the hydrogen. In either
|
|
||||||
# case, we set the angle factor to 0
|
|
||||||
f_angle = 0.0
|
|
||||||
weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume)
|
|
||||||
exception, value = version.check_exceptions(group1, group2)
|
|
||||||
if exception:
|
|
||||||
# Do nothing, value should have been assigned.
|
|
||||||
pass
|
|
||||||
else:
|
|
||||||
value = version.calculate_side_chain_energy(
|
|
||||||
dist, dpka_max, cutoff, weight, f_angle)
|
|
||||||
return value
|
|
||||||
|
|
||||||
|
|
||||||
def electrostatic_interaction(group1, group2, dist, version):
|
|
||||||
"""Calculate electrostatic interaction betwee two groups.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first interacting group
|
|
||||||
group2: second interacting group
|
|
||||||
dist: distance between groups
|
|
||||||
version: version-specific object with parameters and functions
|
|
||||||
Returns:
|
|
||||||
electrostatic interaction energy or None (if no interaction is appropriate)
|
|
||||||
"""
|
|
||||||
# check if we should do coulomb interaction at all
|
|
||||||
if not version.check_coulomb_pair(group1, group2, dist):
|
|
||||||
return None
|
|
||||||
weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume)
|
|
||||||
value = version.calculate_coulomb_energy(dist, weight)
|
|
||||||
return value
|
|
||||||
|
|
||||||
|
|
||||||
def check_coulomb_pair(parameters, group1, group2, dist):
|
|
||||||
"""Checks if this Coulomb interaction should be done.
|
|
||||||
|
|
||||||
NOTE - this is a propka2.0 hack
|
|
||||||
TODO - figure out why a similar function exists in version.py
|
|
||||||
|
|
||||||
Args:
|
|
||||||
parameters: parameters for Coulomb calculations
|
|
||||||
group1: first interacting group
|
|
||||||
group2: second interacting group
|
|
||||||
dist: distance between groups
|
|
||||||
Returns:
|
|
||||||
Boolean
|
|
||||||
"""
|
|
||||||
num_volume = group1.num_volume + group2.num_volume
|
|
||||||
do_coulomb = True
|
|
||||||
# check if both groups are titratable (ions are taken care of in
|
|
||||||
# determinants::set_ion_determinants)
|
|
||||||
if not (group1.titratable and group2.titratable):
|
|
||||||
do_coulomb = False
|
|
||||||
# check if the distance is not too big
|
|
||||||
if dist > parameters.coulomb_cutoff2:
|
|
||||||
do_coulomb = False
|
|
||||||
# check that desolvation is ok
|
|
||||||
if num_volume < parameters.Nmin:
|
|
||||||
do_coulomb = False
|
|
||||||
return do_coulomb
|
|
||||||
|
|
||||||
|
|
||||||
def coulomb_energy(dist, weight, parameters):
|
|
||||||
"""Calculates the Coulomb interaction pKa shift based on Coulomb's law.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
dist: distance for electrostatic interaction
|
|
||||||
weight: scaling of dielectric constant
|
|
||||||
parameters: parameter object for calculation
|
|
||||||
Returns:
|
|
||||||
pKa shift
|
|
||||||
"""
|
|
||||||
diel = UNK_DIELECTRIC1 - (UNK_DIELECTRIC1 - UNK_DIELECTRIC2)*weight
|
|
||||||
dist = max(dist, parameters.coulomb_cutoff1)
|
|
||||||
scale = (
|
|
||||||
(dist - parameters.coulomb_cutoff2)
|
|
||||||
/ (parameters.coulomb_cutoff1 - parameters.coulomb_cutoff2))
|
|
||||||
scale = max(0.0, scale)
|
|
||||||
scale = min(1.0, scale)
|
|
||||||
dpka = UNK_PKA_SCALING1/(diel*dist)*scale
|
|
||||||
return abs(dpka)
|
|
||||||
|
|
||||||
|
|
||||||
def backbone_reorganization(_, conformation):
|
|
||||||
"""Perform calculations related to backbone reorganizations.
|
|
||||||
|
|
||||||
NOTE - this was described in the code as "adding test stuff"
|
|
||||||
NOTE - this function does not appear to be used
|
|
||||||
TODO - figure out why a similar function exists in version.py
|
|
||||||
|
|
||||||
Args:
|
|
||||||
_: not used
|
|
||||||
conformation: specific molecule conformation
|
|
||||||
"""
|
|
||||||
titratable_groups = conformation.get_backbone_reorganisation_groups()
|
|
||||||
bbc_groups = conformation.get_backbone_co_groups()
|
|
||||||
|
|
||||||
for titratable_group in titratable_groups:
|
|
||||||
weight = titratable_group.buried
|
|
||||||
dpka = 0.00
|
|
||||||
for bbc_group in bbc_groups:
|
|
||||||
center = [titratable_group.x, titratable_group.y, titratable_group.z]
|
|
||||||
atom2 = bbc_group.get_interaction_atoms(titratable_group)[0]
|
|
||||||
dist, f_angle, _ = angle_distance_factors(atom2=atom2,
|
|
||||||
atom3=bbc_group.atom,
|
|
||||||
center=center)
|
|
||||||
if dist < UNK_BACKBONE_DISTANCE1 and f_angle > UNK_FANGLE_MIN:
|
|
||||||
value = (
|
|
||||||
1.0 - (dist-UNK_BACKBONE_DISTANCE2)
|
|
||||||
/ (UNK_BACKBONE_DISTANCE1-UNK_BACKBONE_DISTANCE2))
|
|
||||||
dpka += UNK_PKA_SCALING2 * min(1.0, value)
|
|
||||||
titratable_group.energy_local = dpka*weight
|
|
||||||
|
|
||||||
|
|
||||||
def check_exceptions(version, group1, group2):
|
|
||||||
"""Checks for atypical behavior in interactions between two groups.
|
|
||||||
Checks are made based on group type.
|
|
||||||
|
|
||||||
TODO - figure out why a similar function exists in version.py
|
|
||||||
|
|
||||||
Args:
|
|
||||||
version: version object
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
res_type1 = group1.type
|
|
||||||
res_type2 = group2.type
|
|
||||||
if (res_type1 == "COO") and (res_type2 == "ARG"):
|
|
||||||
exception, value = check_coo_arg_exception(group1, group2, version)
|
|
||||||
elif (res_type1 == "ARG") and (res_type2 == "COO"):
|
|
||||||
exception, value = check_coo_arg_exception(group2, group1, version)
|
|
||||||
elif (res_type1 == "COO") and (res_type2 == "COO"):
|
|
||||||
exception, value = check_coo_coo_exception(group1, group2, version)
|
|
||||||
elif (res_type1 == "CYS") and (res_type2 == "CYS"):
|
|
||||||
exception, value = check_cys_cys_exception(group1, group2, version)
|
|
||||||
elif ((res_type1 == "COO") and (res_type2 == "HIS")
|
|
||||||
or (res_type1 == "HIS") and (res_type2 == "COO")):
|
|
||||||
exception, value = check_coo_his_exception(group1, group2, version)
|
|
||||||
elif ((res_type1 == "OCO") and (res_type2 == "HIS")
|
|
||||||
or (res_type1 == "HIS") and (res_type2 == "OCO")):
|
|
||||||
exception, value = check_oco_his_exception(group1, group2, version)
|
|
||||||
elif ((res_type1 == "CYS") and (res_type2 == "HIS")
|
|
||||||
or (res_type1 == "HIS") and (res_type2 == "CYS")):
|
|
||||||
exception, value = check_cys_his_exception(group1, group2, version)
|
|
||||||
else:
|
|
||||||
# do nothing, no exception for this pair
|
|
||||||
exception = False
|
|
||||||
value = None
|
|
||||||
return exception, value
|
|
||||||
|
|
||||||
|
|
||||||
def check_coo_arg_exception(group_coo, group_arg, version):
|
|
||||||
"""Check for COO-ARG interaction atypical behavior.
|
|
||||||
|
|
||||||
Uses the two shortest unique distances (involving 2+2 atoms)
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group_coo: COO group
|
|
||||||
group_arg: ARG group
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = True
|
|
||||||
value_tot = 0.00
|
|
||||||
# needs to be this way since you want to find shortest distance first
|
|
||||||
atoms_coo = []
|
|
||||||
atoms_coo.extend(group_coo.get_interaction_atoms(group_arg))
|
|
||||||
atoms_arg = []
|
|
||||||
atoms_arg.extend(group_arg.get_interaction_atoms(group_coo))
|
|
||||||
for _ in ["shortest", "runner-up"]:
|
|
||||||
# find the closest interaction pair
|
|
||||||
[closest_coo_atom, dist, closest_arg_atom] = get_smallest_distance(atoms_coo,
|
|
||||||
atoms_arg)
|
|
||||||
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_coo_atom,
|
|
||||||
closest_arg_atom)
|
|
||||||
# calculate and sum up interaction energy
|
|
||||||
f_angle = 1.00
|
|
||||||
if group_arg.type in version.parameters.angular_dependent_sidechain_interactions:
|
|
||||||
atom3 = closest_arg_atom.bonded_atoms[0]
|
|
||||||
dist, f_angle, _ = angle_distance_factors(closest_coo_atom,
|
|
||||||
closest_arg_atom,
|
|
||||||
atom3)
|
|
||||||
value = hydrogen_bond_energy(dist, dpka_max, cutoff, f_angle)
|
|
||||||
value_tot += value
|
|
||||||
# remove closest atoms before we attemp to find the runner-up pair
|
|
||||||
atoms_coo.remove(closest_coo_atom)
|
|
||||||
atoms_arg.remove(closest_arg_atom)
|
|
||||||
return exception, value_tot
|
|
||||||
|
|
||||||
|
|
||||||
def check_coo_coo_exception(group1, group2, version):
|
|
||||||
"""Check for COO-COO hydrogen-bond atypical interaction behavior.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = True
|
|
||||||
interact_groups12 = group1.get_interaction_atoms(group2)
|
|
||||||
interact_groups21 = group2.get_interaction_atoms(group1)
|
|
||||||
[closest_atom1, dist, closest_atom2] = get_smallest_distance(interact_groups12,
|
|
||||||
interact_groups21)
|
|
||||||
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,
|
|
||||||
closest_atom2)
|
|
||||||
f_angle = 1.00
|
|
||||||
value = hydrogen_bond_energy(dist, dpka_max, cutoff, f_angle)
|
|
||||||
weight = calculate_pair_weight(version.parameters, group1.num_volume, group2.num_volume)
|
|
||||||
value = value * (1.0 + weight)
|
|
||||||
return exception, value
|
|
||||||
|
|
||||||
|
|
||||||
def check_coo_his_exception(group1, group2, version):
|
|
||||||
"""Check for COO-HIS atypical interaction behavior
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = False
|
|
||||||
if check_buried(group1.num_volume, group2.num_volume):
|
|
||||||
exception = True
|
|
||||||
return exception, version.parameters.COO_HIS_exception
|
|
||||||
|
|
||||||
|
|
||||||
def check_oco_his_exception(group1, group2, version):
|
|
||||||
"""Check for OCO-HIS atypical interaction behavior
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = False
|
|
||||||
if check_buried(group1.num_volume, group2.num_volume):
|
|
||||||
exception = True
|
|
||||||
return exception, version.parameters.OCO_HIS_exception
|
|
||||||
|
|
||||||
|
|
||||||
def check_cys_his_exception(group1, group2, version):
|
|
||||||
"""Check for CYS-HIS atypical interaction behavior
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = False
|
|
||||||
if check_buried(group1.num_volume, group2.num_volume):
|
|
||||||
exception = True
|
|
||||||
return exception, version.parameters.CYS_HIS_exception
|
|
||||||
|
|
||||||
|
|
||||||
def check_cys_cys_exception(group1, group2, version):
|
|
||||||
"""Check for CYS-CYS atypical interaction behavior
|
|
||||||
|
|
||||||
Args:
|
|
||||||
group1: first group for check
|
|
||||||
group2: second group for check
|
|
||||||
version: version object
|
|
||||||
Returns:
|
|
||||||
1. Boolean indicating atypical behavior,
|
|
||||||
2. value associated with atypical interaction (None if Boolean is False)
|
|
||||||
"""
|
|
||||||
exception = False
|
|
||||||
if check_buried(group1.num_volume, group2.num_volume):
|
|
||||||
exception = True
|
|
||||||
return exception, version.parameters.CYS_CYS_exception
|
|
||||||
|
|
||||||
|
|
||||||
def check_buried(num_volume1, num_volume2):
|
|
||||||
"""Check to see if an interaction is buried
|
|
||||||
|
|
||||||
Args:
|
|
||||||
num_volume1: number of buried heavy atoms in volume 1
|
|
||||||
num_volume2: number of buried heavy atoms in volume 2
|
|
||||||
Returns:
|
|
||||||
True if interaction is buried, False otherwise
|
|
||||||
"""
|
|
||||||
if ((num_volume1 + num_volume2 <= COMBINED_NUM_BURIED_MAX)
|
|
||||||
and (num_volume1 <= SEPARATE_NUM_BURIED_MAX
|
|
||||||
or num_volume2 <= SEPARATE_NUM_BURIED_MAX)):
|
|
||||||
return False
|
|
||||||
return True
|
|
||||||
|
|||||||
@@ -8,7 +8,7 @@ import propka.iterative
|
|||||||
import propka.lib
|
import propka.lib
|
||||||
import propka.vector_algebra
|
import propka.vector_algebra
|
||||||
from propka.calculations import squared_distance, get_smallest_distance
|
from propka.calculations import squared_distance, get_smallest_distance
|
||||||
from propka.calculations import angle_distance_factors, hydrogen_bond_energy
|
from propka.energy import angle_distance_factors, hydrogen_bond_energy
|
||||||
from propka.determinant import Determinant
|
from propka.determinant import Determinant
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
531
propka/energy.py
Normal file
531
propka/energy.py
Normal file
@@ -0,0 +1,531 @@
|
|||||||
|
"""Energy calculations."""
|
||||||
|
import math
|
||||||
|
from propka.lib import warning
|
||||||
|
from propka.calculations import squared_distance, get_smallest_distance
|
||||||
|
|
||||||
|
|
||||||
|
# TODO - I have no idea what these constants mean so I labeled them "UNK_"
|
||||||
|
UNK_MIN_DISTANCE = 2.75
|
||||||
|
MIN_DISTANCE_4TH = math.pow(UNK_MIN_DISTANCE, 4)
|
||||||
|
UNK_DIELECTRIC1 = 160
|
||||||
|
UNK_DIELECTRIC2 = 30
|
||||||
|
UNK_PKA_SCALING1 = 244.12
|
||||||
|
UNK_BACKBONE_DISTANCE1 = 6.0
|
||||||
|
UNK_BACKBONE_DISTANCE2 = 3.0
|
||||||
|
UNK_FANGLE_MIN = 0.001
|
||||||
|
UNK_PKA_SCALING2 = 0.8
|
||||||
|
COMBINED_NUM_BURIED_MAX = 900
|
||||||
|
SEPARATE_NUM_BURIED_MAX = 400
|
||||||
|
|
||||||
|
|
||||||
|
def radial_volume_desolvation(parameters, group):
|
||||||
|
"""Calculate desolvation terms for group.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
parameters: parameters for desolvation calculation
|
||||||
|
group: group of atoms for calculation
|
||||||
|
"""
|
||||||
|
all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms()
|
||||||
|
volume = 0.0
|
||||||
|
group.num_volume = 0
|
||||||
|
min_dist_4th = MIN_DISTANCE_4TH
|
||||||
|
for atom in all_atoms:
|
||||||
|
# ignore atoms in the same residue
|
||||||
|
if (atom.res_num == group.atom.res_num
|
||||||
|
and atom.chain_id == group.atom.chain_id):
|
||||||
|
continue
|
||||||
|
sq_dist = squared_distance(group, atom)
|
||||||
|
# desolvation
|
||||||
|
if sq_dist < parameters.desolv_cutoff_squared:
|
||||||
|
# use a default relative volume of 1.0 if the volume of the element
|
||||||
|
# is not found in parameters
|
||||||
|
# insert check for methyl groups
|
||||||
|
if atom.element == 'C' and atom.name not in ['CA', 'C']:
|
||||||
|
dvol = parameters.VanDerWaalsVolume['C4']
|
||||||
|
else:
|
||||||
|
dvol = parameters.VanDerWaalsVolume.get(atom.element, 1.0)
|
||||||
|
dv_inc = dvol/max(min_dist_4th, sq_dist*sq_dist)
|
||||||
|
volume += dv_inc
|
||||||
|
# buried
|
||||||
|
if sq_dist < parameters.buried_cutoff_squared:
|
||||||
|
group.num_volume += 1
|
||||||
|
group.buried = calculate_weight(parameters, group.num_volume)
|
||||||
|
scale_factor = calculate_scale_factor(parameters, group.buried)
|
||||||
|
volume_after_allowance = max(0.00, volume-parameters.desolvationAllowance)
|
||||||
|
group.energy_volume = (
|
||||||
|
group.charge * parameters.desolvationPrefactor
|
||||||
|
* volume_after_allowance * scale_factor)
|
||||||
|
|
||||||
|
|
||||||
|
def calculate_scale_factor(parameters, weight):
|
||||||
|
"""Calculate desolvation scaling factor.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
parameters: parameters for desolvation calculation
|
||||||
|
weight: weight for scaling factor
|
||||||
|
Returns:
|
||||||
|
scaling factor
|
||||||
|
"""
|
||||||
|
scale_factor = 1.0 - (1.0 - parameters.desolvationSurfaceScalingFactor)*(1.0 - weight)
|
||||||
|
return scale_factor
|
||||||
|
|
||||||
|
|
||||||
|
def calculate_weight(parameters, num_volume):
|
||||||
|
"""Calculate the atom-based desolvation weight.
|
||||||
|
|
||||||
|
TODO - figure out why a similar function exists in version.py
|
||||||
|
|
||||||
|
Args:
|
||||||
|
parameters: parameters for desolvation calculation
|
||||||
|
num_volume: number of heavy atoms within desolvation calculation volume
|
||||||
|
Returns:
|
||||||
|
desolvation weight
|
||||||
|
"""
|
||||||
|
weight = (
|
||||||
|
float(num_volume - parameters.Nmin)
|
||||||
|
/ float(parameters.Nmax - parameters.Nmin))
|
||||||
|
weight = min(1.0, weight)
|
||||||
|
weight = max(0.0, weight)
|
||||||
|
return weight
|
||||||
|
|
||||||
|
|
||||||
|
def calculate_pair_weight(parameters, num_volume1, num_volume2):
|
||||||
|
"""Calculate the atom-pair based desolvation weight.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
num_volume1: number of heavy atoms within first desolvation volume
|
||||||
|
num_volume2: number of heavy atoms within second desolvation volume
|
||||||
|
Returns:
|
||||||
|
desolvation weight
|
||||||
|
"""
|
||||||
|
num_volume = num_volume1 + num_volume2
|
||||||
|
num_min = 2*parameters.Nmin
|
||||||
|
num_max = 2*parameters.Nmax
|
||||||
|
weight = float(num_volume - num_min)/float(num_max - num_min)
|
||||||
|
weight = min(1.0, weight)
|
||||||
|
weight = max(0.0, weight)
|
||||||
|
return weight
|
||||||
|
|
||||||
|
|
||||||
|
def hydrogen_bond_energy(dist, dpka_max, cutoffs, f_angle=1.0):
|
||||||
|
"""Calculate hydrogen-bond interaction pKa shift.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
dist: distance for hydrogen bond
|
||||||
|
dpka_max: maximum pKa value shift
|
||||||
|
cutoffs: array with max and min distance values
|
||||||
|
f_angle: angle scaling factor
|
||||||
|
Returns:
|
||||||
|
pKa shift value
|
||||||
|
"""
|
||||||
|
if dist < cutoffs[0]:
|
||||||
|
value = 1.00
|
||||||
|
elif dist > cutoffs[1]:
|
||||||
|
value = 0.00
|
||||||
|
else:
|
||||||
|
value = 1.0 - (dist - cutoffs[0])/(cutoffs[1] - cutoffs[0])
|
||||||
|
dpka = dpka_max*value*f_angle
|
||||||
|
return abs(dpka)
|
||||||
|
|
||||||
|
|
||||||
|
def angle_distance_factors(atom1=None, atom2=None, atom3=None, center=None):
|
||||||
|
"""Calculate distance and angle factors for three atoms for backbone interactions.
|
||||||
|
|
||||||
|
NOTE - you need to use atom1 to be the e.g. ASP atom if distance is reset at
|
||||||
|
return: [O1 -- H2-N3].
|
||||||
|
|
||||||
|
Also generalized to be able to be used for residue 'centers' for C=O COO interactions.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
atom1: first atom for calculation (could be None)
|
||||||
|
atom2: second atom for calculation
|
||||||
|
atom3: third atom for calculation
|
||||||
|
center: center point between atoms 1 and 2
|
||||||
|
Returns
|
||||||
|
[distance factor between atoms 1 and 2,
|
||||||
|
angle factor,
|
||||||
|
distance factor between atoms 2 and 3]
|
||||||
|
"""
|
||||||
|
# The angle factor
|
||||||
|
#
|
||||||
|
# ---closest_atom1/2
|
||||||
|
# .
|
||||||
|
# .
|
||||||
|
# the_hydrogen---closest_atom2/1---
|
||||||
|
dx_32 = atom2.x - atom3.x
|
||||||
|
dy_32 = atom2.y - atom3.y
|
||||||
|
dz_32 = atom2.z - atom3.z
|
||||||
|
dist_23 = math.sqrt(dx_32 * dx_32 + dy_32 * dy_32 + dz_32 * dz_32)
|
||||||
|
dx_32 = dx_32/dist_23
|
||||||
|
dy_32 = dy_32/dist_23
|
||||||
|
dz_32 = dz_32/dist_23
|
||||||
|
if atom1 is None:
|
||||||
|
dx_21 = center[0] - atom2.x
|
||||||
|
dy_21 = center[1] - atom2.y
|
||||||
|
dz_21 = center[2] - atom2.z
|
||||||
|
else:
|
||||||
|
dx_21 = atom1.x - atom2.x
|
||||||
|
dy_21 = atom1.y - atom2.y
|
||||||
|
dz_21 = atom1.z - atom2.z
|
||||||
|
dist_12 = math.sqrt(dx_21 * dx_21 + dy_21 * dy_21 + dz_21 * dz_21)
|
||||||
|
dx_21 = dx_21/dist_12
|
||||||
|
dy_21 = dy_21/dist_12
|
||||||
|
dz_21 = dz_21/dist_12
|
||||||
|
f_angle = dx_21*dx_32 + dy_21*dy_32 + dz_21*dz_32
|
||||||
|
return dist_12, f_angle, dist_23
|
||||||
|
|
||||||
|
|
||||||
|
def hydrogen_bond_interaction(group1, group2, version):
|
||||||
|
"""Calculate energy for hydrogen bond interactions between two groups.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first interacting group
|
||||||
|
group2: second interacting group
|
||||||
|
version: an object that contains version-specific parameters
|
||||||
|
Returns:
|
||||||
|
hydrogen bond interaction energy
|
||||||
|
"""
|
||||||
|
# find the smallest distance between interacting atoms
|
||||||
|
atoms1 = group1.get_interaction_atoms(group2)
|
||||||
|
atoms2 = group2.get_interaction_atoms(group1)
|
||||||
|
[closest_atom1, dist, closest_atom2] = get_smallest_distance(atoms1, atoms2)
|
||||||
|
if None in [closest_atom1, closest_atom2]:
|
||||||
|
warning(
|
||||||
|
'Side chain interaction failed for {0:s} and {1:s}'.format(
|
||||||
|
group1.label, group2.label))
|
||||||
|
return None
|
||||||
|
# get the parameters
|
||||||
|
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,
|
||||||
|
closest_atom2)
|
||||||
|
if (dpka_max is None) or (None in cutoff):
|
||||||
|
return None
|
||||||
|
# check that the closest atoms are close enough
|
||||||
|
if dist >= cutoff[1]:
|
||||||
|
return None
|
||||||
|
# check that bond distance criteria is met
|
||||||
|
min_hbond_dist = version.parameters.min_bond_distance_for_hydrogen_bonds
|
||||||
|
if group1.atom.is_atom_within_bond_distance(group2.atom, min_hbond_dist, 1):
|
||||||
|
return None
|
||||||
|
# set angle factor
|
||||||
|
f_angle = 1.0
|
||||||
|
if group2.type in version.parameters.angular_dependent_sidechain_interactions:
|
||||||
|
if closest_atom2.element == 'H':
|
||||||
|
heavy_atom = closest_atom2.bonded_atoms[0]
|
||||||
|
hydrogen = closest_atom2
|
||||||
|
dist, f_angle, _ = angle_distance_factors(closest_atom1, hydrogen,
|
||||||
|
heavy_atom)
|
||||||
|
else:
|
||||||
|
# Either the structure is corrupt (no hydrogen), or the heavy atom
|
||||||
|
# is closer to the titratable atom than the hydrogen. In either
|
||||||
|
# case, we set the angle factor to 0
|
||||||
|
f_angle = 0.0
|
||||||
|
elif group1.type in version.parameters.angular_dependent_sidechain_interactions:
|
||||||
|
if closest_atom1.element == 'H':
|
||||||
|
heavy_atom = closest_atom1.bonded_atoms[0]
|
||||||
|
hydrogen = closest_atom1
|
||||||
|
dist, f_angle, _ = angle_distance_factors(closest_atom2, hydrogen,
|
||||||
|
heavy_atom)
|
||||||
|
else:
|
||||||
|
# Either the structure is corrupt (no hydrogen), or the heavy atom
|
||||||
|
# is closer to the titratable atom than the hydrogen. In either
|
||||||
|
# case, we set the angle factor to 0
|
||||||
|
f_angle = 0.0
|
||||||
|
weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume)
|
||||||
|
exception, value = version.check_exceptions(group1, group2)
|
||||||
|
if exception:
|
||||||
|
# Do nothing, value should have been assigned.
|
||||||
|
pass
|
||||||
|
else:
|
||||||
|
value = version.calculate_side_chain_energy(
|
||||||
|
dist, dpka_max, cutoff, weight, f_angle)
|
||||||
|
return value
|
||||||
|
|
||||||
|
|
||||||
|
def electrostatic_interaction(group1, group2, dist, version):
|
||||||
|
"""Calculate electrostatic interaction betwee two groups.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first interacting group
|
||||||
|
group2: second interacting group
|
||||||
|
dist: distance between groups
|
||||||
|
version: version-specific object with parameters and functions
|
||||||
|
Returns:
|
||||||
|
electrostatic interaction energy or None (if no interaction is appropriate)
|
||||||
|
"""
|
||||||
|
# check if we should do coulomb interaction at all
|
||||||
|
if not version.check_coulomb_pair(group1, group2, dist):
|
||||||
|
return None
|
||||||
|
weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume)
|
||||||
|
value = version.calculate_coulomb_energy(dist, weight)
|
||||||
|
return value
|
||||||
|
|
||||||
|
|
||||||
|
def check_coulomb_pair(parameters, group1, group2, dist):
|
||||||
|
"""Checks if this Coulomb interaction should be done.
|
||||||
|
|
||||||
|
NOTE - this is a propka2.0 hack
|
||||||
|
TODO - figure out why a similar function exists in version.py
|
||||||
|
|
||||||
|
Args:
|
||||||
|
parameters: parameters for Coulomb calculations
|
||||||
|
group1: first interacting group
|
||||||
|
group2: second interacting group
|
||||||
|
dist: distance between groups
|
||||||
|
Returns:
|
||||||
|
Boolean
|
||||||
|
"""
|
||||||
|
num_volume = group1.num_volume + group2.num_volume
|
||||||
|
do_coulomb = True
|
||||||
|
# check if both groups are titratable (ions are taken care of in
|
||||||
|
# determinants::set_ion_determinants)
|
||||||
|
if not (group1.titratable and group2.titratable):
|
||||||
|
do_coulomb = False
|
||||||
|
# check if the distance is not too big
|
||||||
|
if dist > parameters.coulomb_cutoff2:
|
||||||
|
do_coulomb = False
|
||||||
|
# check that desolvation is ok
|
||||||
|
if num_volume < parameters.Nmin:
|
||||||
|
do_coulomb = False
|
||||||
|
return do_coulomb
|
||||||
|
|
||||||
|
|
||||||
|
def coulomb_energy(dist, weight, parameters):
|
||||||
|
"""Calculates the Coulomb interaction pKa shift based on Coulomb's law.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
dist: distance for electrostatic interaction
|
||||||
|
weight: scaling of dielectric constant
|
||||||
|
parameters: parameter object for calculation
|
||||||
|
Returns:
|
||||||
|
pKa shift
|
||||||
|
"""
|
||||||
|
diel = UNK_DIELECTRIC1 - (UNK_DIELECTRIC1 - UNK_DIELECTRIC2)*weight
|
||||||
|
dist = max(dist, parameters.coulomb_cutoff1)
|
||||||
|
scale = (
|
||||||
|
(dist - parameters.coulomb_cutoff2)
|
||||||
|
/ (parameters.coulomb_cutoff1 - parameters.coulomb_cutoff2))
|
||||||
|
scale = max(0.0, scale)
|
||||||
|
scale = min(1.0, scale)
|
||||||
|
dpka = UNK_PKA_SCALING1/(diel*dist)*scale
|
||||||
|
return abs(dpka)
|
||||||
|
|
||||||
|
|
||||||
|
def backbone_reorganization(_, conformation):
|
||||||
|
"""Perform calculations related to backbone reorganizations.
|
||||||
|
|
||||||
|
NOTE - this was described in the code as "adding test stuff"
|
||||||
|
NOTE - this function does not appear to be used
|
||||||
|
TODO - figure out why a similar function exists in version.py
|
||||||
|
|
||||||
|
Args:
|
||||||
|
_: not used
|
||||||
|
conformation: specific molecule conformation
|
||||||
|
"""
|
||||||
|
titratable_groups = conformation.get_backbone_reorganisation_groups()
|
||||||
|
bbc_groups = conformation.get_backbone_co_groups()
|
||||||
|
|
||||||
|
for titratable_group in titratable_groups:
|
||||||
|
weight = titratable_group.buried
|
||||||
|
dpka = 0.00
|
||||||
|
for bbc_group in bbc_groups:
|
||||||
|
center = [titratable_group.x, titratable_group.y, titratable_group.z]
|
||||||
|
atom2 = bbc_group.get_interaction_atoms(titratable_group)[0]
|
||||||
|
dist, f_angle, _ = angle_distance_factors(atom2=atom2,
|
||||||
|
atom3=bbc_group.atom,
|
||||||
|
center=center)
|
||||||
|
if dist < UNK_BACKBONE_DISTANCE1 and f_angle > UNK_FANGLE_MIN:
|
||||||
|
value = (
|
||||||
|
1.0 - (dist-UNK_BACKBONE_DISTANCE2)
|
||||||
|
/ (UNK_BACKBONE_DISTANCE1-UNK_BACKBONE_DISTANCE2))
|
||||||
|
dpka += UNK_PKA_SCALING2 * min(1.0, value)
|
||||||
|
titratable_group.energy_local = dpka*weight
|
||||||
|
|
||||||
|
|
||||||
|
def check_exceptions(version, group1, group2):
|
||||||
|
"""Checks for atypical behavior in interactions between two groups.
|
||||||
|
Checks are made based on group type.
|
||||||
|
|
||||||
|
TODO - figure out why a similar function exists in version.py
|
||||||
|
|
||||||
|
Args:
|
||||||
|
version: version object
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
res_type1 = group1.type
|
||||||
|
res_type2 = group2.type
|
||||||
|
if (res_type1 == "COO") and (res_type2 == "ARG"):
|
||||||
|
exception, value = check_coo_arg_exception(group1, group2, version)
|
||||||
|
elif (res_type1 == "ARG") and (res_type2 == "COO"):
|
||||||
|
exception, value = check_coo_arg_exception(group2, group1, version)
|
||||||
|
elif (res_type1 == "COO") and (res_type2 == "COO"):
|
||||||
|
exception, value = check_coo_coo_exception(group1, group2, version)
|
||||||
|
elif (res_type1 == "CYS") and (res_type2 == "CYS"):
|
||||||
|
exception, value = check_cys_cys_exception(group1, group2, version)
|
||||||
|
elif ((res_type1 == "COO") and (res_type2 == "HIS")
|
||||||
|
or (res_type1 == "HIS") and (res_type2 == "COO")):
|
||||||
|
exception, value = check_coo_his_exception(group1, group2, version)
|
||||||
|
elif ((res_type1 == "OCO") and (res_type2 == "HIS")
|
||||||
|
or (res_type1 == "HIS") and (res_type2 == "OCO")):
|
||||||
|
exception, value = check_oco_his_exception(group1, group2, version)
|
||||||
|
elif ((res_type1 == "CYS") and (res_type2 == "HIS")
|
||||||
|
or (res_type1 == "HIS") and (res_type2 == "CYS")):
|
||||||
|
exception, value = check_cys_his_exception(group1, group2, version)
|
||||||
|
else:
|
||||||
|
# do nothing, no exception for this pair
|
||||||
|
exception = False
|
||||||
|
value = None
|
||||||
|
return exception, value
|
||||||
|
|
||||||
|
|
||||||
|
def check_coo_arg_exception(group_coo, group_arg, version):
|
||||||
|
"""Check for COO-ARG interaction atypical behavior.
|
||||||
|
|
||||||
|
Uses the two shortest unique distances (involving 2+2 atoms)
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group_coo: COO group
|
||||||
|
group_arg: ARG group
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = True
|
||||||
|
value_tot = 0.00
|
||||||
|
# needs to be this way since you want to find shortest distance first
|
||||||
|
atoms_coo = []
|
||||||
|
atoms_coo.extend(group_coo.get_interaction_atoms(group_arg))
|
||||||
|
atoms_arg = []
|
||||||
|
atoms_arg.extend(group_arg.get_interaction_atoms(group_coo))
|
||||||
|
for _ in ["shortest", "runner-up"]:
|
||||||
|
# find the closest interaction pair
|
||||||
|
[closest_coo_atom, dist, closest_arg_atom] = get_smallest_distance(atoms_coo,
|
||||||
|
atoms_arg)
|
||||||
|
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_coo_atom,
|
||||||
|
closest_arg_atom)
|
||||||
|
# calculate and sum up interaction energy
|
||||||
|
f_angle = 1.00
|
||||||
|
if group_arg.type in version.parameters.angular_dependent_sidechain_interactions:
|
||||||
|
atom3 = closest_arg_atom.bonded_atoms[0]
|
||||||
|
dist, f_angle, _ = angle_distance_factors(closest_coo_atom,
|
||||||
|
closest_arg_atom,
|
||||||
|
atom3)
|
||||||
|
value = hydrogen_bond_energy(dist, dpka_max, cutoff, f_angle)
|
||||||
|
value_tot += value
|
||||||
|
# remove closest atoms before we attemp to find the runner-up pair
|
||||||
|
atoms_coo.remove(closest_coo_atom)
|
||||||
|
atoms_arg.remove(closest_arg_atom)
|
||||||
|
return exception, value_tot
|
||||||
|
|
||||||
|
|
||||||
|
def check_coo_coo_exception(group1, group2, version):
|
||||||
|
"""Check for COO-COO hydrogen-bond atypical interaction behavior.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = True
|
||||||
|
interact_groups12 = group1.get_interaction_atoms(group2)
|
||||||
|
interact_groups21 = group2.get_interaction_atoms(group1)
|
||||||
|
[closest_atom1, dist, closest_atom2] = get_smallest_distance(interact_groups12,
|
||||||
|
interact_groups21)
|
||||||
|
[dpka_max, cutoff] = version.get_hydrogen_bond_parameters(closest_atom1,
|
||||||
|
closest_atom2)
|
||||||
|
f_angle = 1.00
|
||||||
|
value = hydrogen_bond_energy(dist, dpka_max, cutoff, f_angle)
|
||||||
|
weight = calculate_pair_weight(version.parameters, group1.num_volume, group2.num_volume)
|
||||||
|
value = value * (1.0 + weight)
|
||||||
|
return exception, value
|
||||||
|
|
||||||
|
|
||||||
|
def check_coo_his_exception(group1, group2, version):
|
||||||
|
"""Check for COO-HIS atypical interaction behavior
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = False
|
||||||
|
if check_buried(group1.num_volume, group2.num_volume):
|
||||||
|
exception = True
|
||||||
|
return exception, version.parameters.COO_HIS_exception
|
||||||
|
|
||||||
|
|
||||||
|
def check_oco_his_exception(group1, group2, version):
|
||||||
|
"""Check for OCO-HIS atypical interaction behavior
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = False
|
||||||
|
if check_buried(group1.num_volume, group2.num_volume):
|
||||||
|
exception = True
|
||||||
|
return exception, version.parameters.OCO_HIS_exception
|
||||||
|
|
||||||
|
|
||||||
|
def check_cys_his_exception(group1, group2, version):
|
||||||
|
"""Check for CYS-HIS atypical interaction behavior
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = False
|
||||||
|
if check_buried(group1.num_volume, group2.num_volume):
|
||||||
|
exception = True
|
||||||
|
return exception, version.parameters.CYS_HIS_exception
|
||||||
|
|
||||||
|
|
||||||
|
def check_cys_cys_exception(group1, group2, version):
|
||||||
|
"""Check for CYS-CYS atypical interaction behavior
|
||||||
|
|
||||||
|
Args:
|
||||||
|
group1: first group for check
|
||||||
|
group2: second group for check
|
||||||
|
version: version object
|
||||||
|
Returns:
|
||||||
|
1. Boolean indicating atypical behavior,
|
||||||
|
2. value associated with atypical interaction (None if Boolean is False)
|
||||||
|
"""
|
||||||
|
exception = False
|
||||||
|
if check_buried(group1.num_volume, group2.num_volume):
|
||||||
|
exception = True
|
||||||
|
return exception, version.parameters.CYS_CYS_exception
|
||||||
|
|
||||||
|
|
||||||
|
def check_buried(num_volume1, num_volume2):
|
||||||
|
"""Check to see if an interaction is buried
|
||||||
|
|
||||||
|
Args:
|
||||||
|
num_volume1: number of buried heavy atoms in volume 1
|
||||||
|
num_volume2: number of buried heavy atoms in volume 2
|
||||||
|
Returns:
|
||||||
|
True if interaction is buried, False otherwise
|
||||||
|
"""
|
||||||
|
if ((num_volume1 + num_volume2 <= COMBINED_NUM_BURIED_MAX)
|
||||||
|
and (num_volume1 <= SEPARATE_NUM_BURIED_MAX
|
||||||
|
or num_volume2 <= SEPARATE_NUM_BURIED_MAX)):
|
||||||
|
return False
|
||||||
|
return True
|
||||||
346
propka/hydrogens.py
Normal file
346
propka/hydrogens.py
Normal file
@@ -0,0 +1,346 @@
|
|||||||
|
"""Calculations related to hydrogen placement."""
|
||||||
|
import math
|
||||||
|
from propka.lib import info
|
||||||
|
from propka.protonate import Protonate
|
||||||
|
from propka.bonds import BondMaker
|
||||||
|
from propka.atom import Atom
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def setup_bonding_and_protonation(molecular_container):
|
||||||
|
"""Set up bonding and protonation for a molecule.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
parameters: not used
|
||||||
|
molecular_container: molecule container.
|
||||||
|
"""
|
||||||
|
# make bonds
|
||||||
|
my_bond_maker = setup_bonding(molecular_container)
|
||||||
|
# set up ligand atom names
|
||||||
|
set_ligand_atom_names(molecular_container)
|
||||||
|
# apply information on pi electrons
|
||||||
|
my_bond_maker.add_pi_electron_information(molecular_container)
|
||||||
|
# Protonate atoms
|
||||||
|
if molecular_container.options.protonate_all:
|
||||||
|
protonator = Protonate(verbose=False)
|
||||||
|
protonator.protonate(molecular_container)
|
||||||
|
|
||||||
|
|
||||||
|
def setup_bonding(molecular_container):
|
||||||
|
"""Set up bonding for a molecular container.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
molecular_container: the molecular container in question
|
||||||
|
Returns:
|
||||||
|
BondMaker object
|
||||||
|
"""
|
||||||
|
my_bond_maker = BondMaker()
|
||||||
|
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
|
||||||
|
return my_bond_maker
|
||||||
|
|
||||||
|
|
||||||
|
def setup_bonding_and_protonation_30_style(molecular_container):
|
||||||
|
"""Set up bonding for a molecular container.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
molecular_container: the molecular container in question
|
||||||
|
Returns:
|
||||||
|
BondMaker object
|
||||||
|
"""
|
||||||
|
# Protonate atoms
|
||||||
|
protonate_30_style(molecular_container)
|
||||||
|
# make bonds
|
||||||
|
bond_maker = BondMaker()
|
||||||
|
bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
|
||||||
|
return bond_maker
|
||||||
|
|
||||||
|
|
||||||
|
def protonate_30_style(molecular_container):
|
||||||
|
"""Protonate the molecule.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
molecular_container: molecule
|
||||||
|
"""
|
||||||
|
for name in molecular_container.conformation_names:
|
||||||
|
info('Now protonating', name)
|
||||||
|
# split atom into residues
|
||||||
|
curres = -1000000
|
||||||
|
residue = []
|
||||||
|
o_atom = None
|
||||||
|
c_atom = None
|
||||||
|
for atom in molecular_container.conformations[name].atoms:
|
||||||
|
if atom.res_num != curres:
|
||||||
|
curres = atom.res_num
|
||||||
|
if len(residue) > 0:
|
||||||
|
#backbone
|
||||||
|
[o_atom, c_atom] = add_backbone_hydrogen(
|
||||||
|
residue, o_atom, c_atom)
|
||||||
|
#arginine
|
||||||
|
if residue[0].res_name == 'ARG':
|
||||||
|
add_arg_hydrogen(residue)
|
||||||
|
#histidine
|
||||||
|
if residue[0].res_name == 'HIS':
|
||||||
|
add_his_hydrogen(residue)
|
||||||
|
#tryptophan
|
||||||
|
if residue[0].res_name == 'TRP':
|
||||||
|
add_trp_hydrogen(residue)
|
||||||
|
#amides
|
||||||
|
if residue[0].res_name in ['GLN', 'ASN']:
|
||||||
|
add_amd_hydrogen(residue)
|
||||||
|
residue = []
|
||||||
|
if atom.type == 'atom':
|
||||||
|
residue.append(atom)
|
||||||
|
|
||||||
|
|
||||||
|
def set_ligand_atom_names(molecular_container):
|
||||||
|
"""Set names for ligands in molecular container.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
molecular_container: molecular container for ligand names
|
||||||
|
"""
|
||||||
|
for name in molecular_container.conformation_names:
|
||||||
|
molecular_container.conformations[name].set_ligand_atom_names()
|
||||||
|
|
||||||
|
|
||||||
|
def add_arg_hydrogen(residue):
|
||||||
|
"""Adds Arg hydrogen atoms to residues according to the 'old way'.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
residue: arginine residue to protonate
|
||||||
|
Returns:
|
||||||
|
list of hydrogen atoms
|
||||||
|
"""
|
||||||
|
#info('Adding arg H',residue)
|
||||||
|
for atom in residue:
|
||||||
|
if atom.name == "CD":
|
||||||
|
cd_atom = atom
|
||||||
|
elif atom.name == "CZ":
|
||||||
|
cz_atom = atom
|
||||||
|
elif atom.name == "NE":
|
||||||
|
ne_atom = atom
|
||||||
|
elif atom.name == "NH1":
|
||||||
|
nh1_atom = atom
|
||||||
|
elif atom.name == "NH2":
|
||||||
|
nh2_atom = atom
|
||||||
|
h1_atom = protonate_sp2(cd_atom, ne_atom, cz_atom)
|
||||||
|
h1_atom.name = "HE"
|
||||||
|
h2_atom = protonate_direction(nh1_atom, ne_atom, cz_atom)
|
||||||
|
h2_atom.name = "HN1"
|
||||||
|
h3_atom = protonate_direction(nh1_atom, ne_atom, cd_atom)
|
||||||
|
h3_atom.name = "HN2"
|
||||||
|
h4_atom = protonate_direction(nh2_atom, ne_atom, cz_atom)
|
||||||
|
h4_atom.name = "HN3"
|
||||||
|
h5_atom = protonate_direction(nh2_atom, ne_atom, h1_atom)
|
||||||
|
h5_atom.name = "HN4"
|
||||||
|
return [h1_atom, h2_atom, h3_atom, h4_atom, h5_atom]
|
||||||
|
|
||||||
|
|
||||||
|
def add_his_hydrogen(residue):
|
||||||
|
"""Adds His hydrogen atoms to residues according to the 'old way'.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
residue: histidine residue to protonate
|
||||||
|
"""
|
||||||
|
for atom in residue:
|
||||||
|
if atom.name == "CG":
|
||||||
|
cg_atom = atom
|
||||||
|
elif atom.name == "ND1":
|
||||||
|
nd_atom = atom
|
||||||
|
elif atom.name == "CD2":
|
||||||
|
cd_atom = atom
|
||||||
|
elif atom.name == "CE1":
|
||||||
|
ce_atom = atom
|
||||||
|
elif atom.name == "NE2":
|
||||||
|
ne_atom = atom
|
||||||
|
hd_atom = protonate_sp2(cg_atom, nd_atom, ce_atom)
|
||||||
|
hd_atom.name = "HND"
|
||||||
|
he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom)
|
||||||
|
he_atom.name = "HNE"
|
||||||
|
|
||||||
|
|
||||||
|
def add_trp_hydrogen(residue):
|
||||||
|
"""Adds Trp hydrogen atoms to residues according to the 'old way'.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
residue: tryptophan residue to protonate
|
||||||
|
"""
|
||||||
|
cd_atom = None
|
||||||
|
ne_atom = None
|
||||||
|
for atom in residue:
|
||||||
|
if atom.name == "CD1":
|
||||||
|
cd_atom = atom
|
||||||
|
elif atom.name == "NE1":
|
||||||
|
ne_atom = atom
|
||||||
|
elif atom.name == "CE2":
|
||||||
|
ce_atom = atom
|
||||||
|
if (cd_atom is None) or (ne_atom is None) or (ce_atom is None):
|
||||||
|
str_ = "Unable to find all atoms for {0:s} {1:s}".format(
|
||||||
|
residue[0].res_name, residue[0].res_num)
|
||||||
|
raise ValueError(str_)
|
||||||
|
he_atom = protonate_sp2(cd_atom, ne_atom, ce_atom)
|
||||||
|
he_atom.name = "HNE"
|
||||||
|
|
||||||
|
|
||||||
|
def add_amd_hydrogen(residue):
|
||||||
|
"""Adds Gln & Asn hydrogen atoms to residues according to the 'old way'.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
residue: glutamine or asparagine residue to protonate
|
||||||
|
"""
|
||||||
|
c_atom = None
|
||||||
|
o_atom = None
|
||||||
|
n_atom = None
|
||||||
|
for atom in residue:
|
||||||
|
if ((atom.res_name == "GLN" and atom.name == "CD")
|
||||||
|
or (atom.res_name == "ASN" and atom.name == "CG")):
|
||||||
|
c_atom = atom
|
||||||
|
elif ((atom.res_name == "GLN" and atom.name == "OE1")
|
||||||
|
or (atom.res_name == "ASN" and atom.name == "OD1")):
|
||||||
|
o_atom = atom
|
||||||
|
elif ((atom.res_name == "GLN" and atom.name == "NE2")
|
||||||
|
or (atom.res_name == "ASN" and atom.name == "ND2")):
|
||||||
|
n_atom = atom
|
||||||
|
if (c_atom is None) or (o_atom is None) or (n_atom is None):
|
||||||
|
str_ = "Unable to find all atoms for {0:s} {1:s}".format(
|
||||||
|
residue[0].res_name, residue[0].res_num)
|
||||||
|
raise ValueError(str_)
|
||||||
|
h1_atom = protonate_direction(n_atom, o_atom, c_atom)
|
||||||
|
h1_atom.name = "HN1"
|
||||||
|
h2_atom = protonate_average_direction(n_atom, c_atom, o_atom)
|
||||||
|
h2_atom.name = "HN2"
|
||||||
|
|
||||||
|
|
||||||
|
def add_backbone_hydrogen(residue, o_atom, c_atom):
|
||||||
|
"""Adds hydrogen backbone atoms to residues according to the old way.
|
||||||
|
|
||||||
|
dR is wrong for the N-terminus (i.e. first residue) but it doesn't affect
|
||||||
|
anything at the moment. Could be improved, but works for now.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
residue: residue to protonate
|
||||||
|
o_atom: backbone oxygen atom
|
||||||
|
c_atom: backbone carbon atom
|
||||||
|
Returns:
|
||||||
|
[new backbone oxygen atom, new backbone carbon atom]
|
||||||
|
"""
|
||||||
|
new_c_atom = None
|
||||||
|
new_o_atom = None
|
||||||
|
n_atom = None
|
||||||
|
for atom in residue:
|
||||||
|
if atom.name == "N":
|
||||||
|
n_atom = atom
|
||||||
|
if atom.name == "C":
|
||||||
|
new_c_atom = atom
|
||||||
|
if atom.name == "O":
|
||||||
|
new_o_atom = atom
|
||||||
|
if None in [c_atom, o_atom, n_atom]:
|
||||||
|
return [new_o_atom, new_c_atom]
|
||||||
|
if n_atom.res_name == "PRO":
|
||||||
|
# PRO doesn't have an H-atom; do nothing
|
||||||
|
pass
|
||||||
|
else:
|
||||||
|
h_atom = protonate_direction(n_atom, o_atom, c_atom)
|
||||||
|
h_atom.name = "H"
|
||||||
|
return [new_o_atom, new_c_atom]
|
||||||
|
|
||||||
|
|
||||||
|
def protonate_direction(x1_atom, x2_atom, x3_atom):
|
||||||
|
"""Protonates an atom, x1_atom, given a direction.
|
||||||
|
|
||||||
|
New direction for x1_atom proton is (x2_atom -> x3_atom).
|
||||||
|
|
||||||
|
Args:
|
||||||
|
x1_atom: atom to be protonated
|
||||||
|
x2_atom: atom for direction
|
||||||
|
x3_atom: other atom for direction
|
||||||
|
Returns:
|
||||||
|
new hydrogen atom
|
||||||
|
"""
|
||||||
|
dx = (x3_atom.x - x2_atom.x)
|
||||||
|
dy = (x3_atom.y - x2_atom.y)
|
||||||
|
dz = (x3_atom.z - x2_atom.z)
|
||||||
|
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
||||||
|
x = x1_atom.x + dx/length
|
||||||
|
y = x1_atom.y + dy/length
|
||||||
|
z = x1_atom.z + dz/length
|
||||||
|
h_atom = make_new_h(x1_atom, x, y, z)
|
||||||
|
h_atom.name = "H"
|
||||||
|
return h_atom
|
||||||
|
|
||||||
|
|
||||||
|
def protonate_average_direction(x1_atom, x2_atom, x3_atom):
|
||||||
|
"""Protonates an atom, x1_atom, given a direction.
|
||||||
|
|
||||||
|
New direction for x1_atom is (x1_atom/x2_atom -> x3_atom).
|
||||||
|
Note, this one uses the average of x1_atom & x2_atom (N & O) unlike
|
||||||
|
the previous N - C = O
|
||||||
|
|
||||||
|
Args:
|
||||||
|
x1_atom: atom to be protonated
|
||||||
|
x2_atom: atom for direction
|
||||||
|
x3_atom: other atom for direction
|
||||||
|
Returns:
|
||||||
|
new hydrogen atom
|
||||||
|
"""
|
||||||
|
dx = (x3_atom.x + x1_atom.x)*0.5 - x2_atom.x
|
||||||
|
dy = (x3_atom.y + x1_atom.y)*0.5 - x2_atom.y
|
||||||
|
dz = (x3_atom.z + x1_atom.z)*0.5 - x2_atom.z
|
||||||
|
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
||||||
|
x = x1_atom.x + dx/length
|
||||||
|
y = x1_atom.y + dy/length
|
||||||
|
z = x1_atom.z + dz/length
|
||||||
|
h_atom = make_new_h(x1_atom, x, y, z)
|
||||||
|
h_atom.name = "H"
|
||||||
|
return h_atom
|
||||||
|
|
||||||
|
|
||||||
|
def protonate_sp2(x1_atom, x2_atom, x3_atom):
|
||||||
|
"""Protonates a SP2 atom, given a list of atoms
|
||||||
|
|
||||||
|
Args:
|
||||||
|
x1_atom: atom to set direction
|
||||||
|
x2_atom: atom to be protonated
|
||||||
|
x3_atom: other atom to set direction
|
||||||
|
Returns:
|
||||||
|
new hydrogen atom
|
||||||
|
"""
|
||||||
|
dx = (x1_atom.x + x3_atom.x)*0.5 - x2_atom.x
|
||||||
|
dy = (x1_atom.y + x3_atom.y)*0.5 - x2_atom.y
|
||||||
|
dz = (x1_atom.z + x3_atom.z)*0.5 - x2_atom.z
|
||||||
|
length = math.sqrt(dx*dx + dy*dy + dz*dz)
|
||||||
|
x = x2_atom.x - dx/length
|
||||||
|
y = x2_atom.y - dy/length
|
||||||
|
z = x2_atom.z - dz/length
|
||||||
|
h_atom = make_new_h(x2_atom, x, y, z)
|
||||||
|
h_atom.name = "H"
|
||||||
|
return h_atom
|
||||||
|
|
||||||
|
|
||||||
|
def make_new_h(atom, x, y, z):
|
||||||
|
"""Add a new hydrogen to an atom at the specified position.
|
||||||
|
|
||||||
|
Args:
|
||||||
|
atom: atom to protonate
|
||||||
|
x: x position of hydrogen
|
||||||
|
y: y position of hydrogen
|
||||||
|
z: z position of hydrogen
|
||||||
|
Returns:
|
||||||
|
new hydrogen atom
|
||||||
|
"""
|
||||||
|
new_h = Atom()
|
||||||
|
new_h.set_property(
|
||||||
|
numb=None, name='H{0:s}'.format(atom.name[1:]),
|
||||||
|
res_name=atom.res_name, chain_id=atom.chain_id,
|
||||||
|
res_num=atom.res_num, x=x, y=y, z=z, occ=None, beta=None)
|
||||||
|
new_h.element = 'H'
|
||||||
|
new_h.bonded_atoms = [atom]
|
||||||
|
new_h.charge = 0
|
||||||
|
new_h.steric_number = 0
|
||||||
|
new_h.number_of_lone_pairs = 0
|
||||||
|
new_h.number_of_protons_to_add = 0
|
||||||
|
new_h.num_pi_elec_2_3_bonds = 0
|
||||||
|
atom.bonded_atoms.append(new_h)
|
||||||
|
atom.conformation_container.add_atom(new_h)
|
||||||
|
return new_h
|
||||||
|
|
||||||
|
|
||||||
@@ -3,7 +3,13 @@
|
|||||||
TODO - this module unnecessarily confuses the code. Can we eliminate it?
|
TODO - this module unnecessarily confuses the code. Can we eliminate it?
|
||||||
"""
|
"""
|
||||||
from propka.lib import info
|
from propka.lib import info
|
||||||
import propka.calculations as calcs
|
from propka.hydrogens import setup_bonding_and_protonation, setup_bonding
|
||||||
|
from propka.hydrogens import setup_bonding_and_protonation_30_style
|
||||||
|
from propka.energy import radial_volume_desolvation, calculate_pair_weight
|
||||||
|
from propka.energy import hydrogen_bond_energy, hydrogen_bond_interaction
|
||||||
|
from propka.energy import electrostatic_interaction, check_coulomb_pair
|
||||||
|
from propka.energy import coulomb_energy, check_exceptions
|
||||||
|
from propka.energy import backbone_reorganization
|
||||||
|
|
||||||
|
|
||||||
class Version:
|
class Version:
|
||||||
@@ -79,8 +85,7 @@ class Version:
|
|||||||
|
|
||||||
def setup_bonding_and_protonation(self, molecular_container):
|
def setup_bonding_and_protonation(self, molecular_container):
|
||||||
"""Setup bonding and protonation using assigned model."""
|
"""Setup bonding and protonation using assigned model."""
|
||||||
return self.molecular_preparation_method(
|
return self.molecular_preparation_method(molecular_container)
|
||||||
self.parameters, molecular_container)
|
|
||||||
|
|
||||||
def setup_bonding(self, molecular_container):
|
def setup_bonding(self, molecular_container):
|
||||||
"""Setup bonding using assigned model."""
|
"""Setup bonding using assigned model."""
|
||||||
@@ -94,18 +99,18 @@ class VersionA(Version):
|
|||||||
"""Initialize object with parameters."""
|
"""Initialize object with parameters."""
|
||||||
# set the calculation rutines used in this version
|
# set the calculation rutines used in this version
|
||||||
super().__init__(parameters)
|
super().__init__(parameters)
|
||||||
self.molecular_preparation_method = calcs.setup_bonding_and_protonation
|
self.molecular_preparation_method = setup_bonding_and_protonation
|
||||||
self.prepare_bonds = calcs.setup_bonding
|
self.prepare_bonds = setup_bonding
|
||||||
self.desolvation_model = calcs.radial_volume_desolvation
|
self.desolvation_model = radial_volume_desolvation
|
||||||
self.weight_pair_method = calcs.calculate_pair_weight
|
self.weight_pair_method = calculate_pair_weight
|
||||||
self.sidechain_interaction_model = calcs.hydrogen_bond_energy
|
self.sidechain_interaction_model = hydrogen_bond_energy
|
||||||
self.hydrogen_bond_interaction_model = calcs.hydrogen_bond_interaction
|
self.hydrogen_bond_interaction_model = hydrogen_bond_interaction
|
||||||
self.electrostatic_interaction_model = calcs.electrostatic_interaction
|
self.electrostatic_interaction_model = electrostatic_interaction
|
||||||
self.check_coulomb_pair_method = calcs.check_coulomb_pair
|
self.check_coulomb_pair_method = check_coulomb_pair
|
||||||
self.coulomb_interaction_model = calcs.coulomb_energy
|
self.coulomb_interaction_model = coulomb_energy
|
||||||
self.backbone_interaction_model = calcs.hydrogen_bond_energy
|
self.backbone_interaction_model = hydrogen_bond_energy
|
||||||
self.backbone_reorganisation_method = calcs.backbone_reorganization
|
self.backbone_reorganisation_method = backbone_reorganization
|
||||||
self.exception_check_method = calcs.check_exceptions
|
self.exception_check_method = check_exceptions
|
||||||
|
|
||||||
def get_hydrogen_bond_parameters(self, atom1, atom2):
|
def get_hydrogen_bond_parameters(self, atom1, atom2):
|
||||||
"""Get hydrogen bond parameters for two atoms.
|
"""Get hydrogen bond parameters for two atoms.
|
||||||
@@ -265,14 +270,14 @@ class Propka30(Version):
|
|||||||
# set the calculation routines used in this version
|
# set the calculation routines used in this version
|
||||||
super().__init__(parameters)
|
super().__init__(parameters)
|
||||||
self.molecular_preparation_method = (
|
self.molecular_preparation_method = (
|
||||||
calcs.setup_bonding_and_protonation_30_style)
|
setup_bonding_and_protonation_30_style)
|
||||||
self.desolvation_model = calcs.radial_volume_desolvation
|
self.desolvation_model = radial_volume_desolvation
|
||||||
self.weight_pair_method = calcs.calculate_pair_weight
|
self.weight_pair_method = calculate_pair_weight
|
||||||
self.sidechain_interaction_model = calcs.hydrogen_bond_energy
|
self.sidechain_interaction_model = hydrogen_bond_energy
|
||||||
self.check_coulomb_pair_method = calcs.check_coulomb_pair
|
self.check_coulomb_pair_method = check_coulomb_pair
|
||||||
self.coulomb_interaction_model = calcs.coulomb_energy
|
self.coulomb_interaction_model = coulomb_energy
|
||||||
self.backbone_reorganisation_method = calcs.backbone_reorganization
|
self.backbone_reorganisation_method = backbone_reorganization
|
||||||
self.exception_check_method = calcs.check_exceptions
|
self.exception_check_method = check_exceptions
|
||||||
|
|
||||||
def get_hydrogen_bond_parameters(self, atom1, atom2):
|
def get_hydrogen_bond_parameters(self, atom1, atom2):
|
||||||
"""Get hydrogen bond parameters for two atoms.
|
"""Get hydrogen bond parameters for two atoms.
|
||||||
|
|||||||
Reference in New Issue
Block a user