diff --git a/propka/calculations.py b/propka/calculations.py index ea2a54b..5bf2449 100644 --- a/propka/calculations.py +++ b/propka/calculations.py @@ -636,7 +636,7 @@ def check_coulomb_pair(parameters, group1, group2, dist): num_volume = group1.Nmass + group2.Nmass do_coulomb = True # check if both groups are titratable (ions are taken care of in - # determinants::setIonDeterminants) + # determinants::set_ion_determinants) if not (group1.titratable and group2.titratable): do_coulomb = False # check if the distance is not too big diff --git a/propka/conformation_container.py b/propka/conformation_container.py index 7b2f370..107a0a8 100644 --- a/propka/conformation_container.py +++ b/propka/conformation_container.py @@ -4,8 +4,8 @@ import propka.ligand from propka.output import make_interaction_map from propka.determinant import Determinant from propka.coupled_groups import NCCG -from propka.determinants import setBackBoneDeterminants, setIonDeterminants -from propka.determinants import setDeterminants +from propka.determinants import set_backbone_determinants, set_ion_determinants +from propka.determinants import set_determinants from propka.group import Group, is_group from propka.lib import info @@ -190,15 +190,15 @@ class ConformationContainer: for group in self.get_titratable_groups() + self.get_ions(): version.calculate_desolvation(group) # calculate backbone interactions - setBackBoneDeterminants(self.get_titratable_groups(), + set_backbone_determinants(self.get_titratable_groups(), self.get_backbone_groups(), version) # setting ion determinants - setIonDeterminants(self, version) + set_ion_determinants(self, version) # calculating the back-bone reorganization/desolvation term version.calculatebackbone_reorganization(self) # setting remaining non-iterative and iterative side-chain & Coulomb # interaction determinants - setDeterminants(self.get_sidechain_groups(), version=version, + set_determinants(self.get_sidechain_groups(), version=version, options=options) # calculating the total pKa values for group in self.groups: diff --git a/propka/determinants.py b/propka/determinants.py index 52da557..7b75d3f 100644 --- a/propka/determinants.py +++ b/propka/determinants.py @@ -1,23 +1,33 @@ +"""Functions to manipulate Determinant objects. -from __future__ import division -from __future__ import print_function - -import math, time - -import propka.iterative, propka.lib, propka.vector_algebra -import propka.calculations +TODO - it is confusing to have both `determinant.py` and `determinants.py`. +Should these be merged? +""" +import math +import propka.iterative +import propka.lib +import propka.vector_algebra +from propka.calculations import squared_distance, get_smallest_distance +from propka.calculations import angle_distance_factors, hydrogen_bond_energy from propka.determinant import Determinant -# TODO - it is confusing to have both `determinant.py` and `determinants.py`. -# Should these be merged? +# Cutoff for angle factor +# TODO - this constant appears elsewhere in the package. +# It should be moved to a configuration file. +FANGLE_MIN = 0.001 -def setDeterminants(propka_groups, version=None, options=None): +def set_determinants(propka_groups, version=None, options=None): + """Add side-chain and coulomb determinants/perturbations to all residues. + + NOTE - backbone determinants are set separately + + Args: + propka_groups: groups to adjust + version: version object + options: options object """ - adding side-chain and coulomb determinants/perturbations to all residues - note, backbone determinants are set separately - """ - iterative_interactions = [] # --- NonIterative section ---# for group1 in propka_groups: @@ -27,172 +37,198 @@ def setDeterminants(propka_groups, version=None, options=None): # do not calculate interactions for coupled groups if group2 in group1.covalently_coupled_groups: break - distance = propka.calculations.distance(group1, group2) - if distance < version.parameters.coulomb_cutoff2: - interaction_type = version.parameters.interaction_matrix.get_value(group1.type,group2.type) + interaction_type = version.parameters.interaction_matrix.get_value(group1.type, + group2.type) if interaction_type == 'I': - propka.iterative.addtoDeterminantList(group1, group2, distance, iterative_interactions, version=version) + propka.iterative.addtoDeterminantList(group1, group2, + distance, + iterative_interactions, + version=version) elif interaction_type == 'N': - addDeterminants(group1, group2, distance, version) - - + add_determinants(group1, group2, distance, version) # --- Iterative section ---# - propka.iterative.addDeterminants(iterative_interactions, version, options=options) + propka.iterative.add_determinants(iterative_interactions, version, + options=options) -def addDeterminants(group1, group2, distance, version): +def add_determinants(group1, group2, distance, version): + """Add determinants and perturbations for distance(R1, R2) < coulomb_cutoff. + + Args: + group1: first group to add + group2: second group to add + distance: distance between groups + version: version object """ - adding determinants/perturbations, distance(R1, R2) < coulomb_cutoff always - """ - # side-chain determinant - addSidechainDeterminants(group1, group2, version) - + add_sidechain_determinants(group1, group2, version) # Coulomb determinant - addCoulombDeterminants(group1, group2, distance, version) + add_coulomb_determinants(group1, group2, distance, version) - return -def addSidechainDeterminants(group1, group2, version=None): +def add_sidechain_determinants(group1, group2, version=None): + """Add side-chain determinants and perturbations. + + NOTE - res_num1 > res_num2 + + Args: + group1: first group to add + group2: second group to add + version: version object """ - adding side-chain determinants/perturbations - Note, res_num1 > res_num2 - """ - hbond_interaction = version.hydrogen_bond_interaction(group1, group2) - if hbond_interaction: - if group1.charge == group2.charge: # acid pair or base pair if group1.model_pka < group2.model_pka: - newDeterminant1 = Determinant(group2, -hbond_interaction) - newDeterminant2 = Determinant(group1, hbond_interaction) + new_determinant1 = Determinant(group2, -hbond_interaction) + new_determinant2 = Determinant(group1, hbond_interaction) else: - newDeterminant1 = Determinant(group2, hbond_interaction) - newDeterminant2 = Determinant(group1, -hbond_interaction) + new_determinant1 = Determinant(group2, hbond_interaction) + new_determinant2 = Determinant(group1, -hbond_interaction) else: - newDeterminant1 = Determinant(group2, hbond_interaction*group1.charge) - newDeterminant2 = Determinant(group1, hbond_interaction*group2.charge) + new_determinant1 = Determinant(group2, hbond_interaction*group1.charge) + new_determinant2 = Determinant(group1, hbond_interaction*group2.charge) + group1.determinants['sidechain'].append(new_determinant1) + group2.determinants['sidechain'].append(new_determinant2) - group1.determinants['sidechain'].append(newDeterminant1) - group2.determinants['sidechain'].append(newDeterminant2) - return +def add_coulomb_determinants(group1, group2, distance, version): + """Add non-iterative Coulomb determinants and perturbations. -def addCoulombDeterminants(group1, group2, distance, version): + Args: + group1: first group to add + group2: second group to add + distance: distance between groups + version: version object """ - adding NonIterative Coulomb determinants/perturbations - """ - - coulomb_interaction = version.electrostatic_interaction(group1, group2, distance) - + coulomb_interaction = version.electrostatic_interaction(group1, group2, + distance) if coulomb_interaction: - Q1 = group1.charge - Q2 = group2.charge - + q1 = group1.charge + q2 = group2.charge # assigning the Coulombic interaction - if Q1 < 0.0 and Q2 < 0.0: - """ both are acids """ - addCoulombAcidPair(group1, group2, coulomb_interaction) - elif Q1 > 0.0 and Q2 > 0.0: - """ both are bases """ - addCoulombBasePair(group1, group2, coulomb_interaction) + if q1 < 0.0 and q2 < 0.0: + # both are acids + add_coulomb_acid_pair(group1, group2, coulomb_interaction) + elif q1 > 0.0 and q2 > 0.0: + # both are bases + add_coulomb_base_pair(group1, group2, coulomb_interaction) else: - """ one of each """ - addCoulombIonPair(group1, group2, coulomb_interaction) - - return + # one of each + add_coulomb_ion_pair(group1, group2, coulomb_interaction) -def addCoulombAcidPair(object1, object2, value): +def add_coulomb_acid_pair(object1, object2, value): + """Add the Coulomb interaction (an acid pair). + + The higher pKa is raised. + + Args: + object1: first part of pair + object2: second part of pair + value: determinant value """ - Adding the Coulomb interaction (an acid pair): - the higher pKa is raised - """ - if object1.model_pka > object2.model_pka: - newDeterminant = Determinant(object2, value) - object1.determinants['coulomb'].append(newDeterminant) + new_determinant = Determinant(object2, value) + object1.determinants['coulomb'].append(new_determinant) else: - newDeterminant = Determinant(object1, value) - object2.determinants['coulomb'].append(newDeterminant) + new_determinant = Determinant(object1, value) + object2.determinants['coulomb'].append(new_determinant) -def addCoulombBasePair(object1, object2, value): - """ - Adding the Coulomb interaction (a base pair): - the lower pKa is lowered +def add_coulomb_base_pair(object1, object2, value): + """Add the Coulomb interaction (a base pair). + + The lower pKa is lowered. + + Args: + object1: first part of pair + object2: second part of pair + value: determinant value """ if object1.model_pka < object2.model_pka: - newDeterminant = Determinant(object2, -value) - object1.determinants['coulomb'].append(newDeterminant) + new_determinant = Determinant(object2, -value) + object1.determinants['coulomb'].append(new_determinant) else: - newDeterminant = Determinant(object1, -value) - object2.determinants['coulomb'].append(newDeterminant) + new_determinant = Determinant(object1, -value) + object2.determinants['coulomb'].append(new_determinant) -def addCoulombIonPair(object1, object2, value): +def add_coulomb_ion_pair(object1, object2, value): + """Add the Coulomb interaction (an acid-base pair). + + The pKa of the acid is lowered & the pKa of the base is raised. + + Args: + object1: first part of pair + object2: second part of pair + value: determinant value """ - Adding the Coulomb interaction (an acid-base pair): - the pKa of the acid is lowered & the pKa of the base is raised - """ - # residue1 - Q1 = object1.charge - newDeterminant = Determinant(object2, Q1*value) - object1.determinants['coulomb'].append(newDeterminant) - + q1 = object1.charge + new_determinant = Determinant(object2, q1*value) + object1.determinants['coulomb'].append(new_determinant) # residue2 - Q2 = object2.charge - newDeterminant = Determinant(object1, Q2*value) - object2.determinants['coulomb'].append(newDeterminant) + q2 = object2.charge + new_determinant = Determinant(object1, q2*value) + object2.determinants['coulomb'].append(new_determinant) +def set_ion_determinants(conformation_container, version): + """Add ion determinants and perturbations. - -def setIonDeterminants(conformation_container, version): - """ - adding ion determinants/perturbations + Args: + conformation_container: conformation to set + version: version object """ for titratable_group in conformation_container.get_titratable_groups(): for ion_group in conformation_container.get_ions(): - squared_distance = propka.calculations.squared_distance(titratable_group, ion_group) - if squared_distance < version.parameters.coulomb_cutoff2_squared: - weight = version.calculate_pair_weight(titratable_group.Nmass, ion_group.Nmass) - # the pKa of both acids and bases are shifted up by negative ions (and vice versa) - value = (-ion_group.charge) * version.calculate_coulomb_energy(math.sqrt(squared_distance), weight) - newDeterminant = Determinant(ion_group, value) - titratable_group.determinants['coulomb'].append(newDeterminant) + dist_sq = squared_distance(titratable_group, ion_group) + if dist_sq < version.parameters.coulomb_cutoff2_squared: + weight = version.calculate_pair_weight(titratable_group.Nmass, + ion_group.Nmass) + # the pKa of both acids and bases are shifted up by negative + # ions (and vice versa) + value = (-ion_group.charge) \ + * version.calculate_coulomb_energy(math.sqrt(dist_sq), + weight) + new_det = Determinant(ion_group, value) + titratable_group.determinants['coulomb'].append(new_det) - return -def setBackBoneDeterminants(titratable_groups, backbone_groups, version): +def set_backbone_determinants(titratable_groups, backbone_groups, version): + """Set determinants between titrable and backbone groups. + Args: + titratable_groups: list of titratable groups + backbone_groups: list of backbone groups + version: version object + """ for titratable_group in titratable_groups: - titratable_group_interaction_atoms = titratable_group.interaction_atoms_for_acids + titratable_group_interaction_atoms \ + = titratable_group.interaction_atoms_for_acids if not titratable_group_interaction_atoms: continue - # find out which backbone groups this titratable is interacting with for backbone_group in backbone_groups: # find the interacting atoms - backbone_interaction_atoms = backbone_group.get_interaction_atoms(titratable_group) + backbone_interaction_atoms \ + = backbone_group.get_interaction_atoms(titratable_group) if not backbone_interaction_atoms: continue - # find the smallest distance - [backbone_atom, distance, titratable_atom] = propka.calculations.get_smallest_distance(backbone_interaction_atoms, - titratable_group_interaction_atoms) + [backbone_atom, distance, titratable_atom] \ + = get_smallest_distance(backbone_interaction_atoms, \ + titratable_group_interaction_atoms) # get the parameters - parameters = version.get_backbone_hydrogen_bond_parameters(backbone_atom, titratable_atom) + parameters = version.get_backbone_hydrogen_bond_parameters(backbone_atom, + titratable_atom) if not parameters: continue - [dpKa_max, [cutoff1, cutoff2]] = parameters - - + [dpka_max, [cutoff1, cutoff2]] = parameters if distance < cutoff2: # calculate angle factor f_angle = 1.0 @@ -206,19 +242,20 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): # || # C if backbone_group.type == 'BBC': - if titratable_group.type in version.parameters.angular_dependent_sidechain_interactions: + if titratable_group.type \ + in version.parameters.angular_dependent_sidechain_interactions: if titratable_atom.element == 'H': - heavy_atom = titratable_atom.bonded_atoms[0] + heavy_atom = titratable_atom.bonded_atoms[0] hydrogen_atom = titratable_atom - [d1, f_angle, d2] = propka.calculations.angle_distance_factors(atom1=heavy_atom, - atom2=hydrogen_atom, - atom3=backbone_atom) + [_, f_angle, _] = angle_distance_factors(atom1=heavy_atom, + atom2=hydrogen_atom, + atom3=backbone_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 + # 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 - # for BBN groups, the hydrogen is on the backbone group # # Titra. @@ -229,23 +266,21 @@ def setBackBoneDeterminants(titratable_groups, backbone_groups, version): # / \ if backbone_group.type == 'BBN': if backbone_atom.element == 'H': - backbone_N = backbone_atom.bonded_atoms[0] - backbone_H = backbone_atom - [d1, f_angle, d2] = propka.calculations.angle_distance_factors(atom1=titratable_atom, - atom2=backbone_H, - atom3=backbone_N) + backbone_n = backbone_atom.bonded_atoms[0] + backbone_h = backbone_atom + [_, f_angle, _] = angle_distance_factors(atom1=titratable_atom, + atom2=backbone_h, + atom3=backbone_n) 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 + # 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 - - - if f_angle > 0.001: - value = titratable_group.charge * propka.calculations.hydrogen_bond_energy(distance, dpKa_max, [cutoff1,cutoff2], f_angle) - - newDeterminant = Determinant(backbone_group, value) - titratable_group.determinants['backbone'].append(newDeterminant) - - - return + if f_angle > FANGLE_MIN: + value = titratable_group.charge * hydrogen_bond_energy(distance, + dpka_max, + [cutoff1, cutoff2], + f_angle) + new_determinant = Determinant(backbone_group, value) + titratable_group.determinants['backbone'].append(new_determinant) + \ No newline at end of file diff --git a/propka/iterative.py b/propka/iterative.py index 3614d93..132d42b 100644 --- a/propka/iterative.py +++ b/propka/iterative.py @@ -166,7 +166,7 @@ def addIterativeIonPair(object1, object2, interaction, version): object2.determinants['sidechain'].append(interaction) -def addDeterminants(iterative_interactions, version, options=None): +def add_determinants(iterative_interactions, version, options=None): """ The iterative pKa scheme. Later it is all added in 'calculateTotalPKA' """