From e38ae9cfce0d98f84f26295f1766eda81b91f015 Mon Sep 17 00:00:00 2001 From: zhiyiwu Date: Wed, 29 Jan 2020 16:16:48 +0000 Subject: [PATCH 1/5] Replace distance search with capped_distance --- .../analysis/hbonds/wbridge_analysis.py | 284 ++++++++---------- 1 file changed, 129 insertions(+), 155 deletions(-) diff --git a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py index 9f203f0b596..35ce2d34bca 100644 --- a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py @@ -569,6 +569,7 @@ def analysis(current, output, u, **kwargs): from ..base import AnalysisBase from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch +from MDAnalysis.lib.distances import capped_distance, calc_angles from MDAnalysis import NoDataError, MissingDataWarning, SelectionError from MDAnalysis.lib import distances @@ -789,9 +790,9 @@ def __init__(self, universe, selection1='protein', # set up the donors/acceptors lists if donors is None: - donors = [] + donors = () if acceptors is None: - acceptors = [] + acceptors = () self.forcefield = forcefield self.donors = tuple(set(self.DEFAULT_DONORS[forcefield]).union(donors)) self.acceptors = tuple(set(self.DEFAULT_ACCEPTORS[forcefield]).union(acceptors)) @@ -820,93 +821,90 @@ def _log_parameters(self): acceptor names", self.forcefield) def _update_selection(self): - self._s1 = self.u.select_atoms(self.selection1) - self._s2 = self.u.select_atoms(self.selection2) + self._s1 = self.u.select_atoms(self.selection1).ix + self._s2 = self.u.select_atoms(self.selection2).ix - if self.filter_first and self._s1: + if self.filter_first and len(self._s1): self.logger_debug('Size of selection 1 before filtering:' ' {} atoms'.format(len(self._s1))) - ns_selection_1 = AtomNeighborSearch(self._s1, box=self.box) - self._s1 = ns_selection_1.search(self._s2, self.selection_distance) + ns_selection_1 = AtomNeighborSearch(self.u.atoms[self._s1], box=self.box) + self._s1 = ns_selection_1.search(self.u.atoms[self._s2], self.selection_distance).ix self.logger_debug("Size of selection 1: {0} atoms".format(len(self._s1))) - if not self._s1: + if len(self._s1) == 0: logger.warning('Selection 1 "{0}" did not select any atoms.' .format(str(self.selection1)[:80])) return - if self.filter_first and self._s2: + if self.filter_first and len(self._s2): self.logger_debug('Size of selection 2 before filtering:' ' {} atoms'.format(len(self._s2))) - ns_selection_2 = AtomNeighborSearch(self._s2, box=self.box) - self._s2 = ns_selection_2.search(self._s1, self.selection_distance) + ns_selection_2 = AtomNeighborSearch(self.u.atoms[self._s2], box=self.box) + self._s2 = ns_selection_2.search(self.u.atoms[self._s1], self.selection_distance).ix self.logger_debug('Size of selection 2: {0} atoms'.format(len(self._s2))) - if not self._s2: + if len(self._s2) == 0: logger.warning('Selection 2 "{0}" did not select any atoms.' .format(str(self.selection2)[:80])) return if self.selection1_type in ('donor', 'both'): - self._s1_donors = self._s1.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - self._s1_donors_h = {} - for i, atom in enumerate(self._s1_donors): - tmp = self._get_bonded_hydrogens(atom) - if tmp: - self._s1_donors_h[i] = tmp + self._s1_donors = self.u.atoms[self._s1].select_atoms( + 'name {0}'.format(' '.join(self.donors))).ix + for i, atom_ix in enumerate(self._s1_donors): + tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix + for atom in tmp: + self._s1_donors_h[atom] = atom_ix self.logger_debug("Selection 1 donors: {0}".format(len(self._s1_donors))) self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_donors_h))) if self.selection1_type in ('acceptor', 'both'): - self._s1_acceptors = self._s1.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) + self._s1_acceptors = self.u.atoms[self._s1].select_atoms( + 'name {0}'.format(' '.join(self.acceptors))).ix self.logger_debug("Selection 1 acceptors: {0}".format(len(self._s1_acceptors))) - if not self._s2: + if len(self._s2) == 0: return None if self.selection1_type in ('donor', 'both'): - self._s2_acceptors = self._s2.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) + self._s2_acceptors = self.u.atoms[self._s2].select_atoms( + 'name {0}'.format(' '.join(self.acceptors))).ix self.logger_debug("Selection 2 acceptors: {0:d}".format(len(self._s2_acceptors))) if self.selection1_type in ('acceptor', 'both'): - self._s2_donors = self._s2.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - self._s2_donors_h = {} - for i, d in enumerate(self._s2_donors): - tmp = self._get_bonded_hydrogens(d) - if tmp: - self._s2_donors_h[i] = tmp + self._s2_donors = self.u.atoms[self._s2].select_atoms( + 'name {0}'.format(' '.join(self.donors))).ix + for i, atom_ix in enumerate(self._s2_donors): + tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix + for atom in tmp: + self._s2_donors_h[atom] = atom_ix self.logger_debug("Selection 2 donors: {0:d}".format(len(self._s2_donors))) self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_donors_h))) def _update_water_selection(self): - self._water = self.u.select_atoms(self.water_selection) + self._water = self.u.select_atoms(self.water_selection).ix self.logger_debug('Size of water selection before filtering:' ' {} atoms'.format(len(self._water))) - if self._water and self.filter_first: - filtered_s1 = AtomNeighborSearch(self._water, box=self.box).search(self._s1, + if len(self._water) and self.filter_first: + filtered_s1 = AtomNeighborSearch(self.u.atoms[self._water], box=self.box).search(self.u.atoms[self._s1], self.selection_distance) if filtered_s1: - self._water = AtomNeighborSearch(filtered_s1, box=self.box).search(self._s2, - self.selection_distance) + self._water = AtomNeighborSearch(filtered_s1, box=self.box).search(self.u.atoms[self._s2], + self.selection_distance).ix self.logger_debug("Size of water selection: {0} atoms".format(len(self._water))) - if not self._water: + if len(self._water) == 0: logger.warning("Water selection '{0}' did not select any atoms." .format(str(self.water_selection)[:80])) else: - self._water_donors = self._water.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - self._water_donors_h = {} - for i, d in enumerate(self._water_donors): - tmp = self._get_bonded_hydrogens(d) - if tmp: - self._water_donors_h[i] = tmp + self._water_donors = self.u.atoms[self._water].select_atoms( + 'name {0}'.format(' '.join(self.donors))).ix + for i, atom_ix in enumerate(self._water_donors): + tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix + for atom in tmp: + self._water_donors_h[atom] = atom_ix self.logger_debug("Water donors: {0}".format(len(self._water_donors))) self.logger_debug("Water donor hydrogens: {0}".format(len(self._water_donors_h))) - self._water_acceptors = self._water.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) + self._water_acceptors = self.u.atoms[self._water].select_atoms( + 'name {0}'.format(' '.join(self.acceptors))).ix self.logger_debug("Water acceptors: {0}".format(len(self._water_acceptors))) def _get_bonded_hydrogens(self, atom): @@ -952,23 +950,23 @@ def _prepare(self): # plus order of water bridge times the length of OH bond plus hydrogne bond distance self.selection_distance = (2 * 2 + self.order * (2 + self.distance)) - self._s1_donors = {} + self._s1_donors = [] self._s1_donors_h = {} - self._s1_acceptors = {} + self._s1_acceptors = [] - self._s2_donors = {} + self._s2_donors = [] self._s2_donors_h = {} - self._s2_acceptors = {} + self._s2_acceptors = [] self.box = self.u.dimensions if self.pbc else None self._update_selection() - self._water_donors = {} + self._water_donors = [] self._water_donors_h = {} - self._water_acceptors = {} + self._water_acceptors = [] self.timesteps = [] - if self._s1 and self._s2: + if len(self._s1) and len(self._s2): self._update_water_selection() logger.info("WaterBridgeAnalysis: no atoms found in the selection.") @@ -988,28 +986,47 @@ def _prepare(self): self.start, self.stop, self.step) def _donor2acceptor(self, donors, donor_hs, acceptor): + if len(donors) == 0 or len(acceptor) == 0: + return [] + if self.distance_type != 'heavy': + donors_idx = list(donor_hs.keys()) + else: + donors_idx = donors result = [] - ns_acceptors = AtomNeighborSearch(acceptor, self.box) - for i, donor_h_set in donor_hs.items(): - d = donors[i] - for h in donor_h_set: - res = ns_acceptors.search(h, self.distance) - for a in res: - donor_atom = h if self.distance_type != 'heavy' else d - dist = distances.calc_bonds(donor_atom.position, - a.position, box=self.box) - if dist <= self.distance: - angle = distances.calc_angles(d.position, h.position, - a.position, box=self.box) - angle = np.rad2deg(angle) - if angle >= self.angle: - self.logger_debug( - "D: {0!s} <-> A: {1!s} {2:f} A, {3:f} DEG" \ - .format(h.index, a.index, dist, angle)) - result.append((h.index, d.index, a.index, - (h.resname, h.resid, h.name), - (a.resname, a.resid, a.name), - dist, angle)) + # Code modified from p-j-smith + pairs, distances = capped_distance(self.u.atoms[donors_idx].positions, + self.u.atoms[acceptor].positions, + max_cutoff=self.distance, box=self.box, + return_distances=True) + if self.distance_type != 'heavy': + tmp_donors = [donor_hs[donors_idx[idx]] for idx in pairs[:,0]] + tmp_hydrogens = [donors_idx[idx] for idx in pairs[:,0]] + tmp_acceptors = [acceptor[idx] for idx in pairs[:,1]] + else: + tmp_donors = [] + tmp_hydrogens = [] + tmp_acceptors = [] + for idx in len(pairs[:,0]): + for h in donor_hs[donors[pairs[idx,0]]]: + tmp_donors.append(donors[pairs[idx,0]]) + tmp_hydrogens.append(h) + tmp_acceptors.append(acceptor[pairs[idx,1]]) + + angles = np.rad2deg( + calc_angles( + self.u.atoms[tmp_donors].positions, + self.u.atoms[tmp_hydrogens].positions, + self.u.atoms[tmp_acceptors].positions, + box=self.box + ) + ) + hbond_indices = np.where(angles > self.angle)[0] + for index in hbond_indices: + h = tmp_hydrogens[index] + d = tmp_donors[index] + a = tmp_acceptors[index] + result.append((h, d, a, self._expand_index(h), self._expand_index(a), + distances[index], angles[index])) return result def _single_frame(self): @@ -1018,7 +1035,7 @@ def _single_frame(self): if self.update_selection: self._update_selection() - if self._s1 and self._s2: + if len(self._s1) and len(self._s2): if self.update_water_selection: self._update_water_selection() else: @@ -1032,24 +1049,30 @@ def _single_frame(self): if self.selection1_type in ('donor', 'both'): # check for direct hbond from s1 to s2 - if self._s2_acceptors: - self.logger_debug("Selection 1 Donors <-> Selection 2 Acceptors") - results = self._donor2acceptor(self._s1_donors, self._s1_donors_h, self._s2_acceptors) - for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - water_pool[(a_resname, a_resid)] = None - selection_1.append((h_index, d_index, a_index, None, dist, angle)) - selection_2.append((a_resname, a_resid)) - if self.order > 0 and self._water_acceptors: + self.logger_debug("Selection 1 Donors <-> Selection 2 Acceptors") + results = self._donor2acceptor(self._s1_donors, self._s1_donors_h,self._s2_acceptors) + for line in results: + h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line + water_pool[(a_resname, a_resid)] = None + selection_1.append((h_index, d_index, a_index, None, dist, angle)) + selection_2.append((a_resname, a_resid)) + if self.order > 0: self.logger_debug("Selection 1 Donors <-> Water Acceptors") results = self._donor2acceptor(self._s1_donors, self._s1_donors_h, self._water_acceptors) for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - next_round_water.add((a_resname, a_resid)) + h_index, d_index, a_index, (h_resname, h_resid, h_name), ( + a_resname, a_resid, a_name), dist, angle = line selection_1.append((h_index, d_index, a_index, None, dist, angle)) - if (self.selection1_type in ('acceptor', 'both') and - self._s1_acceptors): + self.logger_debug("Water Donors <-> Selection 2 Acceptors") + results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._s2_acceptors) + for line in results: + h_index, d_index, a_index, (h_resname, h_resid, h_name), ( + a_resname, a_resid, a_name), dist, angle = line + water_pool[(h_resname, h_resid)].append((h_index, d_index, a_index, None, dist, angle)) + selection_2.append((a_resname, a_resid)) + + if self.selection1_type in ('acceptor', 'both'): self.logger_debug("Selection 2 Donors <-> Selection 1 Acceptors") results = self._donor2acceptor(self._s2_donors, self._s2_donors_h, self._s1_acceptors) for line in results: @@ -1059,78 +1082,29 @@ def _single_frame(self): selection_2.append((h_resname, h_resid)) if self.order > 0: + self.logger_debug("Selection 2 Donors <-> Water Acceptors") + results = self._donor2acceptor(self._s2_donors, self._s2_donors_h, self._water_acceptors) + for line in results: + h_index, d_index, a_index, (h_resname, h_resid, h_name), ( + a_resname, a_resid, a_name), dist, angle = line + water_pool[(a_resname, a_resid)].append((a_index, None, h_index, d_index, dist, angle)) + selection_2.append((h_resname, h_resid)) + self.logger_debug("Selection 1 Acceptors <-> Water Donors") results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._s1_acceptors) for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - next_round_water.add((h_resname, h_resid)) + h_index, d_index, a_index, (h_resname, h_resid, h_name), ( + a_resname, a_resid, a_name), dist, angle = line selection_1.append((a_index, None, h_index, d_index, dist, angle)) - for i in range(self.order): - # Narrow down the water selection - selection_resn_id = list(next_round_water) - if (not selection_resn_id) or (not self._water): - logger.warning("No water forming hydrogen bonding with selection 1.") - break - selection_resn_id = ['(resname {} and resid {})'.format( - resname, resid) for resname, resid in selection_resn_id] - water_bridges = self._water.select_atoms(' or '.join(selection_resn_id)) - self.logger_debug("Size of water bridge selection for the {} order of water bridge: {} atoms".format(i+1, - len(water_bridges))) - - water_bridges_donors = water_bridges.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - water_bridges_donors_h = {} - for j, d in enumerate(water_bridges_donors): - tmp = self._get_bonded_hydrogens(d) - if tmp: - water_bridges_donors_h[j] = tmp - self.logger_debug("water bridge donors: {0}".format(len(water_bridges_donors))) - self.logger_debug("water bridge donor hydrogens: {0}".format(len(water_bridges_donors_h))) - water_bridges_acceptors = water_bridges.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) - self.logger_debug("water bridge acceptors: {0}".format(len(water_bridges_acceptors))) - - if i < self.order: - next_round_water = set([]) - # first check if the water can touch the other end - # Finding the hydrogen bonds between water bridge and selection 2 - if self._s2_acceptors: - self.logger_debug("Order {} water donor <-> Selection 2 Acceptors".format(i + 1)) - results = self._donor2acceptor(water_bridges_donors, water_bridges_donors_h, self._s2_acceptors) - for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - water_pool[(h_resname, h_resid)].append((h_index, d_index, a_index, None, dist, angle)) - selection_2.append((a_resname, a_resid)) - - # Finding the hydrogen bonds between water bridge and selection 2 - if water_bridges_acceptors: - self.logger_debug("Selection 2 Donors <-> Order {} water".format(i + 1)) - results = self._donor2acceptor(self._s2_donors, self._s2_donors_h, water_bridges_acceptors) - for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - water_pool[(a_resname, a_resid)].append((a_index, None, h_index, d_index, dist, angle)) - selection_2.append((h_resname, h_resid)) - - # find the water water hydrogen bond - if water_bridges_acceptors: - self.logger_debug("Order {} water acceptor <-> Order {} water donor".format(i + 1, i + 2)) - results = self._donor2acceptor(self._water_donors, self._water_donors_h, water_bridges_acceptors) - for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - water_pool[(a_resname, a_resid)].append((a_index, None, h_index, d_index, dist, angle)) - next_round_water.add((h_resname, h_resid)) - - if water_bridges_donors_h: - self.logger_debug("Order {} water donor <-> Order {} water acceptor".format(i + 1, i + 2)) - results = self._donor2acceptor(water_bridges_donors, water_bridges_donors_h, self._water_acceptors) - for line in results: - h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line - water_pool[(h_resname, h_resid)].append((h_index, d_index, a_index, None, dist, angle)) - next_round_water.add((a_resname, a_resid)) - - # eliminate the water which has been tested - next_round_water = next_round_water.difference(set(water_pool.keys())) + if self.order > 1: + self.logger_debug("Water donor <-> Water Acceptors") + results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._water_acceptors) + for line in results: + h_index, d_index, a_index, (h_resname, h_resid, h_name), ( + a_resname, a_resid, a_name), dist, angle = line + water_pool[(a_resname, a_resid)].append((a_index, None, h_index, d_index, dist, angle)) + water_pool[(h_resname, h_resid)].append((h_index, d_index, a_index, None, dist, angle)) # solve the connectivity network # The following code attempt to generate a water network which is formed by the class dict. From 15bbdfe994031e9d947b6ebe73467e6bf64c2af2 Mon Sep 17 00:00:00 2001 From: zhiyiwu Date: Thu, 30 Jan 2020 10:25:27 +0000 Subject: [PATCH 2/5] Add heavy test --- .../analysis/hbonds/wbridge_analysis.py | 64 ++++++++++--------- .../MDAnalysisTests/analysis/test_wbridge.py | 8 +++ 2 files changed, 43 insertions(+), 29 deletions(-) diff --git a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py index 35ce2d34bca..e9e7a5e0608 100644 --- a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py @@ -821,6 +821,16 @@ def _log_parameters(self): acceptor names", self.forcefield) def _update_selection(self): + self._s1_donors = [] + self._s1_h_donors = {} + self._s1_donors_h = defaultdict(list) + self._s1_acceptors = [] + + self._s2_donors = [] + self._s2_h_donors = {} + self._s2_donors_h = defaultdict(list) + self._s2_acceptors = [] + self._s1 = self.u.select_atoms(self.selection1).ix self._s2 = self.u.select_atoms(self.selection2).ix @@ -854,9 +864,10 @@ def _update_selection(self): for i, atom_ix in enumerate(self._s1_donors): tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix for atom in tmp: - self._s1_donors_h[atom] = atom_ix + self._s1_h_donors[atom] = atom_ix + self._s1_donors_h[atom_ix].append(atom) self.logger_debug("Selection 1 donors: {0}".format(len(self._s1_donors))) - self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_donors_h))) + self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_h_donors))) if self.selection1_type in ('acceptor', 'both'): self._s1_acceptors = self.u.atoms[self._s1].select_atoms( 'name {0}'.format(' '.join(self.acceptors))).ix @@ -874,9 +885,10 @@ def _update_selection(self): for i, atom_ix in enumerate(self._s2_donors): tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix for atom in tmp: - self._s2_donors_h[atom] = atom_ix + self._s2_h_donors[atom] = atom_ix + self._s2_donors_h[atom_ix].append(atom) self.logger_debug("Selection 2 donors: {0:d}".format(len(self._s2_donors))) - self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_donors_h))) + self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_h_donors))) def _update_water_selection(self): self._water = self.u.select_atoms(self.water_selection).ix @@ -900,9 +912,10 @@ def _update_water_selection(self): for i, atom_ix in enumerate(self._water_donors): tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix for atom in tmp: - self._water_donors_h[atom] = atom_ix + self._water_h_donors[atom] = atom_ix + self._water_donors_h[atom_ix].append(atom) self.logger_debug("Water donors: {0}".format(len(self._water_donors))) - self.logger_debug("Water donor hydrogens: {0}".format(len(self._water_donors_h))) + self.logger_debug("Water donor hydrogens: {0}".format(len(self._water_h_donors))) self._water_acceptors = self.u.atoms[self._water].select_atoms( 'name {0}'.format(' '.join(self.acceptors))).ix self.logger_debug("Water acceptors: {0}".format(len(self._water_acceptors))) @@ -950,19 +963,12 @@ def _prepare(self): # plus order of water bridge times the length of OH bond plus hydrogne bond distance self.selection_distance = (2 * 2 + self.order * (2 + self.distance)) - self._s1_donors = [] - self._s1_donors_h = {} - self._s1_acceptors = [] - - self._s2_donors = [] - self._s2_donors_h = {} - self._s2_acceptors = [] - self.box = self.u.dimensions if self.pbc else None self._update_selection() self._water_donors = [] - self._water_donors_h = {} + self._water_h_donors = {} + self._water_donors_h = defaultdict(list) self._water_acceptors = [] self.timesteps = [] @@ -985,13 +991,13 @@ def _prepare(self): logger.info("Starting analysis (frame index start=%d stop=%d, step=%d)", self.start, self.stop, self.step) - def _donor2acceptor(self, donors, donor_hs, acceptor): + def _donor2acceptor(self, donors, h_donors, acceptor): if len(donors) == 0 or len(acceptor) == 0: return [] if self.distance_type != 'heavy': - donors_idx = list(donor_hs.keys()) + donors_idx = list(h_donors.keys()) else: - donors_idx = donors + donors_idx = list(donors.keys()) result = [] # Code modified from p-j-smith pairs, distances = capped_distance(self.u.atoms[donors_idx].positions, @@ -999,16 +1005,16 @@ def _donor2acceptor(self, donors, donor_hs, acceptor): max_cutoff=self.distance, box=self.box, return_distances=True) if self.distance_type != 'heavy': - tmp_donors = [donor_hs[donors_idx[idx]] for idx in pairs[:,0]] + tmp_donors = [h_donors[donors_idx[idx]] for idx in pairs[:,0]] tmp_hydrogens = [donors_idx[idx] for idx in pairs[:,0]] tmp_acceptors = [acceptor[idx] for idx in pairs[:,1]] else: tmp_donors = [] tmp_hydrogens = [] tmp_acceptors = [] - for idx in len(pairs[:,0]): - for h in donor_hs[donors[pairs[idx,0]]]: - tmp_donors.append(donors[pairs[idx,0]]) + for idx in range(len(pairs[:,0])): + for h in donors[donors_idx[pairs[idx,0]]]: + tmp_donors.append(donors_idx[pairs[idx,0]]) tmp_hydrogens.append(h) tmp_acceptors.append(acceptor[pairs[idx,1]]) @@ -1050,7 +1056,7 @@ def _single_frame(self): if self.selection1_type in ('donor', 'both'): # check for direct hbond from s1 to s2 self.logger_debug("Selection 1 Donors <-> Selection 2 Acceptors") - results = self._donor2acceptor(self._s1_donors, self._s1_donors_h,self._s2_acceptors) + results = self._donor2acceptor(self._s1_donors_h, self._s1_h_donors,self._s2_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line water_pool[(a_resname, a_resid)] = None @@ -1058,14 +1064,14 @@ def _single_frame(self): selection_2.append((a_resname, a_resid)) if self.order > 0: self.logger_debug("Selection 1 Donors <-> Water Acceptors") - results = self._donor2acceptor(self._s1_donors, self._s1_donors_h, self._water_acceptors) + results = self._donor2acceptor(self._s1_donors_h, self._s1_h_donors, self._water_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), ( a_resname, a_resid, a_name), dist, angle = line selection_1.append((h_index, d_index, a_index, None, dist, angle)) self.logger_debug("Water Donors <-> Selection 2 Acceptors") - results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._s2_acceptors) + results = self._donor2acceptor(self._water_donors_h, self._water_h_donors, self._s2_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), ( a_resname, a_resid, a_name), dist, angle = line @@ -1074,7 +1080,7 @@ def _single_frame(self): if self.selection1_type in ('acceptor', 'both'): self.logger_debug("Selection 2 Donors <-> Selection 1 Acceptors") - results = self._donor2acceptor(self._s2_donors, self._s2_donors_h, self._s1_acceptors) + results = self._donor2acceptor(self._s2_donors_h, self._s2_h_donors, self._s1_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), (a_resname, a_resid, a_name), dist, angle = line water_pool[(h_resname, h_resid)] = None @@ -1083,7 +1089,7 @@ def _single_frame(self): if self.order > 0: self.logger_debug("Selection 2 Donors <-> Water Acceptors") - results = self._donor2acceptor(self._s2_donors, self._s2_donors_h, self._water_acceptors) + results = self._donor2acceptor(self._s2_donors_h, self._s2_h_donors, self._water_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), ( a_resname, a_resid, a_name), dist, angle = line @@ -1091,7 +1097,7 @@ def _single_frame(self): selection_2.append((h_resname, h_resid)) self.logger_debug("Selection 1 Acceptors <-> Water Donors") - results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._s1_acceptors) + results = self._donor2acceptor(self._water_donors_h, self._water_h_donors, self._s1_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), ( a_resname, a_resid, a_name), dist, angle = line @@ -1099,7 +1105,7 @@ def _single_frame(self): if self.order > 1: self.logger_debug("Water donor <-> Water Acceptors") - results = self._donor2acceptor(self._water_donors, self._water_donors_h, self._water_acceptors) + results = self._donor2acceptor(self._water_donors_h, self._water_h_donors, self._water_acceptors) for line in results: h_index, d_index, a_index, (h_resname, h_resid, h_name), ( a_resname, a_resid, a_name), dist, angle = line diff --git a/testsuite/MDAnalysisTests/analysis/test_wbridge.py b/testsuite/MDAnalysisTests/analysis/test_wbridge.py index c613105b53c..ec2904ac720 100644 --- a/testsuite/MDAnalysisTests/analysis/test_wbridge.py +++ b/testsuite/MDAnalysisTests/analysis/test_wbridge.py @@ -335,6 +335,14 @@ def test_donor_accepter(self, universe_DA): network = wb._network[0] assert_equal(list(network.keys())[0][:4], (1, 0, 2, None)) + def test_donor_accepter_heavy(self, universe_DA): + '''Test zeroth order donor to acceptor hydrogen bonding''' + wb = WaterBridgeAnalysis(universe_DA, 'protein and (resid 1)', + 'protein and (resid 4)', order=0, update_selection=True, debug=True, distance_type='heavy') + wb.run(verbose=False) + network = wb._network[0] + assert_equal(list(network.keys())[0][:4], (1, 0, 2, None)) + def test_donor_accepter_pbc(self, universe_DA_PBC): '''Test zeroth order donor to acceptor hydrogen bonding in PBC conditions''' wb = WaterBridgeAnalysis(universe_DA_PBC, 'protein and (resid 1)', From fb55e7b24f22128009e71361395ef9052c43cf30 Mon Sep 17 00:00:00 2001 From: zhiyiwu Date: Thu, 30 Jan 2020 12:12:36 +0000 Subject: [PATCH 3/5] update water update --- package/MDAnalysis/analysis/hbonds/wbridge_analysis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py index e9e7a5e0608..8a611364896 100644 --- a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py @@ -891,6 +891,11 @@ def _update_selection(self): self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_h_donors))) def _update_water_selection(self): + self._water_donors = [] + self._water_h_donors = {} + self._water_donors_h = defaultdict(list) + self._water_acceptors = [] + self._water = self.u.select_atoms(self.water_selection).ix self.logger_debug('Size of water selection before filtering:' ' {} atoms'.format(len(self._water))) @@ -966,11 +971,6 @@ def _prepare(self): self.box = self.u.dimensions if self.pbc else None self._update_selection() - self._water_donors = [] - self._water_h_donors = {} - self._water_donors_h = defaultdict(list) - self._water_acceptors = [] - self.timesteps = [] if len(self._s1) and len(self._s2): self._update_water_selection() From 21ae01bf2bb63657e52e27069757782e0cd8f775 Mon Sep 17 00:00:00 2001 From: xiki-tempula Date: Sun, 9 Feb 2020 12:50:19 +0000 Subject: [PATCH 4/5] Update CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 2a06b554e6e..63e1fb01bf9 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -65,6 +65,7 @@ Enhancements * Added radius_cut_q as a method to contacts.Contacts (PR #2458) * Added ParmEdParser, ParmEdReader and ParmEdConverter to convert between a parmed.Structure and MDAnalysis Universe (PR #2404) + * Improve the distance search in water bridge analysis with capped_distance (PR #2480) Changes * The fasta2select now always assumes that the gap character in a sequence From 5f4915c961843022ca6e6fe8d7ca3f2ada67a570 Mon Sep 17 00:00:00 2001 From: xiki-tempula Date: Sun, 9 Feb 2020 13:52:10 +0000 Subject: [PATCH 5/5] Update CHANGELOG --- package/CHANGELOG | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 63e1fb01bf9..12674bc2e1a 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,7 +14,7 @@ The rules for this file: ------------------------------------------------------------------------------ mm/dd/yy richardjgowers, kain88-de, lilyminium, p-j-smith, bdice, joaomcteixeira, - PicoCentauri, davidercruz, jbarnoud, RMeli, IAlibay + PicoCentauri, davidercruz, jbarnoud, RMeli, IAlibay, xiki-tempula * 0.21.0 @@ -22,7 +22,7 @@ Fixes * Removes the MassWeight option from gnm (PR #2479). * AtomGroup.guess_bonds now uses periodic boundary information when available (Issue #2350) - * Chainreader and continuous option work correctly when readers work for more + * Chainreader and continuous option work correctly when readers work for more than one format (Issue #2353) * PDBQTParser now gives icode TopologyAttrs (Issue #2361) * GROReader and GROWriter now can handle boxes with box vectors >1000nm @@ -34,9 +34,9 @@ Fixes * The expected frames are made available when a trajectory slice is sliced itself with an incomplete slice and/or with a negative step (Issue #2413) * TXYZ parser uses strings for the atom types like other parsers (Issue #2435) - * ITPParser now parses ITP *and* TOP files from GROMACS, reads #include files, and - substitutes #define variables both from the file and when passed in as a keyword - argument. The directives parsed into bonds, angles, impropers, and dihedrals now + * ITPParser now parses ITP *and* TOP files from GROMACS, reads #include files, and + substitutes #define variables both from the file and when passed in as a keyword + argument. The directives parsed into bonds, angles, impropers, and dihedrals now match TPRParser. (PR #2408) Enhancements @@ -57,13 +57,13 @@ Enhancements * Added match_atoms keyword to analysis.align (Issue #2285, PR #2380) * Added wrap/unwrap transformations (PR #2038) * Read TPR files from Gromacs 2020 (Issue #2412 and #2428) - * Added analysis.align.AverageStructure to get the average structure of an + * Added analysis.align.AverageStructure to get the average structure of an out-of-memory trajectory (Issue #2039) - * Added _add_TopologyObjects, _delete_TopologyObjects, and public convenience + * Added _add_TopologyObjects, _delete_TopologyObjects, and public convenience methods to Universe. Added type and order checking to _Connection topologyattrs. (PR #2382) * Added radius_cut_q as a method to contacts.Contacts (PR #2458) - * Added ParmEdParser, ParmEdReader and ParmEdConverter to + * Added ParmEdParser, ParmEdReader and ParmEdConverter to convert between a parmed.Structure and MDAnalysis Universe (PR #2404) * Improve the distance search in water bridge analysis with capped_distance (PR #2480) @@ -96,7 +96,7 @@ Enhancements * improved atom element guessing in topology.guessers to check for elements after the first element (#2313) * added the zero-based index selection keyword (Issue #1959) - * added position averaging transformation that makes use of the + * added position averaging transformation that makes use of the transformations API (PR #2208) * added find_hydrogen_donors to analysis.bonds.hbond_autocorrel to automatically determine donors where possible (#2181) @@ -130,7 +130,7 @@ Enhancements * added functionality to write files in compressed form(gz,bz2). (Issue #2216, PR #2221) * survival probability additions: residues, intermittency, step with performance, - (PR #2226) + (PR #2226) * added unwrap keyword to center (PR #2275) * added unwrap keyword to center_of_geometry (PR #2279) * added unwrap keyword to center_of_mass (PR #2286) @@ -207,7 +207,7 @@ Deprecations method to be used, this is done automatically in run() Testsuite - * Raised minimum pytest version to 3.3.0, warns now uses the `match` + * Raised minimum pytest version to 3.3.0, warns now uses the `match` argument instead of `match_expr` (Issue #2329) * Add tests for the new water bridge analysis (PR #2087) @@ -230,7 +230,7 @@ Fixes * Fixed reading PDB files with DOS line ending (#2128) Deprecations - * Default ``filename`` directory of align.AlignTraj is deprecated and + * Default ``filename`` directory of align.AlignTraj is deprecated and will change in 1.0 to the current directory. 10/09/18 tylerjereddy, richardjgowers, palnabarun, orbeckst, kain88-de, zemanj, @@ -1802,7 +1802,7 @@ Testsuite MDAnalysisTestData package with a higher release number is available. Testsuite - + * Split off test data trajectories and structures from MDAnalaysis/tests/data into separate package. (Issue 28) * Numbering matches the earliest MDAnalysis release for which the data is