Source code for propka.iterative

"""
Working with Determinants
=========================

Iterative functions for pKa calculations. These appear to mostly
involve :class:`propka.determinant.Determinant` instances.

"""
import logging
from typing import Dict, Iterable, List, Optional, Sequence, Tuple
from propka.determinant import Determinant
from propka.group import Group
from propka.version import Version


_LOGGER = logging.getLogger(__name__)


# TODO - these are undocumented constants
UNK_MIN_VALUE = 0.005

Interaction = list


[docs] def add_to_determinant_list(group1: Group, group2: Group, distance: float, iterative_interactions: List[Interaction], version: Version): """Add iterative determinantes to the list. [[R1, R2], [side-chain, coulomb], [A1, A2]], ... NOTE - sign is determined when the interaction is added to the iterative object! NOTE - distance < coulomb_cutoff here Args: group1: first group in pair group2: second group in pair distance: distance between groups iterative_interactions: interaction list to modify version: version object """ hbond_value = version.hydrogen_bond_interaction(group1, group2) coulomb_value = version.electrostatic_interaction(group1, group2, distance) # adding the interaction to 'iterative_interactions' if hbond_value or coulomb_value: pair = [group1, group2] values = [hbond_value, coulomb_value] while None in values: values[values.index(None)] = 0.0 annihilation = [0., 0.] interaction = [pair, values, annihilation] iterative_interactions.append(interaction)
[docs] def add_iterative_acid_pair(object1: "Iterative", object2: "Iterative", interaction: Interaction): """Add the Coulomb 'iterative' interaction (an acid pair). The higher pKa is raised with QQ+HB The lower pKa is lowered with HB Args: object1: first object in pair object2: second object in pair interaction: list with [values, annihilation] """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] coulomb_value = values[1] diff = coulomb_value + 2*hbond_value comp1 = object1.pka_old + annihilation[0] + diff comp2 = object2.pka_old + annihilation[1] + diff annihilation[0] = 0.0 annihilation[1] = 0.0 if comp1 > comp2: # side-chain determinant = [object2, hbond_value] object1.determinants['sidechain'].append(determinant) determinant = [object1, -hbond_value] object2.determinants['sidechain'].append(determinant) # Coulomb determinant = [object2, coulomb_value] object1.determinants['coulomb'].append(determinant) annihilation[0] = -diff else: # side-chain determinant = [object1, hbond_value] object2.determinants['sidechain'].append(determinant) determinant = [object2, -hbond_value] object1.determinants['sidechain'].append(determinant) # Coulomb determinant = [object1, coulomb_value] object2.determinants['coulomb'].append(determinant) annihilation[1] = -diff
[docs] def add_iterative_base_pair(object1: "Iterative", object2: "Iterative", interaction: Interaction): """Add the Coulomb 'iterative' interaction (a base pair). The lower pKa is lowered Args: object1: first object in pair object2: second object in pair interaction: list with [values, annihilation] """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] coulomb_value = values[1] diff = coulomb_value + 2*hbond_value diff = -diff comp1 = object1.pka_old + annihilation[0] + diff comp2 = object2.pka_old + annihilation[1] + diff annihilation[0] = 0.0 annihilation[1] = 0.0 if comp1 < comp2: # side-chain determinant = [object2, -hbond_value] object1.determinants['sidechain'].append(determinant) determinant = [object1, hbond_value] object2.determinants['sidechain'].append(determinant) # Coulomb determinant = [object2, -coulomb_value] object1.determinants['coulomb'].append(determinant) annihilation[0] = -diff else: # side-chain determinant = [object1, -hbond_value] object2.determinants['sidechain'].append(determinant) determinant = [object2, hbond_value] object1.determinants['sidechain'].append(determinant) # Coulomb determinant = [object1, -coulomb_value] object2.determinants['coulomb'].append(determinant) annihilation[1] = -diff
[docs] def add_iterative_ion_pair(object1: "Iterative", object2: "Iterative", interaction: Interaction, version: Version): """Add the Coulomb 'iterative' interaction (an acid-base pair) the pKa of the acid is lowered & the pKa of the base is raised Args: object1: first object in pair object2: second object in pair interaction: list with [values, annihilation] version: version object """ values = interaction[1] annihilation = interaction[2] hbond_value = values[0] coulomb_value = values[1] q1 = object1.q q2 = object2.q comp1 = object1.pka_old + annihilation[0] + q1*coulomb_value comp2 = object2.pka_old + annihilation[1] + q2*coulomb_value if (object1.res_name not in version.parameters.exclude_sidechain_interactions): comp1 += q1*hbond_value if (object2.res_name not in version.parameters.exclude_sidechain_interactions): comp2 += q2*hbond_value if q1 == -1.0 and comp1 < comp2: # pKa(acid) < pKa(base) add_term = True elif q1 == 1.0 and comp1 > comp2: # pKa(base) > pKa(acid) add_term = True else: add_term = False annihilation[0] = 0.00 annihilation[1] = 0.00 if add_term: # Coulomb if coulomb_value > UNK_MIN_VALUE: # residue1 interaction = [object2, q1*coulomb_value] annihilation[0] += -q1*coulomb_value object1.determinants['coulomb'].append(interaction) # residue2 interaction = [object1, q2*coulomb_value] annihilation[1] += -q2*coulomb_value object2.determinants['coulomb'].append(interaction) # Side-chain if hbond_value > UNK_MIN_VALUE: # residue1 if (object1.res_name not in version.parameters.exclude_sidechain_interactions): interaction = [object2, q1*hbond_value] annihilation[0] += -q1*hbond_value object1.determinants['sidechain'].append(interaction) # residue2 if (object2.res_name not in version.parameters.exclude_sidechain_interactions): interaction = [object1, q2*hbond_value] annihilation[1] += -q2*hbond_value object2.determinants['sidechain'].append(interaction)
[docs] def add_determinants(iterative_interactions: List[Interaction], version: Version): """Add determinants iteratively. The iterative pKa scheme. Later it is all added in 'calculateTotalPKA' Args: iterative_interactions: list of iterative interactions version: version object _: options object """ # --- setup --- iteratives: List[Iterative] = [] done_group = [] # create iterative objects with references to their real group counterparts for interaction in iterative_interactions: pair = interaction[0] for group in pair: if group in done_group: # do nothing - already have an iterative object for this group pass else: new_iterative = Iterative(group) iteratives.append(new_iterative) done_group.append(group) # Initialize iterative scheme _LOGGER.debug( "\n --- pKa iterations ({0:d} groups, {1:d} interactions) " "---".format( len(iteratives), len(iterative_interactions) ) ) converged = False iteration = 0 # set non-iterative pka values as first step for iter_ in iteratives: iter_.pka_iter.append(iter_.pka_noniterative) # --- starting pKa iterations --- while not converged: # initialize pka_new iteration += 1 for itres in iteratives: itres.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []} itres.pka_new = itres.pka_noniterative # Adding interactions to temporary determinant container for interaction in iterative_interactions: pair = interaction[0] object1, object2 = find_iterative(pair, iteratives) q1 = object1.q q2 = object2.q if q1 < 0.0 and q2 < 0.0: # both are acids add_iterative_acid_pair(object1, object2, interaction) elif q1 > 0.0 and q2 > 0.0: # both are bases add_iterative_base_pair(object1, object2, interaction) else: # one of each add_iterative_ion_pair(object1, object2, interaction, version) # Calculating pka_new values for itres in iteratives: for type_ in ['sidechain', 'backbone', 'coulomb']: for determinant in itres.determinants[type_]: itres.pka_new += determinant[1] # Check convergence converged = True for itres in iteratives: if itres.pka_new == itres.pka_old: itres.converged = True else: itres.converged = False converged = False # reset pka_old & storing pka_new in pka_iter for itres in iteratives: assert itres.pka_new is not None itres.pka_old = itres.pka_new itres.pka_iter.append(itres.pka_new) if iteration == 10: _LOGGER.info( "did not converge in {0:d} iterations".format(iteration) ) break # printing pKa iterations # formerly was conditioned on if options.verbosity >= 2 - now unnecessary str_ = ' ' for index in range(iteration+1): str_ += "{0:>8d}".format(index) _LOGGER.debug(str_) for itres in iteratives: str_ = "{0:s} ".format(itres.label) for pka in itres.pka_iter: str_ += "{0:>8.2f}".format(pka) if not itres.converged: str_ += " *" _LOGGER.debug(str_) # creating real determinants and adding them to group object for itres in iteratives: for type_ in ['sidechain', 'backbone', 'coulomb']: for interaction in itres.determinants[type_]: value: float = interaction[1] if value > UNK_MIN_VALUE or value < -UNK_MIN_VALUE: group = interaction[0] new_det = Determinant(group, value) itres.group.determinants[type_].append(new_det)
[docs] def find_iterative( pair: Sequence[Group], iteratives: Iterable["Iterative"], ) -> Tuple["Iterative", "Iterative"]: """Find the 'iteratives' that correspond to the groups in 'pair'. Args: pair: groups to match iteratives: list of iteratives to search Returns: 1. first matched iterative 2. second matched iterative """ iterative0: Optional[Iterative] = None iterative1: Optional[Iterative] = None for iterative in iteratives: if iterative.group == pair[0]: iterative0 = iterative elif iterative.group == pair[1]: iterative1 = iterative if iterative0 is None or iterative1 is None: raise LookupError("iteratives not found") return iterative0, iterative1
[docs] class Iterative: """Iterative class - pKa values and references of iterative groups. NOTE - this class has a fake determinant list, true determinants are made after the iterations are finished. """ def __init__(self, group: Group): """Initialize object with group. Args: group: group to use for initialization. """ self.label = group.label self.atom = group.atom self.res_name = group.residue_type self.q = group.charge self.pka_old: Optional[float] = None self.pka_new: Optional[float] = None self.pka_iter: List[float] = [] self.pka_noniterative = 0.00 self.determinants: Dict[str, list] = { 'sidechain': [], 'backbone': [], 'coulomb': [] } self.group = group self.converged = True # Calculate the Non-Iterative part of pKa from the group object # Side chain side_chain = 0.00 for determinant in group.determinants['sidechain']: value = determinant.value side_chain += value # Back bone back_bone = 0.00 for determinant in group.determinants['backbone']: value = determinant.value back_bone += value # Coulomb coulomb = 0.00 for determinant in group.determinants['coulomb']: value = determinant.value coulomb += value self.pka_noniterative = group.model_pka self.pka_noniterative += group.energy_volume self.pka_noniterative += group.energy_local self.pka_noniterative += side_chain self.pka_noniterative += back_bone self.pka_noniterative += coulomb self.pka_old = self.pka_noniterative def __eq__(self, other) -> bool: """Needed to use objects in sets.""" assert isinstance(other, (Iterative, Group)), type(other) if self.atom.type == 'atom': # In case of protein atoms we trust the labels return self.label == other.label else: # For heterogene atoms we also need to check the residue number return ( self.label == other.label and self.atom.res_num == other.atom.res_num) def __hash__(self): """Needed to use objects in sets.""" return id(self)