diff --git a/vermouth/rcsu/go_structure_bias.py b/vermouth/rcsu/go_structure_bias.py index 619706442..d315b0fb0 100644 --- a/vermouth/rcsu/go_structure_bias.py +++ b/vermouth/rcsu/go_structure_bias.py @@ -22,6 +22,8 @@ from ..selectors import filter_minimal, select_backbone from ..gmx.topology import NonbondParam from .go_utils import get_go_type_from_attributes +from ..log_helpers import StyleAdapter, get_logger +LOGGER = StyleAdapter(get_logger(__name__)) class ComputeStructuralGoBias(Processor): """ @@ -115,8 +117,12 @@ def _chain_id_to_resnode(self, chain, resid): dict a dict matching the chain,resid to the self.res_graph node """ + if self.__chain_id_to_resnode: - return self.__chain_id_to_resnode[(chain, resid)] + if self.__chain_id_to_resnode.get((chain, resid), None) is not None: + return self.__chain_id_to_resnode[(chain, resid)] + else: + LOGGER.debug(stacklevel=5, msg='chain-resid pair not found in molecule') # for each residue collect the chain and residue in a dict # we use this later for identifying the residues from the @@ -130,7 +136,11 @@ def _chain_id_to_resnode(self, chain, resid): resid_key = self.res_graph.nodes[resnode].get('_old_resid') self.__chain_id_to_resnode[(chain_key, resid_key)] = resnode - return self.__chain_id_to_resnode[(chain, resid)] + if self.__chain_id_to_resnode.get((chain, resid), None) is not None: + return self.__chain_id_to_resnode[(chain, resid)] + else: + LOGGER.debug(stacklevel=5, msg='chain-resid pair not found in molecule') + def contact_selector(self, molecule): """ @@ -157,6 +167,7 @@ def contact_selector(self, molecule): # self.res_dist connected_pairs = dict(nx.all_pairs_shortest_path_length(self.res_graph, cutoff=self.res_dist)) + bad_chains_warning = False for contact in self.system.go_params["go_map"][0]: resIDA, chainA, resIDB, chainB = contact # identify the contact in the residue graph based on @@ -166,38 +177,44 @@ def contact_selector(self, molecule): # make sure that both residues are not connected # note: contacts should be symmetric so we only # check against one - if resB not in connected_pairs[resA]: - # now we lookup the backbone nodes within the residue contact - bb_node_A = next(filter_minimal(self.res_graph.nodes[resA]['graph'], select_backbone)) - bb_node_B = next(filter_minimal(self.res_graph.nodes[resB]['graph'], select_backbone)) - # compute the distance between bb-beads - dist = np.linalg.norm(molecule.nodes[bb_node_A]['position'] - - molecule.nodes[bb_node_B]['position']) - # verify that the distance between BB-beads satisfies the - # cut-off criteria - if self.cutoff_long > dist > self.cutoff_short: - atype_a = next(get_go_type_from_attributes(self.res_graph.nodes[resA]['graph'], - _old_resid=resIDA, - chain=chainA, - prefix=self.moltype)) - atype_b = next(get_go_type_from_attributes(self.res_graph.nodes[resB]['graph'], - _old_resid=resIDB, - chain=chainB, - prefix=self.moltype)) - # Check if symmetric contact has already been processed before - # and if so, we append the contact to the final symmetric contact matrix - # and add the exclusions. Else, we add to the full valid contact_matrix - # and continue searching. - if (atype_b, atype_a, dist) in contact_matrix: - # generate backbone-backbone exclusions - # perhaps one day can be its own function - excl = Interaction(atoms=(bb_node_A, bb_node_B), - parameters=[], meta={"group": "Go model exclusion"}) - molecule.interactions['exclusions'].append(excl) - symmetrical_matrix.append((atype_a, atype_b, dist)) - else: - contact_matrix.append((atype_a, atype_b, dist)) + if (resA is not None) and (resB is not None): + if resB not in connected_pairs[resA]: + # now we lookup the backbone nodes within the residue contact + bb_node_A = next(filter_minimal(self.res_graph.nodes[resA]['graph'], select_backbone)) + bb_node_B = next(filter_minimal(self.res_graph.nodes[resB]['graph'], select_backbone)) + # compute the distance between bb-beads + dist = np.linalg.norm(molecule.nodes[bb_node_A]['position'] - + molecule.nodes[bb_node_B]['position']) + # verify that the distance between BB-beads satisfies the + # cut-off criteria + if self.cutoff_long > dist > self.cutoff_short: + atype_a = next(get_go_type_from_attributes(self.res_graph.nodes[resA]['graph'], + _old_resid=resIDA, + chain=chainA, + prefix=self.moltype)) + atype_b = next(get_go_type_from_attributes(self.res_graph.nodes[resB]['graph'], + _old_resid=resIDB, + chain=chainB, + prefix=self.moltype)) + # Check if symmetric contact has already been processed before + # and if so, we append the contact to the final symmetric contact matrix + # and add the exclusions. Else, we add to the full valid contact_matrix + # and continue searching. + if (atype_b, atype_a, dist) in contact_matrix: + # generate backbone-backbone exclusions + # perhaps one day can be its own function + excl = Interaction(atoms=(bb_node_A, bb_node_B), + parameters=[], meta={"group": "Go model exclusion"}) + molecule.interactions['exclusions'].append(excl) + symmetrical_matrix.append((atype_a, atype_b, dist)) + else: + contact_matrix.append((atype_a, atype_b, dist)) + else: + if bad_chains_warning == False: + LOGGER.warning("Mismatch between chain IDs in pdb and contact map. This probably means the " + "chain IDs are missing in the pdb and the contact map has all chains = Z.") + bad_chains_warning = True return symmetrical_matrix def compute_go_interaction(self, contacts): diff --git a/vermouth/rcsu/go_vs_includes.py b/vermouth/rcsu/go_vs_includes.py index 0aa622c60..2e2249d92 100644 --- a/vermouth/rcsu/go_vs_includes.py +++ b/vermouth/rcsu/go_vs_includes.py @@ -52,7 +52,7 @@ def run_molecule(self, molecule): if not self.system: raise ValueError('This processor requires a system.') - print(self.backbone, self.atomname) + self.add_virtual_sites(molecule, prefix=moltype, backbone=self.backbone,